From 8934f4466b4d98bc7dcc7397e1a382396690d427 Mon Sep 17 00:00:00 2001 From: Nicole Jeffery Date: Mon, 25 Mar 2024 15:47:08 -0500 Subject: [PATCH] Bug fix in step therm2 for BGC Corrects division by zero in ice_algae and add_new_ice_bgc Remove small bgc values and track flux properly for conservation Remove small area bgc and account for fluxes Clean write statements nonBGC in ice bgc --- columnphysics/icepack_algae.F90 | 101 +++++++++++----------------- columnphysics/icepack_itd.F90 | 32 ++++++++- columnphysics/icepack_therm_itd.F90 | 11 +-- columnphysics/icepack_zbgc.F90 | 23 ++++--- 4 files changed, 90 insertions(+), 77 deletions(-) diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90 index 7db6b851..f841902f 100644 --- a/columnphysics/icepack_algae.F90 +++ b/columnphysics/icepack_algae.F90 @@ -251,9 +251,6 @@ subroutine zbio (dt, nblyr, & logical (kind=log_kind) :: & write_flux_diag - real (kind=dbl_kind) :: & - a_ice - character(len=*),parameter :: subname='(zbio)' zbgc_snown(:) = c0 @@ -269,7 +266,7 @@ subroutine zbio (dt, nblyr, & call bgc_carbon_sum(nblyr, hbri_old, trcrn(:), carbonInitial,n_doc,n_dic,n_algae,n_don) - if (aice_old > c0) then + if (aice_old > puny) then hsnow_i = vsno_old/aice_old do mm = 1,nbtrcr call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, & @@ -437,8 +434,6 @@ subroutine zbio (dt, nblyr, & if (write_flux_diag) then if (aicen > c0) then - if (n_cat .eq. 1) a_ice = c0 - a_ice = a_ice + aicen write(warnstr,*) subname, 'after merge_bgc_fluxes, n_cat:', n_cat call icepack_warnings_add(warnstr) do mm = 1,nbtrcr @@ -450,7 +445,7 @@ subroutine zbio (dt, nblyr, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'flux_bio_atm(mm)', flux_bio_atm(mm) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'flux_bio_atm(mm)*a_ice', flux_bio_atm(mm)*a_ice + write(warnstr,*) subname, 'flux_bio_atm(mm)*aicen', flux_bio_atm(mm)*aicen call icepack_warnings_add(warnstr) enddo endif @@ -1021,7 +1016,8 @@ 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 + flux_bio_tmp ! local parameters @@ -1061,6 +1057,7 @@ subroutine z_biogeochemistry (n_cat, dt, & mobile(:) = c0 conserve_C(:) = .true. nitrification(:) = c0 + flux_bio_tmp(:) = c0 do m = 1, nbtrcr do k = 1, nblyr+1 @@ -1070,7 +1067,13 @@ subroutine z_biogeochemistry (n_cat, dt, & iphin_N(k) = iphin(k) bphin_N(1) = bphi_min - in_init_cons(k,m) = trcrn(bio_index(m) + k-1)* hbri_old + if (abs(trcrn(bio_index(m) + k-1)) < accuracy) then + flux_bio_tmp(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 if (trcrn(bio_index(m) + k-1) < c0 ) then write(warnstr,*) subname,'zbgc initialization error, first ice = ', first_ice @@ -1329,7 +1332,7 @@ 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+source(mm))/dt-flux_bio(mm)) > max(accuracy, accuracy*abs(flux_bio(mm)))) & + if ((abs((sum_initial-sum_tot+source(mm))/dt-flux_bio(mm)) > max(puny, 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 @@ -1344,22 +1347,10 @@ subroutine z_biogeochemistry (n_cat, dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname,'sum_new,sum_old:',sum_new,sum_old call icepack_warnings_add(warnstr) - write(warnstr,*) subname,'mm,biocons(:):',mm,biocons(:) - call icepack_warnings_add(warnstr) - write(warnstr,*) subname,'biomat_low:',biomat_low - call icepack_warnings_add(warnstr) - write(warnstr,*) subname,'Diff(:):',Diff(:) - call icepack_warnings_add(warnstr) - write(warnstr,*) subname,'dmobile(:):',dmobile(:) - call icepack_warnings_add(warnstr) write(warnstr,*) subname,'mobile(mm):',mobile(mm) call icepack_warnings_add(warnstr) - write(warnstr,*) subname,'initcons_stationary(:):',initcons_stationary(:) - call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'trcrn(nt_zbgc_frac+mm-1):',trcrn(nt_zbgc_frac+mm-1) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'in_init_cons(:,mm):',in_init_cons(:,mm) - call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'exp_ret( mm),exp_rel( mm)',exp_ret( mm),exp_rel( mm) call icepack_warnings_add(warnstr) write(warnstr,*) subname,'darcyV,dhtop,dhbot' @@ -1420,22 +1411,23 @@ subroutine z_biogeochemistry (n_cat, dt, & bio_tmp = (biomat_brine(k,m) + react(k,m))*iphin_N(k) if (tr_bgc_C .and. m .eq. nlt_bgc_DIC(1) .and. bio_tmp .le. -accuracy) then ! satisfy DIC demands from ocean - write(warnstr,*) subname, 'DIC demand from ocean' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'm, k, nlt_bgc_DIC(1), bio_tmp, react(k,m):' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, m, k, nlt_bgc_DIC(1), bio_tmp, react(k,m) - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'flux_bio(m), hbri, hbri_old:' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, flux_bio(m), hbri, hbri_old - call icepack_warnings_add(warnstr) + !Uncomment for additional diagnostics + !write(warnstr,*) subname, 'DIC demand from ocean' + !call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, 'm, k, nlt_bgc_DIC(1), bio_tmp, react(k,m):' + !call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, m, k, nlt_bgc_DIC(1), bio_tmp, react(k,m) + !call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, 'flux_bio(m), hbri, hbri_old:' + !call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, flux_bio(m), hbri, hbri_old + !call icepack_warnings_add(warnstr) flux_bio(m) = flux_bio(m) + bio_tmp*dz(k)*hbri/dt bio_tmp = c0 - write(warnstr,*) subname, 'flux_bio(m) Final:' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, flux_bio(m) - call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, 'flux_bio(m) Final:' + !call icepack_warnings_add(warnstr) + !write(warnstr,*) subname, flux_bio(m) + !call icepack_warnings_add(warnstr) end if if (m .eq. nlt_bgc_Nit) then initcons_mobile(k) = max(c0,(biomat_brine(k,m)-nitrification(k) + & @@ -1477,6 +1469,7 @@ 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) @@ -1516,6 +1509,7 @@ subroutine z_biogeochemistry (n_cat, dt, & trcrn(nt_zbgc_frac+m-1) = zbgc_frac_init(m) if (sum_tot > c0) trcrn(nt_zbgc_frac+m-1) = sum_new/sum_tot end if + flux_bio(m) = flux_bio(m) + flux_bio_tmp(m) enddo ! m end subroutine z_biogeochemistry @@ -1785,6 +1779,7 @@ subroutine algal_dyn (dt, n_doc, n_dic, & rFep(:) = c0 DIC_r(:) = c0 DIC_s(:) = c0 + Zoo = c0 Nitin = ltrcrn(nlt_bgc_Nit) op_dep = c0 @@ -2233,9 +2228,9 @@ subroutine algal_dyn (dt, n_doc, n_dic, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, dN,secday,dt,n_doc call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)' + write(warnstr,*) subname, 'reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(n_algae))' call icepack_warnings_add(warnstr) - write(warnstr,*) subname, reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)) + write(warnstr,*) subname, reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(n_algae)) call icepack_warnings_add(warnstr) if (tr_bgc_Am) then write(warnstr,*) subname, 'reactb(nlt_bgc_Am),Am_r, Am_s' @@ -2256,9 +2251,9 @@ subroutine algal_dyn (dt, n_doc, n_dic, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, reactb(nlt_bgc_DOC(k)) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'DOC_r,DOC_s' + write(warnstr,*) subname, 'DOC_r(k),DOC_s(k),k' call icepack_warnings_add(warnstr) - write(warnstr,*) subname, DOC_r(k),DOC_s(k) + write(warnstr,*) subname, DOC_r(k),DOC_s(k),k end do do k = 1,n_dic write(warnstr,*) subname, 'DICin' @@ -2269,9 +2264,9 @@ subroutine algal_dyn (dt, n_doc, n_dic, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, reactb(nlt_bgc_DIC(k)) call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'DIC_r,DIC_s' + write(warnstr,*) subname, 'DIC_r(k),DIC_s(k),k' call icepack_warnings_add(warnstr) - write(warnstr,*) subname, DIC_r(k),DIC_s(k) + write(warnstr,*) subname, DIC_r(k),DIC_s(k),k end do call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Zoo' @@ -2280,24 +2275,6 @@ subroutine algal_dyn (dt, n_doc, n_dic, & call icepack_warnings_add(warnstr) endif endif - ! if (abs(Cerror) > max(reactb(:))*1.0e-5) then - ! conserve_C = .false. - ! write(warnstr,*) subname, 'Conservation error!' - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, 'Cerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc' - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, Cerror,dN, DONin(1),kn_bac(1),secday,dt,n_doc - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, 'reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)' - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, reactb(nlt_bgc_Nit),reactb(nlt_bgc_N(1)),reactb(nlt_bgc_N(2)) - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, 'reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)), DON_r(1),DON_s(1)' - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, reactb(nlt_bgc_Am),reactb(nlt_bgc_DON(1)),DON_r(1),DON_s(1) - ! call icepack_warnings_add(warnstr) - ! write(warnstr,*) subname, 'Zoo:',Zoo - ! endif end subroutine algal_dyn @@ -2349,7 +2326,7 @@ subroutine thin_ice_flux (hin, hin_old, Cin, flux_o_tot, & sum_bio = c0 dh = hin-hin_old - if (dh .le. c0) then ! keep the brine concentration fixed + if (dh .le. c0 .and. hin_old > puny) then ! keep the brine concentration fixed sum_bio = (Cin(1)+Cin(nblyr+1))/hin_old*zspace*p5 Cin(1) = Cin(1)/hin_old*hin Cin(nblyr+1) = Cin(nblyr+1)/hin_old*hin @@ -2357,7 +2334,7 @@ subroutine thin_ice_flux (hin, hin_old, Cin, flux_o_tot, & sum_bio = sum_bio + Cin(k)/hin_old*zspace Cin(k) = Cin(k)/hin_old*hin + dC enddo - else + else ! spread evenly in ice layers dC = dh*ocean_bio do k = 1, nblyr+1 Cin(k) = Cin(k) + dC diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index dc154659..d12f74f1 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -1096,6 +1096,7 @@ subroutine zap_small_areas (dt, ntrcr, & blevels real (kind=dbl_kind) :: xtmp, sicen ! temporary variables + real (kind=dbl_kind) :: dvssl, dvint ! temporary variables real (kind=dbl_kind) , dimension (1):: trcr_skl real (kind=dbl_kind) , dimension (nblyr+1):: bvol @@ -1146,7 +1147,6 @@ subroutine zap_small_areas (dt, ntrcr, & if (skl_bgc .and. nbtrcr > 0) then blevels = 1 bvol(1) = aicen(n)*sk_l - it = 1 do it = 1, nbtrcr trcr_skl(1) = trcrn(bio_index(it),n) call zap_small_bgc(blevels, dflux_bio(it), & @@ -1273,6 +1273,36 @@ subroutine zap_small_areas (dt, ntrcr, & enddo ! it endif + if (skl_bgc .and. nbtrcr > 0) then + blevels = 1 + bvol(1) = (aice-c1)/aice * sk_l + do it = 1, nbtrcr + trcr_skl(1) = trcrn(bio_index(it),n) + call zap_small_bgc(blevels, dflux_bio(it), & + dt, bvol(1:blevels), trcr_skl(blevels)) + enddo + elseif (z_tracers .and. nbtrcr > 0) then + blevels = nblyr + 1 + bvol(:) = (aice-c1)/aice*vicen(n)/real(nblyr,kind=dbl_kind)*trcrn(nt_fbri,n) + bvol(1) = p5*bvol(1) + bvol(blevels) = p5*bvol(blevels) + do it = 1, nbtrcr + call zap_small_bgc(blevels, dflux_bio(it), & + dt, bvol(1:blevels),trcrn(bio_index(it):bio_index(it)+blevels-1,n)) + if (icepack_warnings_aborted(subname)) return + enddo + ! zap snow zaerosols + dvssl = p5*vsnon(n)/real(nslyr,kind=dbl_kind) ! snow surface layer + dvint = vsnon(n) - dvssl ! snow interior + + do it = 1, nbtrcr + xtmp = (trcrn(bio_index(it)+nblyr+1,n)*dvssl + & + trcrn(bio_index(it)+nblyr+2,n)*dvint)*(aice-c1)/aice/dt + dflux_bio(it) = dflux_bio(it) + xtmp + enddo ! it + + endif + if (tr_iso) then do it = 1, n_iso xtmp = (vsnon(n)*trcrn(nt_isosno+it-1,n) & diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 5267af04..829352fe 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -92,10 +92,10 @@ subroutine linear_itd (ncat, hin_max, & ntrcr, trcr_depend, & trcr_base, n_trcr_strata,& nt_strata, & - aicen_init, vicen_init, & - aicen, trcrn, & - vicen, vsnon, & - aice, aice0, & + aicen_init, vicen_init, & + aicen, trcrn, & + vicen, vsnon, & + aice, aice0, & fpond, Tf ) integer (kind=int_kind), intent(in) :: & @@ -1949,7 +1949,8 @@ subroutine add_new_ice (ncat, nilyr, & aicen_init, vicen_init, vi0_init, & aicen, vicen, vin0new, & ntrcr, trcrn, nbtrcr, & - ocean_bio, flux_bio, hsurp) + ocean_bio, flux_bio, hsurp, & + d_an_tot) if (icepack_warnings_aborted(subname)) return endif diff --git a/columnphysics/icepack_zbgc.F90 b/columnphysics/icepack_zbgc.F90 index 5503469d..d33b2989 100644 --- a/columnphysics/icepack_zbgc.F90 +++ b/columnphysics/icepack_zbgc.F90 @@ -92,7 +92,8 @@ subroutine add_new_ice_bgc (dt, nblyr, ncats, & aicen_init, vicen_init, vi0_init, & aicen, vicen, vin0new, & ntrcr, trcrn, nbtrcr, & - ocean_bio, flux_bio, hsurp) + ocean_bio, flux_bio, hsurp, & + d_an_tot) integer (kind=int_kind), intent(in) :: & nblyr , & ! number of bio layers @@ -118,7 +119,8 @@ subroutine add_new_ice_bgc (dt, nblyr, ncats, & aicen_init , & ! initial concentration of ice vicen_init , & ! initial volume per unit area of ice (m) aicen , & ! concentration of ice - vicen ! volume per unit area of ice (m) + vicen , & ! volume per unit area of ice (m) + d_an_tot real (kind=dbl_kind), dimension (:,:), intent(inout) :: & trcrn ! ice tracers @@ -215,9 +217,9 @@ subroutine add_new_ice_bgc (dt, nblyr, ncats, & vbrin(n) = vbrin(n) + vsurp vice_new = vicen_init(n) + vsurp - if (tr_brine .and. vice_new > c0) then + if (tr_brine .and. vice_new > puny) then !c0) then trcrn(nt_fbri,n) = vbrin(n)/vice_new - elseif (tr_brine .and. vice_new <= c0) then + elseif (tr_brine .and. vicen(n) <= c0) then trcrn(nt_fbri,n) = c1 endif @@ -237,13 +239,13 @@ subroutine add_new_ice_bgc (dt, nblyr, ncats, & ! Combine bgc in new ice grown in open water with category 1 ice. !----------------------------------------------------------------- do n = 1, ncats - if (vin0new(n) > c0) then + if (vin0new(n) > c0 .and. d_an_tot(n) > c0) then vbri1 = vbrin(n) vbrin(n) = vbrin(n) + vin0new(n) - if (tr_brine .and. vicen(n) > c0) then + if (tr_brine .and. vicen(n) > puny) then trcrn(nt_fbri,n) = vbrin(n)/vicen(n) - elseif (tr_brine .and. vicen(n) <= c0) then + elseif (tr_brine .and. vicen(n) <= puny) then trcrn(nt_fbri,n) = c1 endif @@ -251,7 +253,7 @@ subroutine add_new_ice_bgc (dt, nblyr, ncats, & ! ice area changes ! add salt throughout, location = 0 - if (nbtrcr > 0 .and. vbrin(n) > c0) then + if (nbtrcr > 0 .and. vbrin(n) > puny) then do m = 1, nbtrcr bio0new = ocean_bio(m)*zbgc_init_frac(m) do k = 1, nblyr+1 @@ -2229,6 +2231,7 @@ subroutine update_vertical_bio_tracers(nbiolyr, trc, h1, h2, trc0, zspace) !rnilyr = real(nilyr,dbl_kind) z2a = c0 z2b = c0 + if (h2 > c0) then ! loop over new grid cells do k2 = 1, nbiolyr @@ -2268,7 +2271,9 @@ subroutine update_vertical_bio_tracers(nbiolyr, trc, h1, h2, trc0, zspace) trc2(k2) = trc2(k2)/zspace(k2)/h2 !(rnilyr * trc2(k2)) / h2 enddo ! k2 - + else + trc2 = trc + endif ! update vertical tracer array with the adjusted tracer trc = trc2