diff --git a/Registry.xml b/Registry.xml
index a32576e843ac..c5cccefb8097 100644
--- a/Registry.xml
+++ b/Registry.xml
@@ -2649,12 +2649,12 @@
packages="pkgTracerVerticalIron"
/>
diff --git a/column/ice_aerosol.F90 b/column/ice_aerosol.F90
index b5f5e8007de9..82a482e49c9c 100644
--- a/column/ice_aerosol.F90
+++ b/column/ice_aerosol.F90
@@ -455,12 +455,13 @@ subroutine update_snow_bgc (dt, nblyr, &
vice_old, vsno_old, &
vicen, vsnon, &
aicen, flux_bio_atm,&
- zbgc_atm, flux_bio)
+ zbgc_atm, flux_bio, &
+ bio_index_o)
- use ice_colpkg_shared, only: hi_ssl, hs_ssl
+ use ice_colpkg_shared, only: hi_ssl, hs_ssl, hs_ssl_min
use ice_constants_colpkg, only: c0, rhos, rhoi, hs_min, puny, &
- c2, c1
- use ice_zbgc_shared, only: kscavz
+ c2, c1, p5
+ use ice_zbgc_shared, only: kscavz
integer (kind=int_kind), intent(in) :: &
nbtrcr, & ! number of distinct snow tracers
@@ -469,7 +470,8 @@ subroutine update_snow_bgc (dt, nblyr, &
ntrcr ! number of tracers
integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
- bio_index
+ bio_index, &
+ bio_index_o & ! provides index of scavenging (kscavz) data array
real (kind=dbl_kind), intent(in) :: &
dt ! time step
@@ -488,9 +490,11 @@ subroutine update_snow_bgc (dt, nblyr, &
vice_old, &
vsno_old
- real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
+ real (kind=dbl_kind),dimension(nbtrcr), intent(out) :: &
zbgc_snow, & ! aerosol contribution from snow to ice
- zbgc_atm, & ! and atm to ice concentration * volume (kg or mmol/m^3*m)
+ zbgc_atm ! and atm to ice concentration * volume (kg or mmol/m^3*m)
+
+ real (kind=dbl_kind),dimension(nbtrcr), intent(inout) :: &
flux_bio ! total ocean tracer flux (mmol/m^2/s)
real (kind=dbl_kind), dimension(nbtrcr), &
@@ -508,6 +512,8 @@ subroutine update_snow_bgc (dt, nblyr, &
real (kind=dbl_kind) :: &
dzssl, dzssl_new, & ! snow ssl thickness
dzint, dzint_new, & ! snow interior thickness
+ hi, & ! ice thickness (m)
+ hilyr, & ! ice layer thickness (m)
hs, & ! snow thickness (m)
dhs_evap, & ! snow thickness change due to evap
dhs_melts, & ! ... due to surface melt
@@ -541,6 +547,11 @@ subroutine update_snow_bgc (dt, nblyr, &
zbgc_atm(:) = c0
hs_old = vsno_old/aice_old
+ if (aice_old .gt. puny) then
+ hs_old = vsno_old/aice_old
+ else
+ hs_old = c0
+ end if
hslyr_old = hs_old/real(nslyr,kind=dbl_kind)
dzssl = min(hslyr_old/c2, hs_ssl)
@@ -549,40 +560,46 @@ subroutine update_snow_bgc (dt, nblyr, &
if (aicen > c0) then
ar = c1/aicen
hs = vsnon*ar
- dhs_melts = -melts
- dhs_snoice = snoice*rhoi/rhos
+ hi = vicen*ar
else ! ice disappeared during time step
- ar = c1
- hs = vsnon/aice_old
- dhs_melts = -melts
- dhs_snoice = snoice*rhoi/rhos
+ ar = c1
+ hs = c0
+ hi = c0
+ if (aice_old > c0) hs = vsnon/aice_old
endif
-
+ hilyr = hi/real(nblyr,kind=dbl_kind)
+ hslyr = hs/real(nslyr,kind=dbl_kind)
+ dzssl_new = min(hslyr/c2, hs_ssl)
+ dhs_melts = -melts
+ dhs_snoice = snoice*rhoi/rhos
dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice &
+ fsnow/rhos*dt)
! trcrn() has units kg/m^3
- if ((vsno_old .le. puny) .or. (vsnon .le. puny)) then
-
+ if (dzssl_new .lt. hs_ssl_min) then ! Put atm BC/dust flux directly into the sea ice
do k=1,nbtrcr
flux_bio(k) = flux_bio(k) + &
(trcrn(bio_index(k)+ nblyr+1)*dzssl+ &
- trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
+ trcrn(bio_index(k)+ nblyr+2)*dzint)/dt
trcrn(bio_index(k) + nblyr+1) = c0
trcrn(bio_index(k) + nblyr+2) = c0
- zbgc_atm(k) = zbgc_atm(k) &
- + flux_bio_atm(k)*dt
+ if (hilyr .lt. hs_ssl_min) then
+ flux_bio(k) = flux_bio(k) + flux_bio_atm(k)
+ else
+ zbgc_atm(k) = zbgc_atm(k) &
+ + flux_bio_atm(k)*dt
+ end if
enddo
- else
-
+ else
+
do k=1,nbtrcr
flux_bio_o(k) = flux_bio(k)
aerosno (k,1) = trcrn(bio_index(k)+ nblyr+1) * dzssl
aerosno (k,2) = trcrn(bio_index(k)+ nblyr+2) * dzint
aerosno0(k,:) = aerosno(k,:)
- aerotot0(k) = aerosno(k,2) + aerosno(k,1)
+ aerotot0(k) = aerosno(k,2) + aerosno(k,1)
enddo
!-------------------------------------------------------------------
@@ -590,7 +607,47 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
dzint = dzint + min(dzssl + dhs_evap, c0)
dzssl = max(dzssl + dhs_evap, c0)
-
+ if (dzssl <= puny) then
+ do k = 1,nbtrcr
+ aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
+ aerosno(k,1) = c0
+ end do
+ end if
+ if (dzint <= puny) then
+ do k = 1,nbtrcr
+ flux_bio(k) = flux_bio(k) + (aerosno(k,2) + aerosno(k,1))/dt
+ aerosno(k,2) = c0
+ aerosno(k,1) = c0
+ end do
+ end if
+ !------------------------------------------------------------------
+ ! snowfall
+ !-------------------------------------------------------------------
+ if (fsnow > c0) then
+ sloss1 = c0
+ dz = min(fsnow/rhos*dt,dzssl)
+ do k = 1, nbtrcr
+ if (dzssl > puny) &
+ sloss1 = aerosno(k,1)*dz/dzssl
+ aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
+ aerosno(k,2) = aerosno(k,2) + sloss1
+ end do
+ dzssl = dzssl - dz + fsnow/rhos*dt
+ dzint = dzint + dz
+ end if
+ if (dzssl <= puny) then
+ do k = 1,nbtrcr
+ aerosno(k,2) = aerosno(k,2) + aerosno(k,1)
+ aerosno(k,1) = c0
+ end do
+ end if
+ if (dzint <= puny) then
+ do k = 1,nbtrcr
+ flux_bio(k) = flux_bio(k) + (aerosno(k,2) + aerosno(k,1))/dt
+ aerosno(k,2) = c0
+ aerosno(k,1) = c0
+ end do
+ end if
!-------------------------------------------------------------------
! surface snow melt
!-------------------------------------------------------------------
@@ -598,39 +655,37 @@ subroutine update_snow_bgc (dt, nblyr, &
do k = 1, nbtrcr
sloss1 = c0
sloss2 = c0
- if (dzssl > puny) &
- sloss1 = kscavz(k)*aerosno(k,1) &
- *min(-dhs_melts,dzssl)/dzssl
- aerosno(k,1) = aerosno(k,1) - sloss1
- if (dzint > puny) &
- sloss2 = kscavz(k)*aerosno(k,2) &
- *max(-dhs_melts-dzssl,c0)/dzint
- aerosno(k,2) = aerosno(k,2) - sloss2
- zbgc_snow(k) = zbgc_snow(k) + (sloss1+sloss2)
- enddo !
+ if (dzssl > puny) &
+ sloss1 = kscavz(bio_index_o(k))*aerosno(k,1) &
+ *min(-dhs_melts,dzssl)/dzssl
+ aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
+ if (dzint > puny) &
+ sloss2 = kscavz(bio_index_o(k))*aerosno(k,2) &
+ *max(-dhs_melts-dzssl,c0)/dzint
+ aerosno(k,2) = max(c0,aerosno(k,2) - sloss2)
+ flux_bio(k) = flux_bio(k) + (sloss1+sloss2)/dt ! all not scavenged ends in ocean
+ enddo
! update snow thickness
dzint=dzint+min(dzssl+dhs_melts, c0)
dzssl=max(dzssl+dhs_melts, c0)
- if ( dzssl <= puny ) then ! ssl melts away
- aerosno(:,2) = aerosno(:,1) + aerosno(:,2)
- aerosno(:,1) = c0
+ if ( dzssl .le. puny ) then ! ssl melts away
+ do k = 1,nbtrcr
+ aerosno(k,2) = aerosno(k,1) + aerosno(k,2)
+ aerosno(k,1) = c0
+ end do
dzssl = max(dzssl, c0)
endif
- if (dzint <= puny ) then ! all snow melts away
- zbgc_snow(:) = zbgc_snow(:) &
- + max(c0,aerosno(:,1) + aerosno(:,2))
- aerosno(:,:) = c0
+ if (dzint .le. puny ) then ! all snow melts away
+ do k = 1,nbtrcr
+ zbgc_snow(k) = zbgc_snow(k) &
+ + aerosno(k,1) + aerosno(k,2)
+ aerosno(k,:) = c0
+ enddo
dzint = max(dzint, c0)
endif
- endif
-
- !-------------------------------------------------------------------
- ! snowfall
- !-------------------------------------------------------------------
- if (fsnow > c0) dzssl = dzssl + fsnow/rhos*dt
-
+ endif ! -dhs_melts > puny
!-------------------------------------------------------------------
! snow-ice formation
!-------------------------------------------------------------------
@@ -638,39 +693,46 @@ subroutine update_snow_bgc (dt, nblyr, &
do k = 1, nbtrcr
sloss1 = c0
sloss2 = c0
- if (dzint > puny) &
- sloss2 = min(dhs_snoice, dzint) &
- *aerosno(k,2)/dzint
- aerosno(k,2) = aerosno(k,2) - sloss2
- if (dzssl > puny) &
+ if (dzint > puny .and. aerosno(k,2) > c0) &
+ sloss2 = min(dhs_snoice, dzint) &
+ *aerosno(k,2)/dzint
+ aerosno(k,2) = max(c0,aerosno(k,2) - sloss2)
+ if (dzssl > puny .and. aerosno(k,1) > c0) &
sloss1 = max(dhs_snoice-dzint, c0) &
*aerosno(k,1)/dzssl
- aerosno(k,1) = aerosno(k,1) - sloss1
+
+ aerosno(k,1) = max(c0,aerosno(k,1) - sloss1)
+ flux_bio(k) = flux_bio(k) &
+ + kscavz(bio_index_o(k)) * (sloss2+sloss1)/dt
zbgc_snow(k) = zbgc_snow(k) &
- + (sloss2+sloss1)
+ + (c1-kscavz(bio_index_o(k)))*(sloss2+sloss1)
enddo
- dzssl = dzssl - max(dhs_snoice-dzint, c0)
+ dzssl = max(c0,dzssl - max(dhs_snoice-dzint, c0))
dzint = max(dzint-dhs_snoice, c0)
- endif
+ endif ! dhs_snowice > puny
!-------------------------------------------------------------------
! aerosol deposition
!-------------------------------------------------------------------
- if (aicen > c0) then
- hs = vsnon * ar
- else
- hs = c0
- endif
- if (hs >= hs_min) then !should this really be hs_min or 0?
- ! should use same hs_min value as in radiation
+ ! if (aicen > c0) then
+ ! hs = vsnon * ar
+ ! else
+ ! hs = c0
+ ! endif
+ ! Spread out the atm dust flux in the snow interior for small snow surface layers
+ if (dzssl .ge. hs_ssl*p5) then
+
do k=1,nbtrcr
aerosno(k,1) = aerosno(k,1) &
+ flux_bio_atm(k)*dt
enddo
- else
+ else
+ dz = (hs_ssl*p5 - dzssl)/(hs_ssl*p5)
do k=1,nbtrcr
- zbgc_atm(k) = zbgc_atm(k) &
- + flux_bio_atm(k)*dt
+ aerosno(k,1) = aerosno(k,1) &
+ + flux_bio_atm(k)*dt*(c1-dz)
+ aerosno(k,2) = aerosno(k,2) &
+ + flux_bio_atm(k)*dt*dz
enddo
endif
@@ -690,30 +752,31 @@ subroutine update_snow_bgc (dt, nblyr, &
endif
if (dzint <= puny) then ! nothing in Snow Int
do k = 1, nbtrcr
- zbgc_snow(k) = zbgc_snow(k) + max(c0,aerosno(k,2))
+ zbgc_snow(k) = zbgc_snow(k) + max(c0,aerosno(k,2)+aerosno(k,1))
+ aerosno(k,1) = c0
aerosno(k,2) = c0
enddo
endif
hslyr = hs/real(nslyr,kind=dbl_kind)
dzssl_new = min(hslyr/c2, hs_ssl)
- dzint_new = hs - dzssl_new
+ dzint_new = max(c0,hs - dzssl_new)
- if (hs > hs_min) then !should this really be hs_min or 0?
+ if (hs > hs_min) then
do k = 1, nbtrcr
dznew = min(dzssl_new-dzssl, c0)
sloss1 = c0
- if (dzssl > puny) &
+ if (dzssl > puny .and. aerosno(k,1) > c0) &
sloss1 = dznew*aerosno(k,1)/dzssl ! not neccesarily a loss
- dznew = max(dzssl_new-dzssl, c0)
- if (dzint > puny) &
- sloss1 = sloss1 + aerosno(k,2)*dznew/dzint
- aerosno(k,1) = aerosno(k,1) + sloss1
- aerosno(k,2) = aerosno(k,2) - sloss1
+ dznew = max(dzssl_new-dzssl, c0)
+ if (dzint > puny .and. aerosno(k,2) > c0) &
+ sloss1 = aerosno(k,2)*dznew/dzint
+ aerosno(k,1) = max(c0,aerosno(k,1) + sloss1)
+ aerosno(k,2) = max(c0,aerosno(k,2) - sloss1)
enddo
else
zbgc_snow(:) = zbgc_snow(:) &
- + max(c0,aerosno(:,1) + aerosno(:,2))
+ + aerosno(:,1) + aerosno(:,2)
aerosno(:,:) = c0
endif
@@ -722,11 +785,16 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
do k = 1, nbtrcr
aerotot(k) = aerosno(k,2) + aerosno(k,1) &
- + zbgc_snow(k) + zbgc_atm(k)
+ + zbgc_snow(k) + zbgc_atm(k)
aero_cons(k) = aerotot(k)-aerotot0(k) &
- - ( flux_bio_atm(k) &
+ - ( flux_bio_atm(k) &
- (flux_bio(k)-flux_bio_o(k))) * dt
- if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
+ if (aerotot0(k) > aerotot(k) .and. aerotot0(k) > c0) then
+ aero_cons(k) = aero_cons(k)/aerotot0(k)
+ else if (aerotot(k) > c0) then
+ aero_cons(k) = aero_cons(k)/aerotot(k)
+ end if
+ if (aero_cons(k) > puny .or. zbgc_snow(k) + zbgc_atm(k) < c0) then
write(warning,*) 'Conservation failure: aerosols in snow'
call add_warning(warning)
write(warning,*) 'test aerosol 1'
@@ -754,14 +822,13 @@ subroutine update_snow_bgc (dt, nblyr, &
!-------------------------------------------------------------------
! reload tracers
!-------------------------------------------------------------------
- if (vsnon > puny) then
+ if (dzssl_new > puny .and. dzint_new > puny .and. vsnon > puny) then
do k = 1,nbtrcr
- trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
+ trcrn(bio_index(k)+nblyr+1)=aerosno(k,1)/dzssl_new
trcrn(bio_index(k)+nblyr+2)=aerosno(k,2)/dzint_new
enddo
else
do k = 1,nbtrcr
- zbgc_snow(k) = (zbgc_snow(k) + aerosno(k,1) + aerosno(k,2))
trcrn(bio_index(k)+nblyr+1)= c0
trcrn(bio_index(k)+nblyr+2)= c0
enddo
diff --git a/column/ice_algae.F90 b/column/ice_algae.F90
index bd97364d26ff..5851d7d83ae6 100644
--- a/column/ice_algae.F90
+++ b/column/ice_algae.F90
@@ -33,11 +33,11 @@ subroutine zbio (dt, nblyr, &
snoice, nbtrcr, &
fsnow, ntrcr, &
trcrn, bio_index, &
- aice_old, &
+ bio_index_o, aice_old, &
vice_old, vsno_old, &
vicen, vsnon, &
- aicen, flux_bio_atm,&
- n_cat, n_algae, &
+ aicen, flux_bio_atm,&
+ n_cat, n_algae, &
n_doc, n_dic, &
n_don, &
n_fed, n_fep, &
@@ -81,7 +81,8 @@ subroutine zbio (dt, nblyr, &
ntrcr ! number of tracers
integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
- bio_index ! references index of bio tracer (nbtrcr) to tracer array (ntrcr)
+ bio_index, & ! references index of bio tracer (nbtrcr) to tracer array (ntrcr)
+ bio_index_o, ! references index of data arrays (eg. kscavz)
real (kind=dbl_kind), intent(in) :: &
dt, & ! time step
@@ -127,13 +128,13 @@ subroutine zbio (dt, nblyr, &
real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
igrid , & ! biology vertical interface points
iTin , & ! salinity vertical interface points
- iphin , & ! Porosity on the igrid
+ iphin , & ! Porosity on the igrid
iDin ! Diffusivity/h on the igrid (1/s)
-
+
real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
- cgrid , & ! CICE vertical coordinate
- icgrid , & ! CICE interface coordinate
- fswthrul ! visible short wave radiation on icgrid (W/m^2)
+ cgrid , & ! CICE vertical coordinate
+ icgrid , & ! CICE interface coordinate
+ fswthrul ! visible short wave radiation on icgrid (W/m^2)
real (kind=dbl_kind), dimension(:), &
intent(in) :: &
@@ -141,26 +142,26 @@ subroutine zbio (dt, nblyr, &
real (kind=dbl_kind), dimension(ntrcr), &
intent(inout) :: &
- trcrn
+ trcrn
- real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
- zfswin ! visible Short wave flux on igrid (W/m^2)
-
- real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
+ real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
+ zfswin ! visible Short wave flux on igrid (W/m^2)
+
+ real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
Zoo ! N losses to the system from reaction terms
- ! (ie. zooplankton/bacteria) (mmol/m^3)
+ ! (ie. zooplankton/bacteria) (mmol/m^3)
- real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
+ real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
!change to inout when updating ocean fields
- ocean_bio ! ocean concentrations (mmol/m^3)
+ ocean_bio ! ocean concentrations (mmol/m^3)
real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
bphin ! Porosity on the bgrid
- real (kind=dbl_kind), intent(inout):: &
+ real (kind=dbl_kind), intent(inout):: &
PP_net , & ! net PP (mg C/m^2/d) times aice
grow_net , & ! net specific growth (m/d) times vice
- upNO , & ! tot nitrate uptake rate (mmol/m^2/d) times aice
+ upNO , & ! tot nitrate uptake rate (mmol/m^2/d) times aice
upNH , & ! tot ammonium uptake rate (mmol/m^2/d) times aice
totalChla ! total chla (mg chla/m^2)
@@ -213,7 +214,7 @@ subroutine zbio (dt, nblyr, &
accuracy = 1.0e-13_dbl_kind
character(len=char_len_long) :: &
- warning
+ warning
real (kind=dbl_kind), dimension (nblyr+1) :: &
zspace ! vertical grid spacing
@@ -255,25 +256,26 @@ subroutine zbio (dt, nblyr, &
vice_old, vsno_old, &
vicen, vsnon, &
aicen, flux_bio_atm, &
- zbgc_atmn, flux_bio_sno)
+ zbgc_atmn, flux_bio_snom &
+ bio_index_o)
call z_biogeochemistry (n_cat, dt, &
nilyr, nslyr, &
nblyr, nbtrcr, &
- n_algae, n_doc, &
+ n_algae, n_doc, &
n_dic, n_don, &
n_fed, n_fep, &
n_zaero, first_ice, &
- aicen, vicen, &
- hice_old, ocean_bio, &
+ aicen, vicen, &
+ hice_old, ocean_bio, &
flux_bion, bphin, &
- iphin, trcrn, &
+ iphin, trcrn, &
iDin, sss, &
fswthrul, grow_alg, &
upNOn, upNHn, &
dh_top, dh_bot, &
dh_top_chl, dh_bot_chl,&
- zfswin, hbri, &
+ zfswin, hbri, &
hbri_old, darcy_V, &
darcy_V_chl, bgrid, &
igrid, icgrid, &
diff --git a/column/ice_colpkg.F90 b/column/ice_colpkg.F90
index 284c55a2eabe..285d43984026 100644
--- a/column/ice_colpkg.F90
+++ b/column/ice_colpkg.F90
@@ -5394,10 +5394,10 @@ subroutine colpkg_biogeochemistry(dt, &
use ice_algae, only: zbio, sklbio
use ice_brine, only: preflushing_changes, compute_microS_mushy, &
- update_hbrine, compute_microS
+ update_hbrine, compute_microS
use ice_colpkg_shared, only: solve_zsal, z_tracers, phi_snow
use ice_colpkg_tracers, only: nt_fbri, tr_brine, &
- nt_bgc_S, nt_qice, nt_sice, nt_zbgc_frac, bio_index
+ nt_bgc_S, nt_qice, nt_sice, nt_zbgc_frac, bio_index, bio_index_o
use ice_constants_colpkg, only: c0, c1, puny, p5
use ice_zsalinity, only: zsalinity
use ice_zbgc_shared, only: zbgc_frac_init
@@ -5419,7 +5419,7 @@ subroutine colpkg_biogeochemistry(dt, &
real (kind=dbl_kind), dimension (:), intent(inout) :: &
bgrid , & ! biology nondimensional vertical grid points
igrid , & ! biology vertical interface points
- cgrid , & ! CICE vertical coordinate
+ 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
@@ -5428,9 +5428,9 @@ subroutine colpkg_biogeochemistry(dt, &
dhbr_bot , & ! brine bottom change
darcy_V , & ! darcy velocity positive up (m/s)
hin_old , & ! old ice thickness
- sice_rho , & ! avg sea ice density (kg/m^3)
- ice_bio_net , & ! depth integrated tracer (mmol/m^2)
- snow_bio_net , & ! depth integrated snow tracer (mmol/m^2)
+ sice_rho , & ! avg sea ice density (kg/m^3)
+ 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
logical (kind=log_kind), dimension (:), intent(inout) :: &
@@ -5445,18 +5445,18 @@ subroutine colpkg_biogeochemistry(dt, &
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
Zoo , & ! N losses accumulated in timestep (ie. zooplankton/bacteria)
! mmol/m^3
- bphi , & ! porosity of layers
+ bphi , & ! porosity of layers
bTiz , & ! layer temperatures interpolated on bio grid (C)
zfswin , & ! Shortwave flux into layers interpolated on bio grid (W/m^2)
- iDi , & ! igrid Diffusivity (m^2/s)
- iki , & ! Ice permeability (m^2)
- trcrn ! tracers
+ iDi , & ! igrid Diffusivity (m^2/s)
+ iki , & ! Ice permeability (m^2)
+ trcrn ! tracers
real (kind=dbl_kind), intent(inout) :: &
grow_net , & ! Specific growth rate (/s) per grid cell
PP_net , & ! Total production (mg C/m^2/s) per grid cell
hbri , & ! brine height, area-averaged for comparison with hi (m)
- zsal_tot , & ! Total ice salinity in per grid cell (g/m^2)
+ zsal_tot , & ! Total ice salinity in per grid cell (g/m^2)
fzsal , & ! Total flux of salt to ocean at time step for conservation
fzsal_g , & ! Total gravity drainage flux
upNO , & ! nitrate uptake rate (mmol/m^2/d) times aice
@@ -5464,7 +5464,7 @@ subroutine colpkg_biogeochemistry(dt, &
totalChla ! ice integrated chla and summed over all algal groups (mg/m^2)
logical (kind=log_kind), intent(inout) :: &
- Rayleigh_criteria ! .true. means Ra_c was reached
+ Rayleigh_criteria ! .true. means Ra_c was reached
real (kind=dbl_kind), dimension (:,:), intent(in) :: &
fswpenln ! visible SW entering ice layers (W m-2)
@@ -5476,11 +5476,11 @@ subroutine colpkg_biogeochemistry(dt, &
meltbn , & ! bottom melt in category n (m)
congeln , & ! congelation ice formation in category n (m)
snoicen , & ! snow-ice formation in category n (m)
- salinz , & ! initial salinity profile (ppt)
- flux_bio_atm, & ! all bio fluxes to ice from atmosphere
+ salinz , & ! initial salinity profile (ppt)
+ flux_bio_atm, & ! all bio fluxes to ice from atmosphere
aicen_init , & ! initial ice concentration, for linear ITD
vicen_init , & ! initial ice volume (m), for linear ITD
- vsnon_init , & ! initial snow volume (m), for aerosol
+ vsnon_init , & ! initial snow volume (m), for aerosol
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
@@ -5496,7 +5496,7 @@ subroutine colpkg_biogeochemistry(dt, &
logical (kind=log_kind), intent(in) :: &
skl_bgc ! if true, solve skeletal biochemistry
- logical (kind=log_kind), intent(inout) :: &
+ logical (kind=log_kind), intent(inout) :: &
l_stop ! if true, abort the model
character (len=*), intent(inout) :: stop_label
@@ -5513,7 +5513,7 @@ subroutine colpkg_biogeochemistry(dt, &
hbr_old , & ! old brine thickness before growh/melt
dhice , & ! change due to sublimation/condensation (m)
kavg , & ! average ice permeability (m^2)
- bphi_o , & ! surface ice porosity
+ bphi_o , & ! surface ice porosity
hbrin , & ! brine height
dh_direct ! surface flooding or runoff
@@ -5525,20 +5525,20 @@ subroutine colpkg_biogeochemistry(dt, &
real (kind=dbl_kind), dimension (nblyr+1) :: &
! Defined on Bio Grid interfaces
- iphin , & ! porosity
+ iphin , & ! porosity
ibrine_sal , & ! brine salinity (ppt)
ibrine_rho , & ! brine_density (kg/m^3)
iTin ! Temperature on the interface grid (oC)
- real (kind=dbl_kind) :: &
+ real (kind=dbl_kind) :: &
sloss ! brine flux contribution from surface runoff (g/m^2)
real (kind=dbl_kind), dimension (ncat) :: &
hbrnInitial, & ! inital brine height
- hbrnFinal ! category initial and final brine heights
+ hbrnFinal ! category initial and final brine heights
! for bgc sk
- real (kind=dbl_kind) :: &
+ real (kind=dbl_kind) :: &
dh_bot_chl , & ! Chlorophyll may or may not flush
dh_top_chl , & ! Chlorophyll may or may not flush
darcy_V_chl
@@ -5712,7 +5712,8 @@ subroutine colpkg_biogeochemistry(dt, &
congeln(n), snoicen(n), &
nbtrcr, fsnow, &
ntrcr, trcrn(1:ntrcr,n), &
- bio_index(1:nbtrcr), aicen_init(n), &
+ bio_index(1:nbtrcr), bio_index_o(:), &
+ aicen_init(n), &
vicen_init(n), vsnon_init(n), &
vicen(n), vsnon(n), &
aicen(n), flux_bio_atm(:), &
@@ -5791,8 +5792,8 @@ end subroutine colpkg_biogeochemistry
subroutine colpkg_init_hbrine(bgrid, igrid, cgrid, &
icgrid, swgrid, nblyr, nilyr, phi_snow)
- use ice_constants_colpkg, only: c1, c1p5, c2, p5, c0, rhoi, rhos
-
+ use ice_constants_colpkg, only: c1, c1p5, c2, p5, c0, rhoi, rhos, p25
+
integer (kind=int_kind), intent(in) :: &
nilyr, & ! number of ice layers
nblyr ! number of bio layers
@@ -5805,9 +5806,9 @@ subroutine colpkg_init_hbrine(bgrid, igrid, cgrid, &
real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
igrid ! biology vertical interface points
-
+
real (kind=dbl_kind), dimension (nilyr+1), intent(out) :: &
- cgrid , & ! CICE vertical coordinate
+ cgrid , & ! CICE vertical coordinate
icgrid , & ! interface grid for CICE (shortwave variable)
swgrid ! grid for ice tracers used in dEdd scheme
@@ -5815,65 +5816,65 @@ subroutine colpkg_init_hbrine(bgrid, igrid, cgrid, &
k , & ! vertical index
n ! thickness category index
- real (kind=dbl_kind) :: &
+ real (kind=dbl_kind) :: &
zspace ! grid spacing for CICE vertical grid
if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
!-----------------------------------------------------------------
- ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
+ ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
!-----------------------------------------------------------------
- bgrid(:) = c0 ! zsalinity grid points
+ bgrid(:) = c0 ! zsalinity grid points
bgrid(nblyr+2) = c1 ! bottom value
- igrid(:) = c0 ! bgc interface grid points
+ igrid(:) = c0 ! bgc interface grid points
igrid(1) = c0 ! ice top
igrid(nblyr+1) = c1 ! ice bottom
-
+
zspace = c1/max(c1,(real(nblyr,kind=dbl_kind)))
do k = 2, nblyr+1
bgrid(k) = zspace*(real(k,kind=dbl_kind) - c1p5)
enddo
-
+
do k = 2, nblyr
igrid(k) = p5*(bgrid(k+1)+bgrid(k))
enddo
!-----------------------------------------------------------------
- ! Calculate CICE cgrid for interpolation ice top (0) to bottom (1)
+ ! Calculate CICE cgrid for interpolation ice top (0) to bottom (1)
!-----------------------------------------------------------------
-
+
cgrid(1) = c0 ! CICE vertical grid top point
zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
-
+
do k = 2, nilyr+1
- cgrid(k) = zspace * (real(k,kind=dbl_kind) - c1p5)
- enddo
+ cgrid(k) = zspace * (real(k,kind=dbl_kind) - c1p5)
+ enddo
!-----------------------------------------------------------------
! Calculate CICE icgrid for ishortwave interpolation top(0) , bottom (1)
!-----------------------------------------------------------------
-
- icgrid(1) = c0
+
+ icgrid(1) = c0
zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
-
+
do k = 2, nilyr+1
icgrid(k) = zspace * (real(k,kind=dbl_kind)-c1)
- enddo
+ enddo
!------------------------------------------------------------------------
! Calculate CICE swgrid for dEdd ice: top of ice (0) , bottom of ice (1)
! Does not include snow
! see ice_shortwave.F90
! swgrid represents the layer index of the delta-eddington ice layer index
- !------------------------------------------------------------------------
+ !------------------------------------------------------------------------
zspace = c1/(real(nilyr,kind=dbl_kind)) ! CICE grid spacing
- swgrid(1) = min(c1/60.0_dbl_kind, zspace/c2)
+ swgrid(1) = min(c1/60.0_dbl_kind, zspace*p25)
swgrid(2) = zspace/c2 !+ swgrid(1)
do k = 3, nilyr+1
swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
- enddo
+ enddo
end subroutine colpkg_init_hbrine
diff --git a/column/ice_colpkg_shared.F90 b/column/ice_colpkg_shared.F90
index 76694ef290d9..cb8125cbd24b 100644
--- a/column/ice_colpkg_shared.F90
+++ b/column/ice_colpkg_shared.F90
@@ -89,7 +89,8 @@ module ice_colpkg_shared
real (kind=dbl_kind), parameter, public :: &
hi_ssl = 0.050_dbl_kind, & ! ice surface scattering layer thickness (m)
- hs_ssl = 0.040_dbl_kind ! snow surface scattering layer thickness (m)
+ hs_ssl = 0.040_dbl_kind, & ! snow surface scattering layer thickness (m)
+ hs_ssl_min = 5.0e-4_dbl_kind ! minimum snow scattering layer thickness for aerosol accumulation (m)
! snicar 5 band system, set in namelist
logical (kind=log_kind), public :: &
diff --git a/column/ice_itd.F90 b/column/ice_itd.F90
index 5438f7b1f05e..a9de8a2ba74f 100644
--- a/column/ice_itd.F90
+++ b/column/ice_itd.F90
@@ -1395,7 +1395,7 @@ subroutine zap_snow(dt, nslyr, &
endif ! tr_aero
if (z_tracers) then
- dvssl = min(p5*vsnon, hs_ssl*aicen) !snow surface layer
+ dvssl = min(p5*vsnon/real(nslyr,kind=dbl_kind), hs_ssl*aicen) !snow surface layer
dvint = vsnon- dvssl !snow interior
do it = 1, nbtrcr
diff --git a/column/ice_shortwave.F90 b/column/ice_shortwave.F90
index e8c12fde79ed..b4a0df15462e 100644
--- a/column/ice_shortwave.F90
+++ b/column/ice_shortwave.F90
@@ -1609,6 +1609,11 @@ subroutine shortwave_dEdd (n_aero, n_zaero, &
write(warning,*) ' aice = ', &
aice
call add_warning(warning)
+ do k = 0, klev
+ write(warning,*) ' zbio(nlt_zaero_sw(1)+k) = ', &
+ zbio(nlt_zaero_sw(1)+k)
+ call add_warning(warning)
+ end do
write(warning,*) ' hs = ', &
hs
call add_warning(warning)
@@ -2545,13 +2550,11 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, &
! aerosol in snow
- if (tr_zaero .and. dEdd_algae) then
+ if (tr_zaero .and. dEdd_algae) then
do k = 0,nslyr
- gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
- wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
- g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
- (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
- w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
+ g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)) / &
+ (w0(k)*tau(k) + wzaer(ns,k))
+ w0(k) = (w0(k)*tau(k) + wzaer(ns,k)) / &
(tau(k) + tzaer(ns,k))
tau(k) = tau(k) + tzaer(ns,k)
enddo
@@ -2725,12 +2728,10 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, &
g(k) = gi_int(ns)
! aerosol in sea ice
if (tr_zaero .and. dEdd_algae) then
- do k = kii, klev
- gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
- wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
- g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
- (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
- w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
+ do k = kii, klev
+ g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k))/ &
+ (w0(k)*tau(k) + wzaer(ns,k))
+ w0(k) = (w0(k)*tau(k) + wzaer(ns,k)) / &
(tau(k) + tzaer(ns,k))
tau(k) = tau(k) + tzaer(ns,k)
enddo
@@ -3353,6 +3354,9 @@ subroutine solution_dEdd &
smr , & ! accumulator for rdif gaussian integration
smt ! accumulator for tdif gaussian integration
+ character(len=char_len_long) :: &
+ warning ! warning message
+
! Delta-Eddington solution expressions
alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
@@ -3455,7 +3459,7 @@ subroutine solution_dEdd &
amg = alp - gam
rdir(k) = apg*rdif_a(k) + amg*(tdif_a(k)*trnlay(k) - c1)
tdir(k) = apg*tdif_a(k) + (amg* rdif_a(k)-apg+c1)*trnlay(k)
-
+
! recalculate rdif,tdif using direct angular integration over rdir,tdir,
! since Delta-Eddington rdif formula is not well-behaved (it is usually
! biased low and can even be negative); use ngmax angles and gaussian
@@ -3768,7 +3772,7 @@ subroutine compute_shortwave_trcr(n_algae, nslyr, &
sw_grid, hin, &
hbri, ntrcr, &
nilyr, nblyr, &
- i_grid, &
+ i_grid, aicen, &
nbtrcr_sw, n_zaero, &
skl_bgc, z_tracers, &
l_stop, stop_label)
@@ -3803,7 +3807,8 @@ subroutine compute_shortwave_trcr(n_algae, nslyr, &
real(kind=dbl_kind), intent(in) :: &
hin , & ! CICE ice thickness
- hbri ! brine height
+ hbri , & ! brine height
+ aicen
logical (kind=log_kind), intent(in) :: &
skl_bgc, & ! skeletal layer bgc
@@ -3839,17 +3844,18 @@ subroutine compute_shortwave_trcr(n_algae, nslyr, &
do k = 1,nilyr+1
icegrid(k) = sw_grid(k)
- enddo
- if (sw_grid(1)*hin*c2 > hi_ssl) then
+ enddo
+ if (sw_grid(1)*hin*c2 > hi_ssl .and. hin > puny) then
icegrid(1) = hi_ssl/c2/hin
endif
+ icegrid(2) = c2*sw_grid(1) + (sw_grid(2) - sw_grid(1)) ! sw_grid(2) = zpace/c2
if (z_tracers) then
if (tr_bgc_N) then
do k = 1, nblyr+1
do n = 1, n_algae
trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) + &
- R_chl2N(n)*F_abs_chl(n)*trcrn(nt_bgc_N(n)+k-1)
+ R_chl2N(n)*F_abs_chl(n)*trcrn(nt_bgc_N(n)+k-1) !*aicen
enddo ! n
enddo ! k
@@ -3888,7 +3894,7 @@ subroutine compute_shortwave_trcr(n_algae, nslyr, &
trtmp(:) = c0
do k = 1, nblyr+1
- trtmp0(nt_zaero(n) + k-1) = trcrn(nt_zaero(n)+k-1)
+ trtmp0(nt_zaero(n) + k-1) = trcrn(nt_zaero(n)+k-1) !*aicen
enddo
top_conc = trtmp0(nt_zaero(n))*min_bgc
@@ -4385,9 +4391,6 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting)
!chlorophyll mass extinction cross section (m^2/mg chla)
- character(len=char_len_long) :: &
- warning ! warning message
-
! SNICAR
! new inputs
integer (kind=int_kind), parameter :: &
@@ -4486,6 +4489,9 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
gi_int_mn_5bd = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind, &
.94_dbl_kind, .94_dbl_kind /)
+ character(len=char_len_long) :: &
+ warning ! warning message
+
!-----------------------------------------------------------------------
! Initialize and tune bare ice/ponded ice iops
@@ -4850,13 +4856,93 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
!aerosol in snow
if (tr_zaero .and. dEdd_algae) then
do k = 0,nslyr
- gzaer_5bd(ns,k) = gzaer_5bd(ns,k)/(wzaer_5bd(ns,k)+puny)
- wzaer_5bd(ns,k) = wzaer_5bd(ns,k)/(tzaer_5bd(ns,k)+puny)
- g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)*wzaer_5bd(ns,k)*tzaer_5bd(ns,k)) / &
- (w0(k)*tau(k) + wzaer_5bd(ns,k)*tzaer_5bd(ns,k))
- w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)*tzaer_5bd(ns,k)) / &
+
+ g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / &
+ (w0(k)*tau(k) + wzaer_5bd(ns,k))
+ w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
(tau(k) + tzaer_5bd(ns,k))
+
+ if (w0(k) > c1) then
+ write(warning,*) 'k = ',k
+ call add_warning(warning)
+ write(warning,*) 'w0(k) = ',w0(k)
+ call add_warning(warning)
+ write(warning,*) 'tau(k) = ',tau(k)
+ call add_warning(warning)
+ write(warning,*) 'g(k) = ',g(k)
+ call add_warning(warning)
+ write(warning,*) 'wzaer_5bd(ns,k) = ',wzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'tzaer_5bd(ns,k) = ',tzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'gzaer_5bd(ns,k) = ',gzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'k_bcexs(k)=', k_bcexs(k)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,1)=',kaer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,1)=',waer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,2)=',kaer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,2)=',waer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,3)=',kaer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,3)=',waer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(1)+k)=',zbio(nlt_zaero_sw(1)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(2)+k)=',zbio(nlt_zaero_sw(2)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(3)+k)=',zbio(nlt_zaero_sw(3)+k)
+ call add_warning(warning)
+ write(warning,*) 'dzk(k)=',dzk(k)
+ call add_warning(warning)
+ w0(k) = c1
+ endif
+
+ if (g(k) .ge. c1) then
+ write(warning,*) 'g(k) >= c1: k = ',k
+ call add_warning(warning)
+ write(warning,*) 'w0(k) = ',w0(k)
+ call add_warning(warning)
+ write(warning,*) 'tau(k) = ',tau(k)
+ call add_warning(warning)
+ write(warning,*) 'g(k) = ',g(k)
+ call add_warning(warning)
+ write(warning,*) 'wzaer_5bd(ns,k) = ',wzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'tzaer_5bd(ns,k) = ',tzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'gzaer_5bd(ns,k) = ',gzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'k_bcexs(k)=', k_bcexs(k)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,1)=',kaer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,1)=',waer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,2)=',kaer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,2)=',waer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,3)=',kaer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,3)=',waer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(1)+k)=',zbio(nlt_zaero_sw(1)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(2)+k)=',zbio(nlt_zaero_sw(2)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(3)+k)=',zbio(nlt_zaero_sw(3)+k)
+ call add_warning(warning)
+ write(warning,*) 'dzk(k)=',dzk(k)
+ call add_warning(warning)
+ endif
+
tau(k) = tau(k) + tzaer_5bd(ns,k)
+
enddo
elseif (tr_aero) then
k = 0 ! snow SSL
@@ -5017,12 +5103,90 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
! aerosol in sea ice
if (tr_zaero .and. dEdd_algae) then
do k = kii, klev
- gzaer_5bd(ns,k) = gzaer_5bd(ns,k)/(wzaer_5bd(ns,k)+puny)
- wzaer_5bd(ns,k) = wzaer_5bd(ns,k)/(tzaer_5bd(ns,k)+puny)
- g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)*wzaer_5bd(ns,k)*tzaer_5bd(ns,k)) / &
- (w0(k)*tau(k) + wzaer_5bd(ns,k)*tzaer_5bd(ns,k))
- w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)*tzaer_5bd(ns,k)) / &
+ g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / &
+ (w0(k)*tau(k) + wzaer_5bd(ns,k))
+ w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
(tau(k) + tzaer_5bd(ns,k))
+
+ if (w0(k) > c1) then
+ write(warning,*) 'k = ',k
+ call add_warning(warning)
+ write(warning,*) 'w0(k) = ',w0(k)
+ call add_warning(warning)
+ write(warning,*) 'tau(k) = ',tau(k)
+ call add_warning(warning)
+ write(warning,*) 'g(k) = ',g(k)
+ call add_warning(warning)
+ write(warning,*) 'wzaer_5bd(ns,k) = ',wzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'tzaer_5bd(ns,k) = ',tzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'gzaer_5bd(ns,k) = ',gzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'k_bcexs(k)=', k_bcexs(k)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,1)=',kaer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,1)=',waer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,2)=',kaer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,2)=',waer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,3)=',kaer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,3)=',waer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(1)+k)=',zbio(nlt_zaero_sw(1)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(2)+k)=',zbio(nlt_zaero_sw(2)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(3)+k)=',zbio(nlt_zaero_sw(3)+k)
+ call add_warning(warning)
+ write(warning,*) 'dzk(k)=',dzk(k)
+ call add_warning(warning)
+ w0(k) = c1
+ endif
+
+ if (g(k) .ge. c1) then
+ write(warning,*) 'g(k) >= c1: k = ',k
+ call add_warning(warning)
+ write(warning,*) 'w0(k) = ',w0(k)
+ call add_warning(warning)
+ write(warning,*) 'tau(k) = ',tau(k)
+ call add_warning(warning)
+ write(warning,*) 'g(k) = ',g(k)
+ call add_warning(warning)
+ write(warning,*) 'wzaer_5bd(ns,k) = ',wzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'tzaer_5bd(ns,k) = ',tzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'gzaer_5bd(ns,k) = ',gzaer_5bd(ns,k)
+ call add_warning(warning)
+ write(warning,*) 'k_bcexs(k)=', k_bcexs(k)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,1)=',kaer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,1)=',waer_bc_tab_5bd(ns,1)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,2)=',kaer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,2)=',waer_bc_tab_5bd(ns,2)
+ call add_warning(warning)
+ write(warning,*) 'kaer_bc_tab_5bd(ns,3)=',kaer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'waer_bc_tab_5bd(ns,3)=',waer_bc_tab_5bd(ns,3)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(1)+k)=',zbio(nlt_zaero_sw(1)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(2)+k)=',zbio(nlt_zaero_sw(2)+k)
+ call add_warning(warning)
+ write(warning,*) 'zbio(nlt_zaero_sw(3)+k)=',zbio(nlt_zaero_sw(3)+k)
+ call add_warning(warning)
+ write(warning,*) 'dzk(k)=',dzk(k)
+ call add_warning(warning)
+ endif
+
tau(k) = tau(k) + tzaer_5bd(ns,k)
enddo
elseif (tr_aero) then
@@ -5216,8 +5380,8 @@ subroutine compute_dEdd_5bd (nilyr, nslyr, klev, klevp, &
enddo
elseif (nsky == 2) then ! diffuse (keep the diffuse incident results)
do k = 0, klevp
- dfdif_snicar(k) = dfdif(k)
- rupdif_snicar(k) = rupdif(k)
+ dfdif_snicar(k) = dfdif(k)
+ rupdif_snicar(k) = rupdif(k)
enddo
endif
enddo ! end direct/diffuse incident nsky
diff --git a/column/ice_therm_itd.F90 b/column/ice_therm_itd.F90
index 162c1fa82e4b..0c8fac09281d 100644
--- a/column/ice_therm_itd.F90
+++ b/column/ice_therm_itd.F90
@@ -977,10 +977,10 @@ subroutine lateral_melt (dt, ncat, &
!-----------------------------------------------------------------
! Biogeochemistry
- !-----------------------------------------------------------------
+ !-----------------------------------------------------------------
if (z_tracers) then ! snow tracers
- dvssl = min(p5*vsnon(n), hs_ssl*aicen(n)) !snow surface layer
+ dvssl = min(p5*vsnon(n)/real(nslyr,kind=dbl_kind), hs_ssl*aicen(n)) !snow surface layer
dvint = vsnon(n)- dvssl !snow interior
do k = 1, nbtrcr
flux_bio(k) = flux_bio(k) &
diff --git a/shared/mpas_seaice_column.F b/shared/mpas_seaice_column.F
index 55a908c878ff..d5e3f5851f15 100644
--- a/shared/mpas_seaice_column.F
+++ b/shared/mpas_seaice_column.F
@@ -8531,13 +8531,15 @@ subroutine get_cice_biogeochemistry_tracer_array_cell(block, tracerObject, trace
verticalHumicsIceCell, &
verticalParticulateIronIceCell, &
verticalDissolvedIronIceCell, &
- verticalAerosolsIceCell
+ verticalAerosolsIceCell, &
+ verticalAerosolsSnowCell
integer :: &
iBioTracers, &
iBioCount, &
iLayers, &
- iIceCount
+ iIceCount, &
+ iSnowCount
call MPAS_pool_get_config(block % configs, "config_use_skeletal_biochemistry", config_use_skeletal_biochemistry)
call MPAS_pool_get_config(block % configs, "config_use_vertical_biochemistry", config_use_vertical_biochemistry)
@@ -8612,6 +8614,7 @@ subroutine get_cice_biogeochemistry_tracer_array_cell(block, tracerObject, trace
call MPAS_pool_get_array(tracers_aggregate, "verticalParticulateIronIceCell", verticalParticulateIronIceCell)
call MPAS_pool_get_array(tracers_aggregate, "verticalDissolvedIronIceCell", verticalDissolvedIronIceCell)
call MPAS_pool_get_array(tracers_aggregate, "verticalAerosolsIceCell", verticalAerosolsIceCell)
+ call MPAS_pool_get_array(tracers_aggregate, "verticalAerosolsSnowCell", verticalAerosolsSnowCell)
call MPAS_pool_get_array(tracers_aggregate, "verticalSalinityCell", verticalSalinityCell)
call MPAS_pool_get_array(tracers_aggregate, "brineFractionCell", brineFractionCell)
@@ -8890,6 +8893,7 @@ subroutine get_cice_biogeochemistry_tracer_array_cell(block, tracerObject, trace
iBioCount = 0
do iBioTracers = 1, nzAerosols
iIceCount = (iBioTracers-1)*nBioLayersP1
+ iSnowCount = (iBioTracers-1)*2
do iLayers = 1,nBioLayersP1
iBioCount = iBioCount + 1
@@ -8901,6 +8905,8 @@ subroutine get_cice_biogeochemistry_tracer_array_cell(block, tracerObject, trace
iBioCount = iBioCount + 1
verticalAerosolsConcCell(iBioCount,iCell) = &
tracerArrayCell(tracerObject % index_verticalAerosolsConc(iBioTracers)+iLayers-1)
+ verticalAerosolsSnowCell(iLayers-nBioLayersP1+iSnowCount,iCell) = &
+ verticalAerosolsConcCell(iBioCount,iCell)
enddo
enddo
endif