Skip to content

Commit

Permalink
Bug fix in step therm2 for BGC
Browse files Browse the repository at this point in the history 
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
  • Loading branch information
@njeffery
njeffery committed Mar 25, 2024
1 parent b1e15e9 commit 8934f44
Show file tree
Hide file tree
Showing 4 changed files with 90 additions and 77 deletions.
101 changes: 39 additions & 62 deletions columnphysics/icepack_algae.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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, &
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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'
Expand Down Expand Up @@ -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) + &
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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'
Expand All @@ -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'
Expand All @@ -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'
Expand All @@ -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

Expand Down Expand Up @@ -2349,15 +2326,15 @@ 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
do k = 2, nblyr
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
Expand Down
32 changes: 31 additions & 1 deletion columnphysics/icepack_itd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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), &
Expand Down Expand Up @@ -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) &
Expand Down
11 changes: 6 additions & 5 deletions columnphysics/icepack_therm_itd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) :: &
Expand Down Expand Up @@ -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

Expand Down
23 changes: 14 additions & 9 deletions columnphysics/icepack_zbgc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -237,21 +239,21 @@ 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

! Diffuse_bio handles concentration changes from ice growth/melt
! 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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down

0 comments on commit 8934f44

Please sign in to comment.