diff --git a/columnphysics/icepack_aerosol.F90 b/columnphysics/icepack_aerosol.F90 index 616065b2..da0715e7 100644 --- a/columnphysics/icepack_aerosol.F90 +++ b/columnphysics/icepack_aerosol.F90 @@ -531,7 +531,7 @@ subroutine update_snow_bgc (dt, nblyr, & end if hslyr_old = hs_old/real(nslyr,kind=dbl_kind) - dzssl = min(hslyr_old/c2, hs_ssl) + dzssl = hslyr_old/c2 dzint = hs_old - dzssl if (aicen > c0) then @@ -546,7 +546,7 @@ subroutine update_snow_bgc (dt, nblyr, & endif hilyr = hi/real(nblyr,kind=dbl_kind) hslyr = hs/real(nslyr,kind=dbl_kind) - dzssl_new = min(hslyr/c2, hs_ssl) + dzssl_new = hslyr/c2 dhs_melts = -melts dhs_snoice = snoice*rhoi/rhos dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice & @@ -738,7 +738,7 @@ subroutine update_snow_bgc (dt, nblyr, & endif hslyr = hs/real(nslyr,kind=dbl_kind) - dzssl_new = min(hslyr/c2, hs_ssl) + dzssl_new = hslyr/c2 dzint_new = max(c0,hs - dzssl_new) if (hs > hs_min) then diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90 index 063a666e..7db6b851 100644 --- a/columnphysics/icepack_algae.F90 +++ b/columnphysics/icepack_algae.F90 @@ -93,6 +93,7 @@ subroutine zbio (dt, nblyr, & n_fed, n_fep, & n_zaero, first_ice, & hice_old, ocean_bio, & + ocean_bio_dh, & bphin, iphin, & iDin, & fswthrul, & @@ -162,7 +163,7 @@ subroutine zbio (dt, nblyr, & flux_bio, & ! total ocean tracer flux (mmol/m^2/s) flux_bion ! category ocean tracer flux (mmol/m^2/s) - real (kind=dbl_kind), intent(inout) :: & + real (kind=dbl_kind), intent(in) :: & hbri_old ! brine height (m) real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & @@ -196,6 +197,10 @@ subroutine zbio (dt, nblyr, & !change to inout when updating ocean fields ocean_bio ! ocean concentrations (mmol/m^3) + real (kind=dbl_kind), dimension (nbtrcr), intent(out) :: & + !change to inout when updating ocean fields + ocean_bio_dh ! ocean concentrations * hbrine * phi (mmol/m^2) + real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & bphin ! Porosity on the bgrid @@ -264,17 +269,15 @@ subroutine zbio (dt, nblyr, & call bgc_carbon_sum(nblyr, hbri_old, trcrn(:), carbonInitial,n_doc,n_dic,n_algae,n_don) - if (write_flux_diag) then - if (aice_old > c0) then - hsnow_i = vsno_old/aice_old - do mm = 1,nbtrcr - call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, & + if (aice_old > c0) then + hsnow_i = vsno_old/aice_old + do mm = 1,nbtrcr + call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, & trcrn(bio_index(mm):bio_index(mm)+nblyr+2), & Tot_BGC_i(mm)) - if (icepack_warnings_aborted(subname)) return - enddo - endif - endif + if (icepack_warnings_aborted(subname)) return + enddo + endif call update_snow_bgc (dt, nblyr, & nslyr, & @@ -301,6 +304,7 @@ subroutine zbio (dt, nblyr, & n_zaero, first_ice, & aicen, vicen, & hice_old, ocean_bio, & + ocean_bio_dh, & flux_bion, bphin, & iphin, trcrn, & iDin, & @@ -847,6 +851,7 @@ subroutine z_biogeochemistry (n_cat, dt, & n_zaero, first_ice, & aicen, vicen, & hice_old, ocean_bio, & + ocean_bio_dh, & flux_bio, bphin, & iphin, trcrn, & iDin, & @@ -908,7 +913,10 @@ subroutine z_biogeochemistry (n_cat, dt, & ocean_bio , & ! ocean concentrations (mmol/m^3) bphin ! Porosity on the bgrid - real (kind=dbl_kind), intent(inout) :: & + real (kind=dbl_kind), dimension (:), intent(out) :: & + ocean_bio_dh ! ocean concentrations * hbrine * phi (mmol/m^2) + + real (kind=dbl_kind), intent(in) :: & hbri_old ! brine height (m) real (kind=dbl_kind), dimension (:,:), intent(out) :: & @@ -1013,7 +1021,7 @@ subroutine z_biogeochemistry (n_cat, dt, & V_alg ! volume of algae (um^3) real (kind=dbl_kind), dimension(nbtrcr) :: & - mobile ! c1 if mobile, c0 otherwise + mobile ! c1 if mobile, c0 otherwise ! local parameters @@ -1048,6 +1056,8 @@ subroutine z_biogeochemistry (n_cat, dt, & dhbot = c0 darcyV = c0 C_top(:) = c0 + C_bot(:) = c0 + ocean_bio_dh(:) = c0 mobile(:) = c0 conserve_C(:) = .true. nitrification(:) = c0 @@ -1060,13 +1070,7 @@ subroutine z_biogeochemistry (n_cat, dt, & iphin_N(k) = iphin(k) bphin_N(1) = bphi_min - if (abs(trcrn(bio_index(m) + k-1)) < accuracy) then - flux_bio(m) = flux_bio(m) + trcrn(bio_index(m) + k-1)* hbri_old * dz(k)/dt - trcrn(bio_index(m) + k-1) = c0 - in_init_cons(k,m) = c0 - else - in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old - endif + in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old if (trcrn(bio_index(m) + k-1) < c0 ) then write(warnstr,*) subname,'zbgc initialization error, first ice = ', first_ice @@ -1170,6 +1174,7 @@ subroutine z_biogeochemistry (n_cat, dt, & endif C_bot(m) = ocean_bio(m)*hbri_old*iphin_N(nblyr+1) + ocean_bio_dh(m) = C_bot(m) enddo ! m @@ -1324,8 +1329,8 @@ subroutine z_biogeochemistry (n_cat, dt, & trcrn(nt_zbgc_frac+mm-1) = zbgc_frac_init(mm) if (sum_tot > c0) trcrn(nt_zbgc_frac+mm-1) = sum_new/sum_tot - if (abs(sum_initial-sum_tot-flux_bio(mm)*dt + source(mm)) > accuracy*max(sum_initial,sum_tot) .or. & - minval(biocons(:)) < c0 .or. minval(initcons_stationary(:)) < c0 & + if ((abs((sum_initial-sum_tot+source(mm))/dt-flux_bio(mm)) > max(accuracy, accuracy*abs(flux_bio(mm)))) & + .or. (minval(biocons(:)) < c0) .or. (minval(initcons_stationary(:)) < c0) & .or. icepack_warnings_aborted()) then write(warnstr,*) subname,'zbgc FCT tracer solution failed, mm:', mm call icepack_warnings_add(warnstr) @@ -1333,9 +1338,9 @@ subroutine z_biogeochemistry (n_cat, dt, & call icepack_warnings_add(warnstr) write(warnstr,*)sum_new,sum_tot,sum_initial,flux_bio(mm),source(mm) call icepack_warnings_add(warnstr) - write(warnstr,*)'error = sum_initial-sum_tot-flux_bio(mm)*dt+source(mm)' + write(warnstr,*)'error = (sum_initial-sum_tot+source(mm))/dt-flux_bio(mm)' call icepack_warnings_add(warnstr) - write(warnstr,*)sum_initial-sum_tot-flux_bio(mm)*dt+source(mm) + write(warnstr,*)(sum_initial-sum_tot+source(mm))/dt-flux_bio(mm) call icepack_warnings_add(warnstr) write(warnstr,*) subname,'sum_new,sum_old:',sum_new,sum_old call icepack_warnings_add(warnstr) @@ -1472,7 +1477,6 @@ subroutine z_biogeochemistry (n_cat, dt, & write(warnstr,*) subname, react(k,m),iphin_N(k),biomat_brine(k,m) call icepack_warnings_add(warnstr) call icepack_warnings_add(subname//' very large bgc value') - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) elseif (bio_tmp < c0) then write(warnstr,*) subname, 'negative bgc' call icepack_warnings_add(warnstr) @@ -2782,7 +2786,8 @@ subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, & C_init_tot , & C_new_tot , & zspace , & !1/nblyr - accuracyC ! centered difference is Order(zspace^2) + accuracyC , & ! centered difference is Order(zspace^2) + var_tmp ! temporary variable character(len=*),parameter :: subname='(check_conservation_FCT)' @@ -2802,8 +2807,8 @@ subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, & C_low(k) = C_new(k) enddo - accuracyC = accuracy*max(c1, C_init_tot, C_new_tot) - fluxbio = (C_init_tot - C_new_tot + source)/dt + accuracyC = 1.0e-11_dbl_kind*max(c1, C_init_tot, C_new_tot) + fluxbio = fluxbio + (C_init_tot - C_new_tot + source)/dt diff_dt =C_new_tot - C_init_tot - (S_top+S_bot+L_bot*C_new(nblyr+1)+L_top*C_new(1))*dt if (minval(C_low) < c0) then @@ -2812,23 +2817,42 @@ subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, & endif if (abs(diff_dt) > accuracyC ) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - write(warnstr,*) subname, 'Conservation of zbgc low order solution failed: diff_dt:',& - diff_dt - write(warnstr,*) subname, 'Total initial tracer', C_init_tot - write(warnstr,*) subname, 'Total final1 tracer', C_new_tot - write(warnstr,*) subname, 'bottom final tracer', C_new(nblyr+1) - write(warnstr,*) subname, 'top final tracer', C_new(1) - write(warnstr,*) subname, 'Near bottom final tracer', C_new(nblyr) - write(warnstr,*) subname, 'Near top final tracer', C_new(2) - write(warnstr,*) subname, 'Top flux*dt into ice:', S_top*dt - write(warnstr,*) subname, 'Bottom flux*dt into ice:', S_bot*dt - write(warnstr,*) subname, 'Remaining bot flux*dt into ice:', L_bot*C_new(nblyr+1)*dt + write(warnstr,*) subname, '' + write(warnstr,*) subname, 'Conservation of zbgc low order solution failed: diff_dt:' + write(warnstr,*) subname, diff_dt + write(warnstr,*) subname, 'Total initial tracer' + write(warnstr,*) subname, C_init_tot + write(warnstr,*) subname, 'Total final1 tracer' + write(warnstr,*) subname, C_new_tot + write(warnstr,*) subname, 'bottom final tracer' + write(warnstr,*) subname, C_new(nblyr+1) + write(warnstr,*) subname, 'top final tracer' + write(warnstr,*) subname, C_new(1) + write(warnstr,*) subname, 'Near bottom final tracer' + write(warnstr,*) subname, C_new(nblyr) + write(warnstr,*) subname, 'Near top final tracer' + write(warnstr,*) subname, C_new(2) + write(warnstr,*) subname, 'Top flux*dt into ice:' + var_tmp = S_top*dt + write(warnstr,*) subname, var_tmp + write(warnstr,*) subname, 'Bottom flux*dt into ice:' + var_tmp = S_bot*dt + write(warnstr,*) subname, var_tmp + write(warnstr,*) subname, 'Remaining bot flux*dt into ice:' + var_tmp = L_bot*C_new(nblyr+1)*dt + write(warnstr,*) subname, var_tmp write(warnstr,*) subname, 'S_bot*dt + L_bot*C_new(nblyr+1)*dt' - write(warnstr,*) subname, S_bot*dt + L_bot*C_new(nblyr+1)*dt - write(warnstr,*) subname, 'fluxbio*dt:', fluxbio*dt - write(warnstr,*) subname, 'fluxbio:', fluxbio - write(warnstr,*) subname, 'Remaining top flux*dt into ice:', L_top*C_new(1)*dt + var_tmp = S_bot*dt + L_bot*C_new(nblyr+1)*dt + write(warnstr,*) subname, var_tmp + write(warnstr,*) subname, 'fluxbio*dt:' + var_tmp = fluxbio*dt + write(warnstr,*) subname, var_tmp + write(warnstr,*) subname, 'fluxbio:' + write(warnstr,*) subname, fluxbio + write(warnstr,*) subname, 'Remaining top flux*dt into ice:' + var_tmp = L_top*C_new(1)*dt + write(warnstr,*) subname, var_tmp + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) endif end subroutine check_conservation_FCT @@ -2869,7 +2893,7 @@ subroutine bgc_column_sum (nblyr, nslyr, hsnow, hbrine, xin, xout) character(len=*),parameter :: subname='(bgc_column_sum)' hslyr = hsnow/real(nslyr,kind=dbl_kind) - dzssl = min(hslyr*p5, hs_ssl) + dzssl = hslyr*p5 dzint = max(c0,hsnow - dzssl) zspace = c1/real(nblyr,kind=dbl_kind) @@ -2891,7 +2915,6 @@ end subroutine bgc_column_sum subroutine bgc_carbon_sum (nblyr, hbrine, xin, xout, n_doc, n_dic, n_algae, n_don) - integer (kind=int_kind), intent(in) :: & nblyr, & ! number of ice layers n_doc, n_dic, n_algae, n_don @@ -2913,6 +2936,8 @@ subroutine bgc_carbon_sum (nblyr, hbrine, xin, xout, n_doc, n_dic, n_algae, n_do integer (kind=int_kind) :: & n, m, iBioCount, iLayer, nBGC ! category/layer index + character(len=*),parameter :: subname='(bgc_carbon_sum)' + zspace(:) = c1/real(nblyr,kind=dbl_kind) zspace(1) = p5*zspace(1) zspace(nblyr+1) = zspace(1) @@ -2994,6 +3019,8 @@ subroutine bgc_carbon_flux (flux_bio_atm, flux_bion, n_doc, & integer (kind=int_kind) :: & m ! biology index + character(len=*),parameter :: subname='(bgc_carbon_flux)' + Tot_Carbon_flux = c0 if (tr_bgc_N) then diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90 index 24986d89..46c7c175 100644 --- a/columnphysics/icepack_brine.F90 +++ b/columnphysics/icepack_brine.F90 @@ -159,7 +159,7 @@ subroutine compute_microS_mushy (nilyr, nblyr, & kperm , & ! average ice permeability (m^2) bphi_min ! surface porosity - real (kind=dbl_kind), intent(inout) :: & + real (kind=dbl_kind), intent(in) :: & hbr_old ! previous timestep brine height (m) real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & @@ -230,7 +230,6 @@ subroutine compute_microS_mushy (nilyr, nblyr, & ! map Sin and qin (cice) profiles to bgc grid surface_S = min_salin - hbr_old = min(hbr_old, maxhbr*hice_old) hinc_old = hice_old call remap_zbgc(nilyr, & diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index b86d1cf7..dc154659 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -1414,8 +1414,8 @@ subroutine zap_snow(dt, nslyr, & endif ! tr_iso if (z_tracers) then - dvssl = min(p5*vsnon/real(nslyr,kind=dbl_kind), hs_ssl*aicen) ! snow surface layer - dvint = vsnon - dvssl ! snow interior + dvssl = p5*vsnon/real(nslyr,kind=dbl_kind) ! snow surface layer + dvint = vsnon - dvssl ! snow interior do it = 1, nbtrcr xtmp = (trcrn(bio_index(it)+nblyr+1)*dvssl + & diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 818faa01..3f015037 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -45,7 +45,7 @@ module icepack_mechred use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice use icepack_tracers, only: nt_apnd, nt_hpnd - use icepack_tracers, only: n_iso + use icepack_tracers, only: n_iso, nblyr, bio_index use icepack_tracers, only: icepack_compute_tracers use icepack_warnings, only: warnstr, icepack_warnings_add @@ -102,7 +102,7 @@ subroutine ridge_ice (dt, ndtd, & mu_rdg, tr_brine, & dardg1dt, dardg2dt, & dvirdgdt, opening, & - fpond, & + fpond, flux_bio, & fresh, fhocn, & faero_ocn, fiso_ocn, & aparticn, krdgn, & @@ -188,6 +188,9 @@ subroutine ridge_ice (dt, ndtd, & real (kind=dbl_kind), dimension(:), intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + flux_bio ! biological and zaerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) @@ -221,6 +224,9 @@ subroutine ridge_ice (dt, ndtd, & real (kind=dbl_kind), dimension (n_aero) :: & maero ! aerosol mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension (nbtrcr) :: & + mbio ! bio mass added to ocean (mmol or kg m-2) + real (kind=dbl_kind), dimension (n_iso) :: & miso ! isotope mass added to ocean (kg m-2) @@ -273,6 +279,7 @@ subroutine ridge_ice (dt, ndtd, & msnow_mlt = c0 esnow_mlt = c0 maero (:) = c0 + mbio (:) = c0 miso (:) = c0 mpond = c0 ardg1 = c0 @@ -401,7 +408,8 @@ subroutine ridge_ice (dt, ndtd, & msnow_mlt, esnow_mlt, & maero, miso, & mpond, Tf, & - aredistn, vredistn) + aredistn, vredistn, & + mbio) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -591,6 +599,11 @@ subroutine ridge_ice (dt, ndtd, & faero_ocn(it) = faero_ocn(it) + maero(it)*dti enddo endif + if (present(flux_bio)) then + do it = 1, nbtrcr + flux_bio(it) = flux_bio(it) + mbio(it)*dti + enddo + endif if (present(fiso_ocn)) then if (tr_iso) then ! check size fiso_ocn vs n_iso ??? @@ -1086,7 +1099,8 @@ subroutine ridge_shift (ntrcr, dt, & msnow_mlt, esnow_mlt, & maero, miso, & mpond, Tf, & - aredistn, vredistn) + aredistn, vredistn, & + mbio) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1163,6 +1177,9 @@ subroutine ridge_shift (ntrcr, dt, & real (kind=dbl_kind), dimension(:), intent(inout) :: & miso ! isotope mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension(:), intent(inout) :: & + mbio ! biology and zaerosol mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension (:), intent(inout), optional :: & aredistn , & ! redistribution function: fraction of new ridge area vredistn ! redistribution function: fraction of new ridge volume @@ -1215,7 +1232,9 @@ subroutine ridge_shift (ntrcr, dt, & hL, hR , & ! left and right limits of integration expL, expR , & ! exponentials involving hL, hR tmpfac , & ! factor by which opening/closing rates are cut - wk1 ! work variable + wk1 , & ! work variable + dzssl , & ! fraction of snow surface biotracers + dzint ! fraction of interior snow biotracers character(len=*),parameter :: subname='(ridge_shift)' @@ -1386,6 +1405,16 @@ subroutine ridge_shift (ntrcr, dt, & enddo endif + if (nbtrcr > 0) then + dzssl = p5/real(nslyr,kind=dbl_kind) + dzint = c1-dzssl + do it = 1, nbtrcr + mbio(it) = mbio(it) + vsrdgn*(c1-fsnowrdg) & + * (trcrn(bio_index(it) + nblyr + 1,n) * dzssl & + + trcrn(bio_index(it) + nblyr + 2,n) * dzint) + enddo + endif + if (tr_pond_topo) then mpond = mpond + ardg1n * trcrn(nt_apnd,n) & * trcrn(nt_hpnd,n) @@ -1875,7 +1904,7 @@ subroutine icepack_step_ridge (dt, ndtd, & mu_rdg, tr_brine, & dardg1dt, dardg2dt, & dvirdgdt, opening, & - fpond, & + fpond, flux_bio, & fresh, fhocn, & faero_ocn, fiso_ocn, & aparticn, krdgn, & diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 19210c88..b5e9a622 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -201,8 +201,9 @@ module icepack_parameters dT_mlt = c1p5 ,&! change in temp for non-melt to melt snow grain ! radius change (C) rsnw_mlt = 1500._dbl_kind,&! maximum melting snow grain radius (10^-6 m) - kalg = 0.60_dbl_kind ! algae absorption coefficient for 0.5 m thick layer + kalg = 0.60_dbl_kind, &! algae absorption coefficient for 0.5 m thick layer ! 0.5 m path of 75 mg Chl a / m2 + R_gC2molC = 12.0107_dbl_kind! g carbon per mol carbon ! weights for albedos ! 4 Jan 2007 BPB Following are appropriate for complete cloud ! in a summer polar atmosphere with 1.5m bare sea ice surface: @@ -561,7 +562,7 @@ subroutine icepack_init_parameters( & phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, & albicev_in, albicei_in, albsnowv_in, & ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, & - kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, & + kalg_in, R_gC2molC_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, & atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, & atmiter_conv_in, calc_dragio_in, & tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, & @@ -571,7 +572,7 @@ subroutine icepack_init_parameters( & bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, & modal_aero_in, use_macromolecules_in, restartbgc_in, skl_bgc_in, & solve_zsal_in, grid_o_in, l_sk_in, & - initbio_frac_in, grid_oS_in, l_skS_in, dEdd_algae_in, & + grid_oS_in, l_skS_in, dEdd_algae_in, & phi_snow_in, T_max_in, fsal_in, & fr_resp_in, algal_vel_in, R_dFe2dust_in, dustFe_sol_in, & op_dep_min_in, fr_graze_s_in, fr_graze_e_in, fr_mort2min_in, & @@ -758,7 +759,8 @@ subroutine icepack_init_parameters( & dT_mlt_in , & ! change in temp for non-melt to melt snow grain ! radius change (C) rsnw_mlt_in , & ! maximum melting snow grain radius (10^-6 m) - kalg_in ! algae absorption coefficient for 0.5 m thick layer + kalg_in , & ! algae absorption coefficient for 0.5 m thick layer + R_gC2molC_in ! g carbon per mol logical (kind=log_kind), intent(in), optional :: & sw_redist_in ! redistribute shortwave @@ -880,7 +882,6 @@ subroutine icepack_init_parameters( & grid_o_in , & ! for bottom flux l_sk_in , & ! characteristic diffusive scale (zsalinity) (m) grid_o_t_in , & ! top grid point length scale - initbio_frac_in, & ! fraction of ocean tracer concentration used to initialize tracer phi_snow_in ! snow porosity at the ice/snow interface real (kind=dbl_kind), intent(in), optional :: & @@ -1178,6 +1179,7 @@ subroutine icepack_init_parameters( & if (present(dT_mlt_in) ) dT_mlt = dT_mlt_in if (present(rsnw_mlt_in) ) rsnw_mlt = rsnw_mlt_in if (present(kalg_in) ) kalg = kalg_in + if (present(R_gC2molC_in) ) R_gC2molC = R_gC2molC_in if (present(kstrength_in) ) kstrength = kstrength_in if (present(krdg_partic_in) ) krdg_partic = krdg_partic_in if (present(krdg_redist_in) ) krdg_redist = krdg_redist_in @@ -1356,7 +1358,6 @@ subroutine icepack_init_parameters( & if (present(grid_o_in) ) grid_o = grid_o_in if (present(l_sk_in) ) l_sk = l_sk_in if (present(grid_o_t_in) ) grid_o_t = grid_o_t_in - if (present(initbio_frac_in) ) initbio_frac = initbio_frac_in if (present(frazil_scav_in) ) frazil_scav = frazil_scav_in if (present(grid_oS_in) ) grid_oS = grid_oS_in if (present(l_skS_in) ) l_skS = l_skS_in @@ -1527,7 +1528,7 @@ subroutine icepack_query_parameters( & albedo_type_out, albicev_out, albicei_out, albsnowv_out, & albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, & rsnw_mlt_out, dEdd_algae_out, & - kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, & + kalg_out, R_gC2molC_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, & atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, & atmiter_conv_out, calc_dragio_out, & tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, & @@ -1732,7 +1733,8 @@ subroutine icepack_query_parameters( & dT_mlt_out , & ! change in temp for non-melt to melt snow grain ! radius change (C) rsnw_mlt_out , & ! maximum melting snow grain radius (10^-6 m) - kalg_out ! algae absorption coefficient for 0.5 m thick layer + kalg_out , & ! algae absorption coefficient for 0.5 m thick layer + R_gC2molC_out ! grams carbon per mol logical (kind=log_kind), intent(out), optional :: & sw_redist_out ! redistribute shortwave @@ -2184,6 +2186,7 @@ subroutine icepack_query_parameters( & if (present(dT_mlt_out) ) dT_mlt_out = dT_mlt if (present(rsnw_mlt_out) ) rsnw_mlt_out = rsnw_mlt if (present(kalg_out) ) kalg_out = kalg + if (present(R_gC2molC_out) ) R_gC2molC_out = R_gC2molC if (present(kstrength_out) ) kstrength_out = kstrength if (present(krdg_partic_out) ) krdg_partic_out = krdg_partic if (present(krdg_redist_out) ) krdg_redist_out = krdg_redist @@ -2494,6 +2497,7 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " dT_mlt = ", dT_mlt write(iounit,*) " rsnw_mlt = ", rsnw_mlt write(iounit,*) " kalg = ", kalg + write(iounit,*) " R_gC2molC = ", R_gC2molC write(iounit,*) " kstrength = ", kstrength write(iounit,*) " krdg_partic= ", krdg_partic write(iounit,*) " krdg_redist= ", krdg_redist diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index b6048dbc..5267af04 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -23,7 +23,7 @@ module icepack_therm_itd use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10 use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity - use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss + use icepack_parameters, only: phi_init, dsin0_frazil, salt_loss use icepack_parameters, only: Tliquidus_max use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist use icepack_parameters, only: kitd, ktherm @@ -31,7 +31,7 @@ module icepack_therm_itd use icepack_parameters, only: cpl_frazil, update_ocn_f, saltflux_option use icepack_parameters, only: icepack_chkoptargflag - use icepack_tracers, only: ntrcr, nbtrcr + use icepack_tracers, only: ntrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice @@ -969,6 +969,7 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), dimension (ncat) :: & aicen_init, & ! initial area fraction vicen_init, & ! volume per unit area of ice (m) + vsnon_init, & ! initial volume of snow (m) G_radialn , & ! rate of lateral melt (m/s) delta_an , & ! change in the ITD rsiden ! delta_an/aicen @@ -1000,7 +1001,7 @@ subroutine lateral_melt (dt, ncat, & dvint = c0 bin1_arealoss = c0 tmp = c0 - vicen_init = c0 + vicen_init = vicen(:) G_radialn = c0 delta_an = c0 rsiden = c0 @@ -1119,6 +1120,7 @@ subroutine lateral_melt (dt, ncat, & ! state variables vicen_init(n) = vicen(n) + vsnon_init(n) = vsnon(n) aicen(n) = aicen(n) * (c1 - rsiden(n)) vicen(n) = vicen(n) * (c1 - rsiden(n)) vsnon(n) = vsnon(n) * (c1 - rsiden(n)) @@ -1228,8 +1230,8 @@ subroutine lateral_melt (dt, ncat, & !----------------------------------------------------------------- if (z_tracers) then ! snow tracers - dvssl = min(p5*vsnon(n)/real(nslyr,kind=dbl_kind), hs_ssl*aicen(n)) ! snow surface layer - dvint = vsnon(n) - dvssl ! snow interior + dvssl = p5*vsnon_init(n)/real(nslyr,kind=dbl_kind) ! snow surface layer + dvint = vsnon_init(n) - dvssl ! snow interior do k = 1, nbtrcr flux_bio(k) = flux_bio(k) & + (trcrn(bio_index(k)+nblyr+1,n)*dvssl & @@ -1243,9 +1245,9 @@ subroutine lateral_melt (dt, ncat, & if (z_tracers) & call lateral_melt_bgc(dt, & ncat, nblyr, & - rside, vicen_init, & !echmod: use rsiden - trcrn, & - flux_bio, nbtrcr) + rsiden, vicen_init, & + trcrn, flux_bio, & + nbtrcr) if (icepack_warnings_aborted(subname)) return endif ! flag @@ -1301,10 +1303,10 @@ end subroutine lateral_melt subroutine add_new_ice (ncat, nilyr, & nfsd, nblyr, & n_aero, dt, & - ntrcr, nltrcr, & + ntrcr, & hin_max, ktherm, & aicen, trcrn, & - vicen, vsnon1, & + vicen, & aice0, aice, & frzmlt, frazil, & frz_onset, yday, & @@ -1312,7 +1314,8 @@ subroutine add_new_ice (ncat, nilyr, & Tf, sss, & salinz, phi_init, & dSin0_frazil, & - bgrid, cgrid, igrid, & + bgrid, cgrid, & + igrid, & nbtrcr, flux_bio, & ocean_bio, & frazil_diag, & @@ -1334,7 +1337,7 @@ subroutine add_new_ice (ncat, nilyr, & nilyr , & ! number of ice layers nblyr , & ! number of bio layers ntrcr , & ! number of tracers - nltrcr, & ! number of zbgc tracers + nbtrcr, & ! number of zbgc tracers n_aero, & ! number of aerosol tracers ktherm ! type of thermodynamics (-1 none, 1 BL99, 2 mushy) @@ -1349,8 +1352,7 @@ subroutine add_new_ice (ncat, nilyr, & aice , & ! total concentration of ice frzmlt, & ! freezing/melting potential (W/m^2) Tf , & ! freezing temperature (C) - sss , & ! sea surface salinity (ppt) - vsnon1 ! category 1 snow volume per ice area (m) + sss ! sea surface salinity (ppt) real (kind=dbl_kind), dimension (:), intent(inout) :: & aicen , & ! concentration of ice @@ -1390,9 +1392,6 @@ subroutine add_new_ice (ncat, nilyr, & real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & cgrid ! CICE vertical coordinate - integer (kind=int_kind), intent(in) :: & - nbtrcr ! number of biology tracers - real (kind=dbl_kind), dimension (:), intent(inout) :: & flux_bio ! tracer flux to ocean from biology (mmol/m^2/s) @@ -1810,7 +1809,6 @@ subroutine add_new_ice (ncat, nilyr, & ncats = 1 ! add new ice to category 1 by default if (tr_fsd) ncats = ncat ! add new ice laterally to all categories - do n = 1, ncats if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n @@ -1879,7 +1877,7 @@ subroutine add_new_ice (ncat, nilyr, & fiso_ocn(it) = fiso_ocn(it) & - frazil_conc*rhoi*vi0new/dt enddo - endif + endif ! if iso if (tr_lvl) then alvl = trcrn(nt_alvl,n) @@ -1898,7 +1896,7 @@ subroutine add_new_ice (ncat, nilyr, & trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n)) endif endif - endif + endif ! vicen > 0 do k = 1, nilyr if (vicen(n) > c0) then @@ -1910,11 +1908,9 @@ subroutine add_new_ice (ncat, nilyr, & trcrn(nt_sice+k-1,n) = & (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n) endif - enddo - - endif ! vi0new > 0 - - enddo ! ncats + enddo ! nilyr + endif ! vin0new > c0 + enddo ! n if (conserv_check) then @@ -1947,14 +1943,13 @@ subroutine add_new_ice (ncat, nilyr, & ! Biogeochemistry !----------------------------------------------------------------- if (tr_brine .or. nbtrcr > 0) then - call add_new_ice_bgc(dt, nblyr, & - ncat, nilyr, nltrcr, & + call add_new_ice_bgc(dt, nblyr, ncats, & + ncat, nilyr, & bgrid, cgrid, igrid, & aicen_init, vicen_init, vi0_init, & - aicen, vicen, vsnon1, & - vi0new, ntrcr, trcrn, & - nbtrcr, sss, ocean_bio,& - flux_bio, hsurp) + aicen, vicen, vin0new, & + ntrcr, trcrn, nbtrcr, & + ocean_bio, flux_bio, hsurp) if (icepack_warnings_aborted(subname)) return endif @@ -1972,10 +1967,10 @@ end subroutine add_new_ice ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL - subroutine icepack_step_therm2 (dt, ncat, nltrcr, & + subroutine icepack_step_therm2 (dt, ncat, & nilyr, nslyr, & hin_max, nblyr, & - aicen, & + aicen, nbtrcr, & vicen, vsnon, & aicen_init, vicen_init, & trcrn, & @@ -2011,8 +2006,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories - nltrcr , & ! number of zbgc tracers nblyr , & ! number of bio layers + nbtrcr , & ! number of bio tracers nilyr , & ! number of ice layers nslyr ! number of snow layers @@ -2190,6 +2185,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & ! Compute fractional ice area in each grid cell. !----------------------------------------------------------------- + flux_bio(:) = c0 + call aggregate_area (ncat, aicen, aice, aice0) if (icepack_warnings_aborted(subname)) return @@ -2233,10 +2230,10 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & call add_new_ice (ncat, nilyr, & nfsd, nblyr, & n_aero, dt, & - ntrcr, nltrcr, & + ntrcr, & hin_max, ktherm, & aicen, trcrn, & - vicen, vsnon(1), & + vicen, & aice0, aice, & frzmlt, frazil, & frz_onset, yday, & diff --git a/columnphysics/icepack_zbgc.F90 b/columnphysics/icepack_zbgc.F90 index 9692b04f..5503469d 100644 --- a/columnphysics/icepack_zbgc.F90 +++ b/columnphysics/icepack_zbgc.F90 @@ -86,23 +86,21 @@ module icepack_zbgc ! Adjust biogeochemical tracers when new frazil ice forms - subroutine add_new_ice_bgc (dt, nblyr, & - ncat, nilyr, nltrcr, & + subroutine add_new_ice_bgc (dt, nblyr, ncats, & + ncat, nilyr, & bgrid, cgrid, igrid, & aicen_init, vicen_init, vi0_init, & - aicen, vicen, vsnon1, & - vi0new, & + aicen, vicen, vin0new, & ntrcr, trcrn, nbtrcr, & - sss, ocean_bio, flux_bio, & - hsurp) + ocean_bio, flux_bio, hsurp) integer (kind=int_kind), intent(in) :: & nblyr , & ! number of bio layers ncat , & ! number of thickness categories nilyr , & ! number of ice layers - nltrcr , & ! number of zbgc tracers nbtrcr , & ! number of biology tracers - ntrcr ! number of tracers in use + ntrcr , & ! number of tracers in use + ncats ! 1 without floe size distribution or ncat real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & bgrid ! biology nondimensional vertical grid points @@ -118,22 +116,18 @@ subroutine add_new_ice_bgc (dt, nblyr, & real (kind=dbl_kind), dimension (:), intent(in) :: & aicen_init , & ! initial concentration of ice - vicen_init , & ! intiial volume per unit area of ice (m) + vicen_init , & ! initial volume per unit area of ice (m) aicen , & ! concentration of ice vicen ! volume per unit area of ice (m) - real (kind=dbl_kind), intent(in) :: & - vsnon1 ! category 1 snow volume per unit area (m) - real (kind=dbl_kind), dimension (:,:), intent(inout) :: & trcrn ! ice tracers - real (kind=dbl_kind), intent(in) :: & - sss !sea surface salinity (ppt) + real (kind=dbl_kind), dimension (:), intent(in) :: & + vin0new ! ice tracers real (kind=dbl_kind), intent(in) :: & - vi0_init , & ! volume of new ice added to cat 1 (intial) - vi0new ! volume of new ice added to cat 1 + vi0_init ! volume of new ice added to cat 1 (intial) real (kind=dbl_kind), intent(in) :: & hsurp ! thickness of new ice added to each cat @@ -148,8 +142,10 @@ subroutine add_new_ice_bgc (dt, nblyr, & integer (kind=int_kind) :: & location , & ! 1 (add frazil to bottom), 0 (add frazil throughout) + m , & ! bio index n , & ! ice category index - k ! ice layer index + k , & ! ice layer index + nbiolayer real (kind=dbl_kind) :: & vbri1 , & ! starting volume of existing brine @@ -164,7 +160,8 @@ subroutine add_new_ice_bgc (dt, nblyr, & vbrin ! trcrn(nt_fbri,n)*vicen(n) real (kind=dbl_kind) :: & - vice_new ! vicen_init + vsurp + vice_new , & ! vicen_init + vsurp + bio0new ! ocean_bio * zbgc_init_fac real (kind=dbl_kind) :: & Tmlts ! melting temperature (oC) @@ -174,9 +171,15 @@ subroutine add_new_ice_bgc (dt, nblyr, & character(len=*),parameter :: subname='(add_new_ice_bgc)' + real (kind=dbl_kind), dimension (nblyr+1) :: & + zspace ! vertical grid spacing !----------------------------------------------------------------- ! brine !----------------------------------------------------------------- + + zspace(:) = c1/real(nblyr,kind=dbl_kind) + zspace(1) = p5*zspace(2) + zspace(nblyr+1) = zspace(1) vbrin(:) = c0 do n = 1, ncat vbrin(n) = vicen_init(n) @@ -211,67 +214,57 @@ subroutine add_new_ice_bgc (dt, nblyr, & vsurp = hsurp * aicen_init(n) vbrin(n) = vbrin(n) + vsurp vice_new = vicen_init(n) + vsurp - if (tr_brine .and. vicen(n) > c0) then - trcrn(nt_fbri,n) = vbrin(n)/vicen(n) - elseif (tr_brine .and. vicen(n) <= c0) then + + if (tr_brine .and. vice_new > c0) then + trcrn(nt_fbri,n) = vbrin(n)/vice_new + elseif (tr_brine .and. vice_new <= c0) then trcrn(nt_fbri,n) = c1 endif - if (nltrcr > 0) then - location = 1 - call adjust_tracer_profile(nbtrcr, dt, ntrcr, & - aicen_init(n), & - vbrin(n), & - vice_new, & - trcrn(:,n), & - vtmp, & - vsurp, sss, & - nilyr, nblyr, & - bgrid, & - cgrid, & - ocean_bio, igrid, & - location) + if (nbtrcr > 0) then + do m = 1, nbtrcr + bio0new = ocean_bio(m)*zbgc_init_frac(m) + nbiolayer = nblyr+1 + call update_vertical_bio_tracers(nbiolayer, trcrn(bio_index(m):bio_index(m) + nblyr,n), & + vtmp, vbrin(n), bio0new,zspace(:)) + enddo !nbtrcr if (icepack_warnings_aborted(subname)) return - endif ! nltrcr + endif ! nbtrcr endif ! hsurp > 0 enddo ! n !----------------------------------------------------------------- ! Combine bgc in new ice grown in open water with category 1 ice. !----------------------------------------------------------------- + do n = 1, ncats + if (vin0new(n) > c0) then - if (vi0new > c0) then - - vbri1 = vbrin(1) - vbrin(1) = vbrin(1) + vi0new - if (tr_brine .and. vicen(1) > c0) then - trcrn(nt_fbri,1) = vbrin(1)/vicen(1) - elseif (tr_brine .and. vicen(1) <= c0) then - trcrn(nt_fbri,1) = c1 + vbri1 = vbrin(n) + vbrin(n) = vbrin(n) + vin0new(n) + if (tr_brine .and. vicen(n) > c0) then + trcrn(nt_fbri,n) = vbrin(n)/vicen(n) + elseif (tr_brine .and. vicen(n) <= c0) then + trcrn(nt_fbri,n) = c1 endif ! Diffuse_bio handles concentration changes from ice growth/melt ! ice area changes ! add salt throughout, location = 0 - if (nltrcr > 0) then - location = 0 - call adjust_tracer_profile(nbtrcr, dt, ntrcr, & - aicen(1), & - vbrin(1), & - vicen(1), & - trcrn(:,1), & - vbri1, & - vi0new, sss, & - nilyr, nblyr, & - bgrid, & - cgrid, & - ocean_bio, igrid, & - location) + if (nbtrcr > 0 .and. vbrin(n) > c0) then + do m = 1, nbtrcr + bio0new = ocean_bio(m)*zbgc_init_frac(m) + do k = 1, nblyr+1 + trcrn(bio_index(m) + k-1,n) = & + (trcrn(bio_index(m) + k-1,n)*vbri1 + bio0new * vin0new(n))/vbrin(n) + enddo + enddo + if (icepack_warnings_aborted(subname)) return - endif ! nltrcr > 0 - endif ! vi0new > 0 + endif ! nbtrcr > 0 + endif ! vin0new(n) > 0 + enddo ! n = 1,ncats if (tr_brine .and. conserv_check) then call column_sum (ncat, vbrin, vbri_final) @@ -291,11 +284,11 @@ end subroutine add_new_ice_bgc ! When sea ice melts laterally, flux bgc to ocean - subroutine lateral_melt_bgc (dt, & - ncat, nblyr, & - rside, vicen, & - trcrn, & - flux_bio, nbltrcr) + subroutine lateral_melt_bgc (dt, & + ncat, nblyr, & + rsiden, vicen_init,& + trcrn, flux_bio, & + nbltrcr) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -306,37 +299,43 @@ subroutine lateral_melt_bgc (dt, & dt ! time step (s) real (kind=dbl_kind), dimension(:), intent(in) :: & - vicen ! volume per unit area of ice (m) + vicen_init! volume per unit area of ice (m) real (kind=dbl_kind), dimension (:,:), intent(in) :: & trcrn ! tracer array - real (kind=dbl_kind), intent(in) :: & - rside ! fraction of ice that melts laterally + real (kind=dbl_kind), dimension(:), intent(in) :: & + rsiden ! fraction of ice that melts laterally real (kind=dbl_kind), dimension(:), intent(inout) :: & flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s) ! local variables + real (kind=dbl_kind) :: & + total_bio_initial, & ! initial column tracer concentration (mmol/m2) + total_bio_final ! final column tracer concentration (mmol/m20 + integer (kind=int_kind) :: & k , & ! layer index m , & ! n ! category index - real (kind=dbl_kind) :: & - zspace ! bio grid spacing + real (kind=dbl_kind), dimension (nblyr+1) :: & + zspace ! vertical grid spacing character(len=*),parameter :: subname='(lateral_melt_bgc)' - zspace = c1/(real(nblyr,kind=dbl_kind)) + zspace(:) = c1/real(nblyr,kind=dbl_kind) + zspace(1) = p5*zspace(2) + zspace(nblyr+1) = zspace(1) do m = 1, nbltrcr do n = 1, ncat do k = 1, nblyr+1 flux_bio(m) = flux_bio(m) + trcrn(nt_fbri,n) & - * vicen(n)*zspace*trcrn(bio_index(m)+k-1,n) & - * rside/dt + * vicen_init(n)*zspace(k)*trcrn(bio_index(m)+k-1,n) & + * rsiden(n)/dt enddo enddo enddo @@ -349,117 +348,107 @@ end subroutine lateral_melt_bgc ! ! author: Nicole Jeffery, LANL - subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & - aicen, vbrin, & - vicen, trcrn, & - vtmp, & - vsurp, sss, & - nilyr, nblyr, & - bgrid, & - cgrid, ocean_bio, & - igrid, location) +! subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & +! aicen, vbrin, & +! vicen, trcrn, & +! vtmp, & +! vsurp, & +! nilyr, nblyr, & +! bgrid, & +! cgrid, ocean_bio, & +! igrid, location)!!! - integer (kind=int_kind), intent(in) :: & - location , & ! 1 (add frazil to bottom), 0 (add frazil throughout) - ntrcr , & ! number of tracers in use - nilyr , & ! number of ice layers - nbtrcr , & ! number of biology tracers - nblyr ! number of biology layers +! integer (kind=int_kind), intent(in) :: & +! location , & ! 1 (add frazil to bottom), 0 (add frazil throughout) +! ntrcr , & ! number of tracers in use +! nilyr , & ! number of ice layers +! nbtrcr , & ! number of biology tracers +! nblyr ! number of biology layers - real (kind=dbl_kind), intent(in) :: & - dt ! timestep (s) +! real (kind=dbl_kind), intent(in) :: & +! dt ! timestep (s)! - real (kind=dbl_kind), intent(in) :: & - aicen , & ! concentration of ice - vicen , & ! volume of ice - sss , & ! ocean salinity (ppt) - ! hsurp , & ! flags new ice added to each cat - vsurp , & ! volume of new ice added to each cat - vtmp ! total volume of new and old ice +! real (kind=dbl_kind), intent(in) :: & +! aicen , & ! concentration of ice +! vicen , & ! volume of ice +! vsurp , & ! volume of new ice added to each cat +! vtmp ! total volume of new and old ice - real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: & - ocean_bio +! real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: & +! ocean_bio - real (kind=dbl_kind), intent(in) :: & - vbrin ! fbri*volume per unit area of ice (m) +! real (kind=dbl_kind), intent(in) :: & +! vbrin ! fbri*volume per unit area of ice (m) - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - igrid ! zbio grid +! real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & +! igrid ! zbio grid - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - bgrid ! zsal grid +! real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & +! bgrid ! zsal grid - real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & - cgrid ! CICE grid +! real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & +! cgrid ! CICE grid - real (kind=dbl_kind), dimension (ntrcr), intent(inout) :: & - trcrn ! ice tracers +! real (kind=dbl_kind), dimension (ntrcr), intent(inout) :: & +! trcrn ! ice tracers ! local variables - real (kind=dbl_kind), dimension (ntrcr+2) :: & - trtmp0, & ! temporary, remapped tracers - trtmp ! temporary, remapped tracers +! real (kind=dbl_kind), dimension (ntrcr+2) :: & +! trtmp0, & ! temporary, remapped tracers +! trtmp ! temporary, remapped tracers - real (kind=dbl_kind) :: & - hin , & ! ice height - hinS_new, & ! brine height - temp_S +! integer (kind=int_kind) :: & +! k, m - integer (kind=int_kind) :: & - k, m +! real (kind=dbl_kind), dimension (nblyr+1) :: & +! C_stationary ! stationary bulk concentration*h (mmol/m^2) - real (kind=dbl_kind), dimension (nblyr+1) :: & - C_stationary ! stationary bulk concentration*h (mmol/m^2) +! real(kind=dbl_kind) :: & +! top_conc , & ! salinity or bgc ocean concentration of frazil +! fluxb , & ! needed for regrid (set to zero here) +! hbri_old , & ! previous timestep brine height +! hbri ! brine height - real (kind=dbl_kind), dimension (nblyr) :: & - S_stationary ! stationary bulk concentration*h (ppt*m) +! character(len=*),parameter :: subname='(adjust_tracer_profile)' - real(kind=dbl_kind) :: & - top_conc , & ! salinity or bgc ocean concentration of frazil - fluxb , & ! needed for regrid (set to zero here) - hbri_old , & ! previous timestep brine height - hbri ! brine height +! trtmp0(:) = c0 +! trtmp(:) = c0 +! fluxb = c0 - character(len=*),parameter :: subname='(adjust_tracer_profile)' + ! if (location == 1 .and. vbrin > c0) then ! add frazil to bottom - trtmp0(:) = c0 - trtmp(:) = c0 - fluxb = c0 +! hbri = vbrin +! hbri_old = vtmp - if (location == 1 .and. vbrin > c0) then ! add frazil to bottom +! do m = 1, nbtrcr +! top_conc = ocean_bio(m)*zbgc_init_frac(m) +! do k = 1, nblyr+1 +! C_stationary(k) = trcrn(bio_index(m) + k-1)* hbri_old +! enddo !k +! call regrid_stationary (C_stationary, hbri_old, & +! hbri, dt, & +! ntrcr, & +! nblyr, top_conc, & +! igrid, fluxb ) +! if (icepack_warnings_aborted(subname)) return +! do k = 1, nblyr+1 +! trcrn(bio_index(m) + k-1) = C_stationary(k)/hbri +! enddo !k +! enddo !m - hbri = vbrin - hbri_old = vtmp +! elseif (vbrin > c0) then ! add frazil throughout location == 0 .and. - do m = 1, nbtrcr - top_conc = ocean_bio(m)*zbgc_init_frac(m) - do k = 1, nblyr+1 - C_stationary(k) = trcrn(bio_index(m) + k-1)* hbri_old - enddo !k - call regrid_stationary (C_stationary, hbri_old, & - hbri, dt, & - ntrcr, & - nblyr, top_conc, & - igrid, fluxb ) - if (icepack_warnings_aborted(subname)) return - do k = 1, nblyr+1 - trcrn(bio_index(m) + k-1) = C_stationary(k)/hbri - enddo !k - enddo !m +! do k = 1, nblyr+1 + ! do m = 1, nbtrcr +! trcrn(bio_index(m) + k-1) = (trcrn(bio_index(m) + k-1) * vtmp & +! + ocean_bio(m)*zbgc_init_frac(m) * vsurp) / vbrin +! enddo +! enddo - elseif (vbrin > c0) then ! add frazil throughout location == 0 .and. +! endif ! location - do k = 1, nblyr+1 - do m = 1, nbtrcr - trcrn(bio_index(m) + k-1) = (trcrn(bio_index(m) + k-1) * vtmp & - + ocean_bio(m)*zbgc_init_frac(m) * vsurp) / vbrin - enddo - enddo - - endif ! location - - end subroutine adjust_tracer_profile +! end subroutine adjust_tracer_profile !======================================================================= !autodocument_start icepack_init_bgc @@ -488,7 +477,7 @@ subroutine icepack_init_bgc(ncat, nblyr, nilyr, & real (kind=dbl_kind), intent(in) :: & sss ! sea surface salinity (ppt) - real (kind=dbl_kind), dimension (:), intent(inout) :: & + real (kind=dbl_kind), dimension (:), intent(in) :: & ocean_bio_all ! fixed order, all values even for tracers false real (kind=dbl_kind), dimension (:), intent(out) :: & @@ -1626,6 +1615,8 @@ subroutine icepack_init_bgc_trcr(nk, nt_fbri, & trcr_base2, & trcr_base3 + character(len=*),parameter :: subname='(icepack_init_bgc_trcr)' + nt_bgc = ntrcr + 1 nbtrcr = nbtrcr + 1 nlt_bgc = nbtrcr @@ -1682,7 +1673,7 @@ subroutine icepack_biogeochemistry(dt, & upNO, upNH, iDi, iki, zfswin, & darcy_V, grow_net, & PP_net, hbri,dhbr_bot, dhbr_top, Zoo,& - fbio_snoice, fbio_atmice, ocean_bio, & + fbio_snoice, fbio_atmice, ocean_bio_dh, ocean_bio, & first_ice, fswpenln, bphi, bTiz, ice_bio_net, & snow_bio_net, totalChla, fswthrun, & bgrid, igrid, icgrid, cgrid, & @@ -1709,7 +1700,6 @@ subroutine icepack_biogeochemistry(dt, & igrid , & ! biology vertical interface points cgrid , & ! CICE vertical coordinate icgrid , & ! interface grid for CICE (shortwave variable) - ocean_bio , & ! contains all the ocean bgc tracer concentrations fbio_snoice , & ! fluxes from snow to ice fbio_atmice , & ! fluxes from atm to ice dhbr_top , & ! brine top change @@ -1718,7 +1708,13 @@ subroutine icepack_biogeochemistry(dt, & hin_old , & ! old ice thickness ice_bio_net , & ! depth integrated tracer (mmol/m^2) snow_bio_net , & ! depth integrated snow tracer (mmol/m^2) - flux_bio ! all bio fluxes to ocean + flux_bio ! all bio fluxes to ocean + + real (kind=dbl_kind), dimension (:), intent(in) :: & + ocean_bio ! contains the ocean bgc tracer concentrations in use (mmol/m^3) + + real (kind=dbl_kind), dimension (:), intent(out) :: & + ocean_bio_dh ! The ocean bgc tracer concentrations in use * brine thickness * phi (mmol/m^2) logical (kind=log_kind), dimension (:), intent(inout) :: & first_ice ! distinguishes ice that disappears (e.g. melts) @@ -1783,7 +1779,8 @@ subroutine icepack_biogeochemistry(dt, & ! local variables integer (kind=int_kind) :: & - n, mm ! thickness category index + k , & ! vertical index + n, mm ! thickness category index real (kind=dbl_kind) :: & hin , & ! new ice thickness @@ -1824,8 +1821,8 @@ subroutine icepack_biogeochemistry(dt, & character(len=*),parameter :: subname='(icepack_biogeochemistry)' zspace(:) = c1/real(nblyr,kind=dbl_kind) - zspace(1) = p5*zspace(1) - zspace(nblyr+1) = p5*zspace(nblyr+1) + zspace(1) = p5*zspace(2) + zspace(nblyr+1) = zspace(1) bioPorosityIceCell(:) = c0 bioSalinityIceCell(:) = c0 @@ -1836,6 +1833,7 @@ subroutine icepack_biogeochemistry(dt, & !----------------------------------------------------------------- ! initialize !----------------------------------------------------------------- + flux_bion(:,n) = c0 hin_old(n) = c0 if (aicen_init(n) > puny) then hin_old(n) = vicen_init(n) & @@ -1943,6 +1941,7 @@ subroutine icepack_biogeochemistry(dt, & n_fed, n_fep, & n_zaero, first_ice(n), & hin_old(n), ocean_bio(1:nbtrcr), & + ocean_bio_dh(1:nbtrcr), & bphi(:,n), iphin, & iDi(:,n), & fswpenln(:,n), & @@ -1988,6 +1987,16 @@ subroutine icepack_biogeochemistry(dt, & first_ice(n) = .false. + else + do mm = 1, nbtrcr + do k = 1, nblyr+1 + flux_bion(mm,n) = flux_bion(mm,n) + trcrn(bio_index(mm) + k-1,n) * & + hin_old(n) * zspace(k)/dt * trcrn(nt_fbri,n) + flux_bio(mm) = flux_bio(mm) + trcrn(bio_index(mm) + k-1,n) * & + vicen_init(n) * zspace(k)/dt * trcrn(nt_fbri,n) + trcrn(bio_index(mm) + k-1,n) = c0 + enddo + enddo endif ! aicen > puny enddo ! ncat @@ -2177,6 +2186,93 @@ subroutine icepack_init_ocean_bio (amm, dmsp, dms, algalN, doc, dic, don, & end subroutine icepack_init_ocean_bio +! +!======================================================================= +! +! Given some added new ice to the base of the existing ice, recalculate +! vertical bio tracer so that new grid cells are all the same size. +! +! author: N. Jeffery, LANL +! date : Mar 3, 2024 +! +! Based on update_vertical_tracers modified for vertical biogeochemistry +! + subroutine update_vertical_bio_tracers(nbiolyr, trc, h1, h2, trc0, zspace) + + integer (kind=int_kind), intent(in) :: & + nbiolyr ! number of bio layers nblyr+1 + + real (kind=dbl_kind), dimension(:), intent(inout) :: & + trc ! vertical tracer + + real (kind=dbl_kind), intent(in) :: & + h1, & ! old thickness + h2, & ! new thickness + trc0 ! tracer value of added ice on ice bottom + + real (kind=dbl_kind), dimension(nbiolyr), intent(in) :: & + zspace + + ! local variables + + real(kind=dbl_kind), dimension(nbiolyr) :: trc2 ! updated tracer temporary + + ! vertical indices for old and new grid + integer :: k1, k2 + + real (kind=dbl_kind) :: & + z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom + z2a, z2b, & ! upper, lower boundary of new cell + overlap , & ! overlap between old and new cell + rnilyr + + !rnilyr = real(nilyr,dbl_kind) + z2a = c0 + z2b = c0 + ! loop over new grid cells + do k2 = 1, nbiolyr + + ! initialize new tracer + trc2(k2) = c0 + + ! calculate upper and lower boundary of new cell + z2a = z2b !((k2 - 1) * h2) * zspace(k2)+z2b ! / rnilyr + z2b = z2b + h2 * zspace(k2) !(k2 * h2) * zspace(k2)+z2a !/ rnilyr + + z1a = c0 + z1b = c0 + ! loop over old grid cells + do k1 = 1, nbiolyr + + ! calculate upper and lower boundary of old cell + z1a = z1b !((k1 - 1) * h1) * zspace(k1)+z1b !/ rnilyr + z1b = z1b + h1 * zspace(k1) !(k1 * h1) * zspace(k1)+z1a !/ rnilyr + + ! calculate overlap between old and new cell + overlap = max(min(z1b, z2b) - max(z1a, z2a), c0) + + ! aggregate old grid cell contribution to new cell + trc2(k2) = trc2(k2) + overlap * trc(k1) + + enddo ! k1 + + ! calculate upper and lower boundary of added new ice at bottom + z1a = h1 + z1b = h2 + + ! calculate overlap between added ice and new cell + overlap = max(min(z1b, z2b) - max(z1a, z2a), c0) + ! aggregate added ice contribution to new cell + trc2(k2) = trc2(k2) + overlap * trc0 + ! renormalize new grid cell + trc2(k2) = trc2(k2)/zspace(k2)/h2 !(rnilyr * trc2(k2)) / h2 + + enddo ! k2 + + ! update vertical tracer array with the adjusted tracer + trc = trc2 + + end subroutine update_vertical_bio_tracers !======================================================================= diff --git a/columnphysics/icepack_zbgc_shared.F90 b/columnphysics/icepack_zbgc_shared.F90 index 5ee9613a..e6099694 100644 --- a/columnphysics/icepack_zbgc_shared.F90 +++ b/columnphysics/icepack_zbgc_shared.F90 @@ -687,8 +687,8 @@ subroutine merge_bgc_fluxes (dt, nblyr, & !----------------------------------------------------------------- ! Merge fluxes !----------------------------------------------------------------- - dvssl = min(p5*vsnon/real(nslyr,kind=dbl_kind), hs_ssl*aicen) ! snow surface layer - dvint = vsnon - dvssl ! snow interior + dvssl = p5*vsnon/real(nslyr,kind=dbl_kind) !snow surface layer + dvint = vsnon - dvssl ! snow interior snow_bio_net(mm) = snow_bio_net(mm) & + trcrn(bio_index(mm)+nblyr+1)*dvssl & + trcrn(bio_index(mm)+nblyr+2)*dvint