diff --git a/columnphysics/icepack_aerosol.F90 b/columnphysics/icepack_aerosol.F90 index 8e6a4933c..d0df24a88 100644 --- a/columnphysics/icepack_aerosol.F90 +++ b/columnphysics/icepack_aerosol.F90 @@ -8,8 +8,8 @@ module icepack_aerosol use icepack_kinds - use icepack_parameters, only: c0, c1, c2, puny, rhoi, rhos, hs_min - use icepack_parameters, only: hi_ssl, hs_ssl + use icepack_parameters, only: c0, c1, c2, p5, puny, rhoi, rhos, hs_min + use icepack_parameters, only: hi_ssl, hs_ssl, hs_ssl_min use icepack_tracers, only: max_aero use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted @@ -27,8 +27,8 @@ module icepack_aerosol !======================================================================= -! Increase aerosol in ice or snow surface due to deposition -! and vertical cycling +! Deposition and vertical cycling of aerosol in ice or snow +! Called from icepack_step_therm1 when tr_aero=T (not used for zbgc tracers) subroutine update_aerosol(dt, & nilyr, nslyr, & @@ -128,6 +128,7 @@ subroutine update_aerosol(dt, & hs = vsnon*ar hi = vicen*ar + ! fluxes were divided by aice for thermo, not yet multiplied by aice dhs_melts = -melts dhi_snoice = snoice dhs_snoice = dhi_snoice*rhoi/rhos @@ -424,8 +425,8 @@ end subroutine update_aerosol !======================================================================= -! Increase aerosol in snow surface due to deposition -! and vertical cycling : after update_aerosol +! Aerosol in snow for vertical biogeochemistry with mushy thermodynamics +! Called from icepack_algae.F90 when z_tracers=T (replaces update_aerosol) subroutine update_snow_bgc (dt, nblyr, & nslyr, & @@ -439,6 +440,9 @@ subroutine update_snow_bgc (dt, nblyr, & vicen, vsnon, & aicen, flux_bio_atm,& zbgc_atm, flux_bio) +!echmod: add bio_index_o when updating ice_algae.F90 +! zbgc_atm, flux_bio, & +! bio_index_o) integer (kind=int_kind), intent(in) :: & nbtrcr, & ! number of distinct snow tracers @@ -448,6 +452,9 @@ subroutine update_snow_bgc (dt, nblyr, & integer (kind=int_kind), dimension (nbtrcr), intent(in) :: & bio_index +!echmod: add bio_index_o when updating ice_algae.F90 +! bio_index, & +! bio_index_o ! provides index of scavenging (kscavz) data array real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -466,9 +473,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), intent(in) :: & @@ -484,6 +493,9 @@ subroutine update_snow_bgc (dt, nblyr, & real (kind=dbl_kind) :: & dzssl, dzssl_new, & ! snow ssl thickness dzint, dzint_new, & ! snow interior thickness + dz, & ! + 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 @@ -506,6 +518,13 @@ subroutine update_snow_bgc (dt, nblyr, & character(len=*),parameter :: subname='(update_snow_bgc)' +!echmod: add bio_index_o to subroutine arguments when updating ice_algae.F90 + integer (kind=int_kind), dimension (nbtrcr) :: & + bio_index_o ! provides index of scavenging (kscavz) data array + do k = 1, nbtrcr + bio_index_o(k) = k ! TEMPORARY + enddo + !------------------------------------------------------------------- ! initialize !------------------------------------------------------------------- @@ -516,6 +535,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) @@ -524,30 +548,36 @@ 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 + 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+1) = c0 trcrn(bio_index(k) + nblyr+2) = c0 - zbgc_atm(k) = zbgc_atm(k) & + 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 @@ -565,6 +595,48 @@ 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 @@ -573,38 +645,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 @@ -613,39 +684,41 @@ 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 + ! 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 + 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 @@ -665,30 +738,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 @@ -701,6 +775,11 @@ subroutine update_snow_bgc (dt, nblyr, & aero_cons(k) = aerotot(k)-aerotot0(k) & - ( flux_bio_atm(k) & - (flux_bio(k)-flux_bio_o(k))) * dt + 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(warnstr,*) subname, 'Conservation failure: aerosols in snow' call icepack_warnings_add(warnstr) @@ -729,18 +808,18 @@ 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+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 endif + !------------------------------------------------------------------- ! check for negative values !------------------------------------------------------------------- diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index d59b38319..40258e7e9 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -186,6 +186,7 @@ module icepack_parameters kappav = 1.4_dbl_kind ,&! vis extnctn coef in ice, wvlngth<700nm (1/m) 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_min = 5.e-4_dbl_kind,&! minimum snow surface scattering layer thickness for aerosol (m) ! baseline albedos for ccsm3 shortwave, set in namelist albicev = 0.78_dbl_kind ,&! visible ice albedo for h > ahmax albicei = 0.36_dbl_kind ,&! near-ir ice albedo for h > ahmax @@ -449,7 +450,7 @@ subroutine icepack_init_parameters( & zref_in, hs_min_in, snowpatch_in, rhosi_in, sk_l_in, & saltmax_in, phi_init_in, min_salin_in, salt_loss_in, & Tliquidus_max_in, & - min_bgc_in, dSin0_frazil_in, hi_ssl_in, hs_ssl_in, & + min_bgc_in, dSin0_frazil_in, hi_ssl_in, hs_ssl_in, hs_ssl_min_in, & awtvdr_in, awtidr_in, awtvdf_in, awtidf_in, & qqqice_in, TTTice_in, qqqocn_in, TTTocn_in, & ktherm_in, conduct_in, fbot_xfer_type_in, calc_Tsfc_in, dts_b_in, & @@ -594,7 +595,8 @@ subroutine icepack_init_parameters( & stefan_boltzmann_in, & ! W/m^2/K^4 kappav_in, & ! vis extnctn coef in ice, wvlngth<700nm (1/m) hi_ssl_in, & ! ice surface scattering layer thickness (m) - hs_ssl_in, & ! visible, direct + hs_ssl_in, & ! snow surface scattering layer thickness (m) + hs_ssl_min_in, & ! minimum snow surface scattering layer thickness for aerosols (m) awtvdr_in, & ! visible, direct ! for history and awtidr_in, & ! near IR, direct ! diagnostics awtvdf_in, & ! visible, diffuse @@ -899,6 +901,7 @@ subroutine icepack_init_parameters( & if (present(dSin0_frazil_in) ) dSin0_frazil = dSin0_frazil_in if (present(hi_ssl_in) ) hi_ssl = hi_ssl_in if (present(hs_ssl_in) ) hs_ssl = hs_ssl_in + if (present(hs_ssl_min_in) ) hs_ssl_min = hs_ssl_min_in if (present(awtvdr_in) ) awtvdr = awtvdr_in if (present(awtidr_in) ) awtidr = awtidr_in if (present(awtvdf_in) ) awtvdf = awtvdf_in @@ -1174,7 +1177,7 @@ subroutine icepack_query_parameters( & zref_out, hs_min_out, snowpatch_out, rhosi_out, sk_l_out, & saltmax_out, phi_init_out, min_salin_out, salt_loss_out, & Tliquidus_max_out, & - min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, & + min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, hs_ssl_min_out, & awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, cpl_frazil_out, & qqqice_out, TTTice_out, qqqocn_out, TTTocn_out, update_ocn_f_out, & Lfresh_out, cprho_out, Cp_out, ustar_min_out, hi_min_out, a_rapid_mode_out, & @@ -1329,7 +1332,8 @@ subroutine icepack_query_parameters( & stefan_boltzmann_out, & ! W/m^2/K^4 kappav_out, & ! vis extnctn coef in ice, wvlngth<700nm (1/m) hi_ssl_out, & ! ice surface scattering layer thickness (m) - hs_ssl_out, & ! visible, direct + hs_ssl_out, & ! snow surface scattering layer thickness (m) + hs_ssl_min_out, & ! minimum snow surface scattering layer thickness for aerosols (m) awtvdr_out, & ! visible, direct ! for history and awtidr_out, & ! near IR, direct ! diagnostics awtvdf_out, & ! visible, diffuse @@ -1666,6 +1670,7 @@ subroutine icepack_query_parameters( & if (present(dSin0_frazil_out) ) dSin0_frazil_out = dSin0_frazil if (present(hi_ssl_out) ) hi_ssl_out = hi_ssl if (present(hs_ssl_out) ) hs_ssl_out = hs_ssl + if (present(hs_ssl_min_out) ) hs_ssl_min_out = hs_ssl_min if (present(awtvdr_out) ) awtvdr_out = awtvdr if (present(awtidr_out) ) awtidr_out = awtidr if (present(awtvdf_out) ) awtvdf_out = awtvdf @@ -1864,6 +1869,7 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " dSin0_frazil = ",dSin0_frazil write(iounit,*) " hi_ssl = ",hi_ssl write(iounit,*) " hs_ssl = ",hs_ssl + write(iounit,*) " hs_ssl_min = ",hs_ssl_min write(iounit,*) " awtvdr = ",awtvdr write(iounit,*) " awtidr = ",awtidr write(iounit,*) " awtvdf = ",awtvdf diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index f6c1a7223..b8150fe3c 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -3497,9 +3497,10 @@ subroutine compute_shortwave_trcr(bgcN, zaero, & do k = 1,nilyr+1 icegrid(k) = sw_grid(k) enddo - if (sw_grid(1)*hin*c2 > hi_ssl) then + 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)) if (z_tracers) then if (tr_bgc_N) then