Skip to content

Commit

Permalink
Fixes to conserve carbon in icepack
Browse files Browse the repository at this point in the history 
Bugfixes in bgc tracers
-corrected bgc ocean flux during lateral melt
-updated add_new_ice_bgc for compatibility with floe size distribution
-replace and greaty simplify adjust_tracer_profile with update_vertical_bio_tracers
-remove bgc tracers for small ice areas (corrects carbon conservation)
-pass carbon molar mass to icepack

Bugfixes in zaerosol snow tracers
-the two BGC snow layers are now fixed fractions of the snow thickness
-added missing snow bgc flux during ridging

Added diagnostics and updates to write statements

Verified carbon conservation in 20-day ocean-ice test.
nonBFB with bgc on. BFB with bgc off.
  • Loading branch information
@njeffery
njeffery committed Mar 12, 2024
1 parent abdd3d4 commit b1e15e9
Show file tree
Hide file tree
Showing 9 changed files with 423 additions and 271 deletions.
6 changes: 3 additions & 3 deletions columnphysics/icepack_aerosol.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 &
Expand Down Expand Up @@ -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
Expand Down
117 changes: 72 additions & 45 deletions columnphysics/icepack_algae.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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) :: &
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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, &
Expand All @@ -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, &
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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) :: &
Expand Down Expand Up @@ -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

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

Expand Down Expand Up @@ -1324,18 +1329,18 @@ 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)
write(warnstr,*)'sum_new,sum_tot,sum_initial,flux_bio(mm),source(mm):'
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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)'

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

Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
3 changes: 1 addition & 2 deletions columnphysics/icepack_brine.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) :: &
Expand Down Expand Up @@ -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, &
Expand Down
4 changes: 2 additions & 2 deletions columnphysics/icepack_itd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 + &
Expand Down
Loading

0 comments on commit b1e15e9

Please sign in to comment.