diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 9ce598c99..30c4bd98d 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -61,12 +61,17 @@ subroutine atmo_boundary_layer (sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, tr_iso, & + Qa_iso, Qref_iso, & uvel, vvel, & Uref ) character (len=3), intent(in) :: & sfctype ! ice or ocean + logical (kind=log_kind), optional, intent(in) :: & + tr_iso + logical (kind=log_kind), intent(in) :: & calc_strair, & ! if true, calculate wind stress components formdrag, & ! if true, calculate form drag @@ -75,6 +80,9 @@ subroutine atmo_boundary_layer (sfctype, & integer (kind=int_kind), intent(in) :: & natmiter ! number of iterations for boundary layer calculations + integer (kind=int_kind), optional, intent(in) :: & + n_iso ! number of isotopes + real (kind=dbl_kind), intent(in) :: & Tsf , & ! surface temperature of ice or ocean potT , & ! air potential temperature (K) @@ -104,6 +112,12 @@ subroutine atmo_boundary_layer (sfctype, & shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat + real (kind=dbl_kind), intent(in), optional, dimension(:) :: & + Qa_iso ! specific isotopic humidity (kg/kg) + + real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & + Qref_iso ! reference specific isotopic humidity (kg/kg) + real (kind=dbl_kind), intent(in) :: & uvel , & ! x-direction ice speed (m/s) vvel ! y-direction ice speed (m/s) @@ -114,7 +128,7 @@ subroutine atmo_boundary_layer (sfctype, & ! local variables integer (kind=int_kind) :: & - k ! iteration index + k,n ! iteration index real (kind=dbl_kind) :: & TsfK , & ! surface temperature in Kelvin (K) @@ -136,6 +150,7 @@ subroutine atmo_boundary_layer (sfctype, & ustar , & ! ustar (m/s) tstar , & ! tstar qstar , & ! qstar + ratio , & ! ratio rdn , & ! sqrt of neutral exchange coefficient (momentum) rhn , & ! sqrt of neutral exchange coefficient (heat) ren , & ! sqrt of neutral exchange coefficient (water) @@ -356,6 +371,19 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif + if (present(tr_iso) .and. present(Qref_iso) .and. & + present(Qa_iso) .and. present(n_iso)) then + if (tr_iso) then + Qref_iso(:) = c0 + do n = 1, n_iso + ratio = c0 + if (Qa_iso(2) > puny) & + ratio = Qa_iso(n)/Qa_iso(2) + Qref_iso(n) = Qa_iso(n) - ratio*delq*fac + enddo + endif + endif + end subroutine atmo_boundary_layer !======================================================================= @@ -800,7 +828,7 @@ end subroutine neutral_drag_coeffs !autodocument_start icepack_atm_boundary ! - subroutine icepack_atm_boundary(sfctype, & + subroutine icepack_atm_boundary(sfctype, & Tsf, potT, & uatm, vatm, & wind, zlvl, & @@ -811,6 +839,8 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, tr_iso, & + Qa_iso, Qref_iso, & uvel, vvel, & Uref) @@ -827,6 +857,12 @@ subroutine icepack_atm_boundary(sfctype, & Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) + integer (kind=int_kind), optional, intent(in) :: & + n_iso ! number of isotopes + + logical (kind=log_kind), optional, intent(in) :: & + tr_iso + real (kind=dbl_kind), intent(inout) :: & Cdn_atm , & ! neutral drag coefficient Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff @@ -844,6 +880,12 @@ subroutine icepack_atm_boundary(sfctype, & shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat + real (kind=dbl_kind), intent(in), optional, dimension(:) :: & + Qa_iso ! specific isotopic humidity (kg/kg) + + real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & + Qref_iso ! reference specific isotopic humidity (kg/kg) + real (kind=dbl_kind), optional, intent(in) :: & uvel , & ! x-direction ice speed (m/s) vvel ! y-direction ice speed (m/s) @@ -858,6 +900,15 @@ subroutine icepack_atm_boundary(sfctype, & real (kind=dbl_kind) :: & worku, workv, workr + integer (kind=int_kind) :: & + niso + + logical (kind=log_kind) :: & + triso + + real (kind=dbl_kind), dimension(:), allocatable :: & + Qaiso, Qriso + character(len=*),parameter :: subname='(icepack_atm_boundary)' worku = c0 @@ -869,6 +920,16 @@ subroutine icepack_atm_boundary(sfctype, & if (present(vvel)) then workv = vvel endif + if (present(n_iso) .and. present(Qa_iso) .and. & + present(Qref_iso) .and. present(tr_iso)) then + allocate(Qaiso(n_iso),Qriso(n_iso)) + Qaiso = Qa_iso + Qriso = Qref_iso + niso = n_iso + triso = tr_iso + else + allocate(Qaiso(1),Qriso(1)) + endif Cdn_atm_ratio_n = c1 @@ -896,8 +957,10 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - worku, workv, & - workr) + n_iso=niso, tr_iso=triso, & + Qa_iso=Qaiso, Qref_iso=Qriso,& + uvel=worku, vvel=workv, & + Uref=workr) if (icepack_warnings_aborted(subname)) return endif ! atmbndy @@ -905,6 +968,12 @@ subroutine icepack_atm_boundary(sfctype, & Uref = workr endif + if (present(Qref_iso)) then + Qref_iso = Qriso + endif + + deallocate(Qaiso,Qriso) + end subroutine icepack_atm_boundary !------------------------------------------------------------ diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index dc1cfb395..fa1730562 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -55,7 +55,11 @@ subroutine merge_fluxes (aicen, & meltt, melts, & meltb, & congel, snoice, & - Uref, Urefn ) + Uref, Urefn, & + Qref_iso, Qrefn_iso, & + fiso_ocn, fiso_ocnn, & + fiso_evap, fiso_evapn ) + ! single category fluxes real (kind=dbl_kind), intent(in) :: & @@ -119,6 +123,16 @@ subroutine merge_fluxes (aicen, & real (kind=dbl_kind), optional, intent(inout):: & Uref ! air speed reference level (m/s) + real (kind=dbl_kind), optional, dimension(:), intent(inout):: & + Qref_iso, & ! isotope air sp hum reference level (kg/kg) + fiso_ocn, & ! isotope fluxes to ocean (kg/m2/s) + fiso_evap ! isotope evaporation (kg/m2/s) + + real (kind=dbl_kind), optional, dimension(:), intent(in):: & + Qrefn_iso, & ! isotope air sp hum reference level (kg/kg) + fiso_ocnn, & ! isotope fluxes to ocean (kg/m2/s) + fiso_evapn ! isotope evaporation (kg/m2/s) + character(len=*),parameter :: subname='(merge_fluxes)' !----------------------------------------------------------------- @@ -148,6 +162,17 @@ subroutine merge_fluxes (aicen, & Tref = Tref + Trefn * aicen Qref = Qref + Qrefn * aicen + ! Isotopes + if (present(Qrefn_iso) .and. present(Qref_iso)) then + Qref_iso(:) = Qref_iso(:) + Qrefn_iso(:) * aicen + endif + if (present(fiso_ocnn) .and. present(fiso_ocn)) then + fiso_ocn(:) = fiso_ocn(:) + fiso_ocnn(:) * aicen + endif + if (present(fiso_evapn) .and. present(fiso_evap)) then + fiso_evap(:) = fiso_evap(:) + fiso_evapn(:) * aicen + endif + ! ocean fluxes if (present(Urefn) .and. present(Uref)) then Uref = Uref + Urefn * aicen diff --git a/columnphysics/icepack_intfc.F90 b/columnphysics/icepack_intfc.F90 index c457a8b09..f6ce6f494 100644 --- a/columnphysics/icepack_intfc.F90 +++ b/columnphysics/icepack_intfc.F90 @@ -42,6 +42,7 @@ module icepack_intfc use icepack_tracers, only: icepack_max_don => max_don use icepack_tracers, only: icepack_max_fe => max_fe use icepack_tracers, only: icepack_max_aero => max_aero + use icepack_tracers, only: icepack_max_iso => max_iso use icepack_tracers, only: icepack_nmodal1 => nmodal1 use icepack_tracers, only: icepack_nmodal2 => nmodal2 use icepack_parameters, only: icepack_nspint => nspint diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 new file mode 100644 index 000000000..c2b30fa5e --- /dev/null +++ b/columnphysics/icepack_isotope.F90 @@ -0,0 +1,680 @@ +!======================================================================= +! +!BOP +! +! !MODULE: icepack_isotope - Isotope tracer within sea ice +! +! !DESCRIPTION: +! +! !REVISION HISTORY: +! SVN:$$ +! +! authors: David Bailey, NCAR +! Jiang Zhu, UW-Madison +! +! 2014: Added i2x evaporation flux +! Added fractionation options +! Fixed bugs +! +! !INTERFACE: +! + module icepack_isotope +! +! !USES: +! + use icepack_kinds + use icepack_parameters, only: c0, c1, c2, p001, p1, p5, puny +! +!EOP +! + implicit none + + character(len=5), parameter :: & + frac = 'cfrac' ! fractionation coefficient calculation method + ! cfrac, constant fractionation + ! nfrac, nonfractionation + ! gfrac, growth-rate dependent for H2_18O + + ! Species indicies - public so thay can be seen by water_tracers + integer, parameter, public :: ispundef = 0 ! Undefined + integer, parameter, public :: isph2o = 1 ! H2O ! "regular" water + integer, parameter, public :: isph216o = 2 ! H216O ! H216O, nearly the same as "regular" water + integer, parameter, public :: isphdo = 3 ! HDO + integer, parameter, public :: isph218o = 4 ! H218O + + integer, parameter, public :: pwtspec = 4 ! number of water species (h2o,hdo,h218o,h216o) + +!======================================================================= + + contains + +!======================================================================= + +!BOP +! +! !ROUTINE: update_isotope +! +! !DESCRIPTION: +! +! Increase isotope in ice or snow surface due to deposition +! +! !REVISION HISTORY: same as module +! +! !INTERFACE: +! + subroutine update_isotope (dt, & + nilyr, nslyr, & + n_iso, & + meltt, melts, & + meltb, congel, & + snoice, evap, & + fsnow, Tsfc, & + Qref_iso, & + isosno, isoice, & + aice_old, & + vice_old, vsno_old, & + vicen, vsnon, aicen, & + fiso_atm, fiso_evapn, & + fiso_ocnn, HDO_ocn, & + H2_16O_ocn, H2_18O_ocn) +! +! !USES: +! +! use water_isotopes, only: wiso_alpi + use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh + use icepack_tracers, only: nt_iso +! +! !INPUT/OUTPUT PARAMETERS: +! + integer (kind=int_kind), intent(in) :: & + nilyr, nslyr, n_iso + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + real (kind=dbl_kind), intent(in) :: & + Tsfc, & + meltt, & + melts, & + meltb, & + congel, & ! congelation ice growth (m/step) + snoice, & ! ice thickness increase (m/step) + evap, & ! surface evaporation + fsnow, & ! snowfall (kg/m^2/s of water) + vicen, & ! volume per unit area of ice (m) + vsnon, & + aicen, & + aice_old, & + vice_old, & + vsno_old, & + HDO_ocn, & + H2_16O_ocn, & + H2_18O_ocn + + real (kind=dbl_kind), dimension(:), intent(in) :: & + fiso_atm, & ! isotopic snowfall (kg/m^2/s of water) + Qref_iso ! isotope reference humidity + + real (kind=dbl_kind), dimension(:), intent(inout) :: & + fiso_ocnn, & ! isotopic freshwater (kg/m^2/s) + fiso_evapn ! evaporative water flux (kg/m^2/s) + + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & + isosno, isoice ! mass of isotopes (kg) + +! +! local variables +! + real (kind=dbl_kind) :: & + dzssl, & + dzint, & + dzssli, & + dzinti, & + evaps, & ! evaporation over snow (m/step) + evapi, & ! evaporation over ice (m/step) + dhs_snoice,& ! snow thickness reduction (m/step) + hi, & + hs + + real (kind=dbl_kind), dimension(n_iso) :: & + isotot, isotot0 ! for diagnostics + + real (kind=dbl_kind) :: & + dzssl_new, & + dzint_new, & + dzssli_new, & + dzinti_new, & + dznew + + real (kind=dbl_kind) :: & + hs_old, hi_old, hslyr_old, hilyr_old, dhs, dhi, & + hslyr, hilyr, sloss1, sloss2, & + TsfK, & ! snow/ice temperature (K) + alphai, & ! ice/vapour fractionation coefficient + ratio, & ! isotopic ratio + work, & ! temporary variable + alpha + + integer (kind=int_kind) :: k + +! These need to be the same as in the DE code. Put in a common place? + real (kind=dbl_kind), parameter :: & + hi_ssl = .050_dbl_kind, & + hs_ssl = .040_dbl_kind + +! initialize + + hs_old=vsno_old/aice_old + hi_old=vice_old/aice_old + hslyr_old=hs_old/real(nslyr,kind=dbl_kind) + hilyr_old=hi_old/real(nilyr,kind=dbl_kind) + + dzssl=min(hslyr_old/c2,hs_ssl) + dzint=hs_old-dzssl + dzssli=min(hilyr_old/c2,hi_ssl) + dzinti=hi_old-dzssli + + if (aicen > puny) then + hs = vsnon/aicen + hi = vicen/aicen + elseif (aice_old > puny) then + hs = vsnon/aice_old + hi = vicen/aice_old + endif + + if (ktherm == 2) then + dhs_snoice = snoice + else + dhs_snoice = snoice*rhoi/rhos + endif + +! if (hs > puny) then +! evaps = evap*dt/rhos +! else +! evapi = evap*dt/rhoi +! endif + evaps = hs-(hs_old-melts-dhs_snoice+& + fsnow/rhos*dt) + evapi = hi-(hi_old-meltt-meltb+congel+snoice) + +! condensation of vapor onto snow and ice + + TsfK = Tsfc + Tffresh + + if (evaps > c0) then ! condensation to snow + do k = 1, n_iso + ratio = c1 ! ratio between 18O(HDO) and 16O in humidity + alphai = c1 ! fractionation coefficient + if (frac.ne.'nfrac' .and. Qref_iso(2)>puny) & + ratio = Qref_iso(k)/Qref_iso(2) + if (Qref_iso(2) <= puny) ratio = c0 + if (frac.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK) + if (frac.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK) + if (frac.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK) + work = alphai*ratio*rhos*evaps*aicen + fiso_evapn(k) = fiso_evapn(k)+work/dt + isosno(k,1) = isosno(k,1)+work + enddo + dzssl = dzssl+evaps + endif + + if (evapi > c0) then ! condensation to ice + do k = 1, n_iso + ratio = c1 ! ratio between 18O(HDO) and 16O in ref humidity + alphai = c1 ! fractionation coefficient + if (frac.ne.'nfrac' .and. Qref_iso(2)>puny) & + ratio = Qref_iso(k)/Qref_iso(2) + if (Qref_iso(2) <= puny) ratio = c0 + if (frac.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK) + if (frac.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK) + if (frac.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK) + work = alphai*ratio*rhoi*evapi*aicen + fiso_evapn(k) = fiso_evapn(k)+work/dt + isoice(k,1) = isoice(k,1)+work + enddo + dzssli = dzssli+evapi + endif + +! basal ice growth and isotope uptake + + if (congel > c0) then + alpha = isoice_alpha(congel/dt,'HDO',frac) + work = alpha*HDO_ocn*rhoi*congel*aicen + isoice(1,2) = isoice(1,2)+work + fiso_ocnn(1) = fiso_ocnn(1)-work/dt + + alpha = isoice_alpha(congel/dt,'H2_16O',frac) + work = alpha*H2_16O_ocn*rhoi*congel*aicen + isoice(2,2) = isoice(2,2)+work + fiso_ocnn(2) = fiso_ocnn(2)-work/dt + + alpha = isoice_alpha(congel/dt,'H2_18O',frac) + work = alpha*H2_18O_ocn*rhoi*congel*aicen + isoice(3,2) = isoice(3,2)+work + fiso_ocnn(3) = fiso_ocnn(3)-work/dt + + dzinti = dzinti+congel + endif + +! sublimation of snow and ice + + if (evaps < c0) then ! snow sublimation (no fractionation) + do k = 1, n_iso + !ratio = c1 ! ratio between 18O(HDO) and 16O in snow ssl + !if (isosno(2,1) > puny) & + ! ratio = isosno(k,1)/isosno(2,1) + !if (ratio > c5) ratio = c1 !! remove latter? + !work = ratio*rhos*evaps*aicen + !fiso_evapn(k) = fiso_evapn(k)+work/dt + + sloss1 = c0 + sloss2 = c0 + if (dzssl > puny) & + sloss1 = isosno(k,1)* & + min(-evaps,dzssl)/dzssl + if (isosno(k,1) >= sloss1) then + isosno(k,1) = isosno(k,1)-sloss1 + else + sloss1 = isosno(k,1) + isosno(k,1) = c0 + endif + if (dzint > puny) & + sloss2 = isosno(k,2)* & + max(-evaps-dzssl,c0)/dzint + if (isosno(k,2) >= sloss2) then + isosno(k,2) = isosno(k,2)-sloss2 + else + sloss2 = isosno(k,2) + isosno(k,2) = c0 + endif +! if (isosno(k,2) < c0) then +! write(nu_diag,*) 'Neg isosno(k,2)',isosno(k,2),sloss2 +! endif + fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt + enddo + + dzint = dzint+min(dzssl+evaps,c0) + dzssl = max(dzssl+evaps,c0) + if ( dzssl <= c0) then ! ssl goes away + fiso_evapn(:) = fiso_evapn(:)-isosno(:,1)/dt + isosno(:,1) = c0 + dzssl = c0 + endif + if (dzint <= c0) then ! int goes away + fiso_evapn(:) = fiso_evapn(:)-isosno(:,2)/dt + isosno(:,2) = c0 + dzint = c0 + endif + endif + + if (evapi < c0) then ! ice sublimation (no fractionation) + do k = 1, n_iso + !!ratio = c1 ! ratio between 18O(HDO) and 16O in ice ssl + !!if (isoice(2,1) > puny) & + !! ratio = isoice(k,1)/isoice(2,1) + !!if (ratio > c5) ratio = c1 ! remove latter? + !!work = ratio*rhoi*evapi*aicen + !!fiso_evapn(k) = fiso_evapn(k)+work/dt + + sloss1 = c0 + sloss2 = c0 + if (dzssli > puny) & + sloss1 = isoice(k,1)* & + min(-evapi,dzssli)/dzssli + if (isoice(k,1) >= sloss1) then + isoice(k,1) = isoice(k,1)-sloss1 + else + sloss1 = isoice(k,1) + isoice(k,1) = c0 + endif + if (dzinti > puny) & + sloss2 = isoice(k,2)* & + max(-evapi-dzssli,c0)/dzinti + if (isoice(k,2) >= sloss2) then + isoice(k,2) = isoice(k,2)-sloss2 + else + sloss2 = isoice(k,2) + isoice(k,2) = c0 + endif + fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt + enddo + + dzinti = dzinti+min(dzssli+evapi,c0) + dzssli = max(dzssli+evapi,c0) + if ( dzssli <= c0) then ! ssl goes away + fiso_evapn(:) = fiso_evapn(:)-isoice(:,1)/dt + isoice(:,1) = c0 + dzssli = c0 + endif + if (dzinti <= c0) then ! int goes away + fiso_evapn(:) = fiso_evapn(:)-isoice(:,2)/dt + isoice(:,2) = c0 + dzinti = c0 + endif + endif + +! surface snow melt + + if (melts > c0) then + do k=1,n_iso + sloss1=c0 + sloss2=c0 + if (dzssl > puny) & + sloss1 = isosno(k,1) & + *min(melts,dzssl)/dzssl + if (isosno(k,1) >= sloss1) then + isosno(k,1) = isosno(k,1)-sloss1 + else + sloss1 = isosno(k,1) + isosno(k,1) = c0 + endif + if (dzint > puny) & + sloss2=isosno(k,2) & + *max(melts-dzssl,c0)/dzint + if (isosno(k,2) >= sloss2) then + isosno(k,2) = isosno(k,2)-sloss2 + else + sloss2 = isosno(k,2) + isosno(k,2) = c0 + endif +! if (isosno(k,2) < c0) then +! write(nu_diag,*) 'Neg isosno(k,2)',isosno(k,2),sloss2 +! endif + fiso_ocnn(k)=fiso_ocnn(k)+(sloss1+sloss2)/dt + enddo ! n_iso + + dzint=dzint+min(dzssl-melts,c0) + dzssl=max(dzssl-melts,c0) + if ( dzssl <= c0) then ! ssl melts away + fiso_ocnn(:)= fiso_ocnn(:)+isosno(:,1)/dt + isosno(:,1) = c0 + dzssl = c0 + endif + if (dzint <= c0) then ! int melts away + fiso_ocnn(:) = fiso_ocnn(:)+isosno(:,2)/dt + isosno(:,2) = c0 + dzint = c0 + endif + endif + +! surface ice melt + if (meltt > c0) then + do k=1,n_iso + sloss1=c0 + sloss2=c0 + if (dzssli > puny) & + sloss1=isoice(k,1) & + *min(meltt,dzssli)/dzssli + if (isoice(k,1) >= sloss1) then + isoice(k,1) = isoice(k,1)-sloss1 + else + sloss1 = isoice(k,1) + isoice(k,1) = c0 + endif + if (dzinti > puny) & + sloss2=isoice(k,2) & + *max(meltt-dzssli,c0)/dzinti + if (isoice(k,2) >= sloss2) then + isoice(k,2) = isoice(k,2)-sloss2 + else + sloss2 = isoice(k,2) + isoice(k,2) = c0 + endif + fiso_ocnn(k)=fiso_ocnn(k)+(sloss1+sloss2)/dt + enddo + + dzinti=dzinti+min(dzssli-meltt,c0) + dzssli=max(dzssli-meltt,c0) + if (dzssli <= c0) then ! ssl ice melts away + fiso_ocnn(:) = fiso_ocnn(:)+isoice(:,1) + isoice(:,1) = c0 + dzssli = c0 + endif + if (dzinti <= c0) then ! int ice melts away + fiso_ocnn(:) = fiso_ocnn(:)+isoice(:,2) + isoice(:,2) = c0 + dzinti = c0 + endif + endif + +! basal ice melt. Assume all isotopes lost in basal melt + + if (meltb > c0) then + do k=1,n_iso + sloss1=c0 + sloss2=c0 + if (dzssli > puny) & + sloss1=max(meltb-dzinti,c0) & + *isoice(k,1)/dzssli + if (isoice(k,1) >= sloss1) then + isoice(k,1) = isoice(k,1)-sloss1 + else + sloss1 = isoice(k,1) + isoice(k,1) = c0 + endif + if (dzinti > puny) & + sloss2=min(meltb,dzinti) & + *isoice(k,2)/dzinti + if (isoice(k,2) >= sloss2) then + isoice(k,2) = isoice(k,2)-sloss2 + else + sloss2 = isoice(k,2) + isoice(k,2) = c0 + endif + fiso_ocnn(k)=fiso_ocnn(k)+(sloss1+sloss2)/dt + enddo + + dzssli = dzssli+min(dzinti-meltb, c0) + dzinti = max(dzinti-meltb, c0) + if (dzssli <= c0) then ! ssl ice melts away + fiso_ocnn(:) = fiso_ocnn(:)+isoice(:,1) + isoice(:,1) = c0 + dzssli = c0 + endif + if (dzinti <= c0) then ! int ice melts away + fiso_ocnn(:) = fiso_ocnn(:)+isoice(:,2) + isoice(:,2) = c0 + dzinti = c0 + endif + endif + +! snowfall and isotope deposition + + if (fsnow > c0) then + isosno(:,1) = isosno(:,1)+ & + fiso_atm(:)*aicen*dt + dzssl = dzssl+fsnow/rhos*dt + endif + +! snoice formation + + if (dhs_snoice > c0) then + do k=1,n_iso + sloss1=c0 + sloss2=c0 + if (dzint > puny) & + sloss2 = min(dhs_snoice,dzint) & + *isosno(k,2)/dzint + if (isosno(k,2) >= sloss2) then + isosno(k,2) = isosno(k,2)-sloss2 + else + sloss2 = isosno(k,2) + isosno(k,2) = c0 + endif +! if (isosno(k,2) < c0) then +! write(nu_diag,*) 'Snow-ice isosno(k,2)',isosno(k,2),sloss2 +! endif + if (dzssl > puny) & + sloss1 = max(dhs_snoice-dzint,c0) & + *isosno(k,1)/dzssl + if (isosno(k,1) >= sloss1) then + isosno(k,1) = isosno(k,1)-sloss1 + else + sloss1 = isosno(k,1) + isosno(k,1) = c0 + endif + isoice(k,1) = isoice(k,1)+sloss1+sloss2 + enddo + + dzssl=dzssl-max(dhs_snoice-dzint,c0) + dzint=max(dzint-dhs_snoice,c0) + dzssli=dzssli+snoice + if ( dzssl <= c0) then ! ssl goes away + fiso_ocnn(:)= fiso_ocnn(:)+isosno(:,1)/dt + isosno(:,1) = c0 + dzssl = c0 + endif + if (dzint <= c0) then ! int goes away + fiso_ocnn(:) = fiso_ocnn(:)+isosno(:,2)/dt + isosno(:,2) = c0 + dzint = c0 + endif + endif + +! redistribute isotope within vertical layers + + hslyr = hs/real(nslyr,kind=dbl_kind) + hilyr = hi/real(nilyr,kind=dbl_kind) + dzssl_new = min(hslyr/c2,hs_ssl) ! new ssl for snow + dzint_new = hs-dzssl_new + dzssli_new = min(hilyr/c2,hi_ssl) ! new ssl for ice + dzinti_new = hi-dzssli_new + + do k=1,n_iso + + dznew=min(dzssl_new-dzssl,c0) + sloss1=c0 + if (dzssl > puny) & + sloss1=dznew*isosno(k,1)/dzssl ! not neccesarily a loss term + dznew=max(dzssl_new-dzssl,c0) + if (dzint > puny) & + sloss1=sloss1+isosno(k,2)*dznew/dzint ! not really a loss term + isosno(k,1) =isosno(k,1)+sloss1 + isosno(k,2) =isosno(k,2)-sloss1 +! if (isosno(k,2) < c0) then +! write(nu_diag,*) 'redistribute isosno(k,2)',isosno(k,2),sloss1 +! write(nu_diag,*) 'dzssl_new,dzssl,dzint',dzssl_new,dzssl,dzint +! endif + + sloss2=c0 + dznew=min(dzssli_new-dzssli,c0) + if (dzssli > puny) & + sloss2=dznew*isoice(k,1)/dzssli + dznew=max(dzssli_new-dzssli,c0) + if (dzinti > puny) & + sloss2=sloss2+isoice(k,2)*dznew/dzinti + isoice(k,1) = isoice(k,1)+sloss2 + isoice(k,2) = isoice(k,2)-sloss2 + + isotot(k)=isosno(k,2)+isosno(k,1) & + +isoice(k,2)+isoice(k,1) +! if ( (isotot(k)-isotot0(k)) & +! - fiso_atm(k)*aicen*dt & +! - fiso_evapn(k)*dt & +! + fiso_ocnn(k)*dt > 1e-6) then +! write(nu_diag,*) 'isotope tracer: ',k +! write(nu_diag,*) 'isotot-isotot0 ',isotot(k)-isotot0(k) +! write(nu_diag,*) 'fiso_atm-fiso_ocnn ',fiso_atm(k)*aicen*dt & +! + fiso_evapn(k)*dt & +! - fiso_ocnn(k)*dt +! endif + + enddo ! n_iso + +! reload tracers + +! scale fiso_ocnn. It will be re-scaled by aicen latter in merge_fluxes + if (aicen > puny) then + fiso_ocnn(:) = fiso_ocnn(:)/aicen + fiso_evapn(:) = fiso_evapn(:)/aicen + endif + + end subroutine update_isotope + +!======================================================================= + function isoice_alpha(growth_rate, sp, frac) + +! !DESCRIPTION: +! +! calculate the fractionation coefficient for sea-ice formation +! +! !REVISION HISTORY: +! +! authors: Jiang Zhu, UW-Madison +! +!----------------------------------------------------------------------- + + real (kind=dbl_kind), intent(in) :: & + growth_rate ! sea-ice formation rate (m/s) + character(*), intent(in) :: & + sp,frac ! species: H2_16O, H2_18O, HDO + ! calculation methods: + ! cfrac, constant fractionation + ! nfrac, nonfractionation + ! gfrac, growth-rate dependent + real (kind=dbl_kind) :: & + isoice_alpha ! return fractionation + + if (frac == 'nfrac') isoice_alpha = c1 + if (sp == 'H2_16O') isoice_alpha = c1 + + ! Lehmann and Siegenthaler, 1991 + !-------------------------------------------------- + if (frac == 'cfrac' .and. sp == 'HDO') & + isoice_alpha = 1.02120_dbl_kind + if (frac == 'cfrac' .and. sp == 'H2_18O') & + isoice_alpha = 1.00291_dbl_kind + + ! Eq.9, Toyota et al., 2013 + ! For HDO, 7.2852 = 0.2120/0.00291 + !-------------------------------------------------- + if (frac == 'gfrac' .and. sp == 'HDO') & + isoice_alpha = c1+7.2852_dbl_kind*1.2280E-3_dbl_kind+ & + 0.7311E-3_dbl_kind*exp(-growth_rate/8.0100E8_dbl_kind)+ & + 0.8441E-3_dbl_kind*exp(-growth_rate/0.7800E6_dbl_kind) + if (frac == 'gfrac' .and. sp == 'H2_18O') & + isoice_alpha = c1+1.2280E-3_dbl_kind+ & + 0.7311E-3_dbl_kind*exp(-growth_rate/8.0100E8_dbl_kind)+ & + 0.8441E-3_dbl_kind*exp(-growth_rate/0.7800E6_dbl_kind) + return + + end function isoice_alpha + +!======================================================================= + function wiso_alpi(isp,tk) + +!----------------------------------------------------------------------- +! Purpose: return ice/vapour fractionation from loop-up tables +! Author: David Noone - Tue Jul 1 12:02:24 MDT 2003 +!----------------------------------------------------------------------- + integer , intent(in) :: isp ! species indes + real(kind=dbl_kind), intent(in) :: tk ! temperature (k) + real(kind=dbl_kind) :: wiso_alpi ! return fractionation + +!From Merlivat & Nief,1967 for HDO, and Majoube, 1971b for H218O: + real(kind=dbl_kind), parameter, dimension(pwtspec) :: & ! ice/vapour + alpai = (/ 0._dbl_kind, 0._dbl_kind, 16289._dbl_kind, 0._dbl_kind /), & + alpbi = (/ 0._dbl_kind, 0._dbl_kind, 0._dbl_kind, 11.839_dbl_kind /), & + alpci = (/ 0._dbl_kind, 0._dbl_kind, -9.45e-2_dbl_kind, -28.224e-3_dbl_kind /) + +!----------------------------------------------------------------------- + if (isp == isph2o) then + wiso_alpi = c1 + return + end if + + wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp)) + +#ifdef NOFRAC + wiso_alpi = c1 +#endif +! + return + end function wiso_alpi + + +!======================================================================= + + end module icepack_isotope + +!======================================================================= diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 59eb856d1..4514498fe 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -29,7 +29,7 @@ module icepack_itd use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin - use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero + use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_iso use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers use icepack_parameters, only: kcatbound, kitd @@ -760,9 +760,9 @@ subroutine cleanup_itd (dt, ntrcr, & aicen, trcrn, & vicen, vsnon, & aice0, aice, & - n_aero, & + n_aero, n_iso, & nbtrcr, nblyr, & - tr_aero, & + tr_aero, tr_iso, & tr_pond_topo, & heat_capacity, & first_ice, & @@ -770,7 +770,8 @@ subroutine cleanup_itd (dt, ntrcr, & n_trcr_strata,nt_strata, & fpond, fresh, & fsalt, fhocn, & - faero_ocn, fzsal, & + faero_ocn, fiso_ocn, & + fzsal, & flux_bio, limit_aice_in) integer (kind=int_kind), intent(in) :: & @@ -780,7 +781,8 @@ subroutine cleanup_itd (dt, ntrcr, & nslyr , & ! number of snow layers ntrcr , & ! number of tracers in use nbtrcr, & ! number of bio tracers in use - n_aero ! number of aerosol tracers + n_aero, & ! number of aerosol tracers + n_iso ! number of isotope tracers real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -813,6 +815,7 @@ subroutine cleanup_itd (dt, ntrcr, & logical (kind=log_kind), intent(in) :: & tr_aero, & ! aerosol flag + tr_iso, & ! isotope flag tr_pond_topo, & ! topo pond flag heat_capacity ! if false, ice and snow have zero heat capacity @@ -834,6 +837,10 @@ subroutine cleanup_itd (dt, ntrcr, & real (kind=dbl_kind), dimension (:), intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension (:), & + intent(inout), optional :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + logical (kind=log_kind), intent(in), optional :: & limit_aice_in ! if false, allow aice to be out of bounds ! may want to allow this for unit tests @@ -854,6 +861,9 @@ subroutine cleanup_itd (dt, ntrcr, & real (kind=dbl_kind), dimension (n_aero) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) + real (kind=dbl_kind), dimension (n_iso) :: & + dfiso_ocn ! zapped isotope flux (kg/m^2/s) + real (kind=dbl_kind), dimension (ntrcr) :: & dflux_bio ! zapped biology flux (mmol/m^2/s) @@ -877,6 +887,7 @@ subroutine cleanup_itd (dt, ntrcr, & dfsalt = c0 dfhocn = c0 dfaero_ocn(:) = c0 + dfiso_ocn(:) = c0 dflux_bio(:) = c0 dfzsal = c0 @@ -931,7 +942,8 @@ subroutine cleanup_itd (dt, ntrcr, & if (limit_aice) then call zap_small_areas (dt, ntrcr, & - ncat, n_aero, & + ncat, & + n_aero, n_iso, & nblyr, & nilyr, nslyr, & aice, aice0, & @@ -939,8 +951,10 @@ subroutine cleanup_itd (dt, ntrcr, & vicen, vsnon, & dfpond, & dfresh, dfsalt, & - dfhocn, dfaero_ocn, & - tr_aero, tr_pond_topo, & + dfhocn, & + dfaero_ocn, dfiso_ocn, & + tr_aero, tr_iso, & + tr_pond_topo, & first_ice, nbtrcr, & dfzsal, dflux_bio ) @@ -966,8 +980,9 @@ subroutine cleanup_itd (dt, ntrcr, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & + dfiso_ocn, tr_iso, & dflux_bio, nbtrcr, & - n_aero) + n_aero, n_iso) if (icepack_warnings_aborted(subname)) return !------------------------------------------------------------------- @@ -987,6 +1002,13 @@ subroutine cleanup_itd (dt, ntrcr, & faero_ocn(it) = faero_ocn(it) + dfaero_ocn(it) enddo endif + if (present(fiso_ocn)) then + if (tr_iso) then + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it) + enddo + endif + endif if (present(flux_bio)) then do it = 1, nbtrcr flux_bio (it) = flux_bio(it) + dflux_bio(it) @@ -1018,7 +1040,8 @@ end subroutine cleanup_itd ! author: William H. Lipscomb, LANL subroutine zap_small_areas (dt, ntrcr, & - ncat, n_aero, & + ncat, & + n_aero, n_iso, & nblyr, & nilyr, nslyr, & aice, aice0, & @@ -1026,8 +1049,10 @@ subroutine zap_small_areas (dt, ntrcr, & vicen, vsnon, & dfpond, & dfresh, dfsalt, & - dfhocn, dfaero_ocn, & - tr_aero, tr_pond_topo, & + dfhocn, & + dfaero_ocn, dfiso_ocn, & + tr_aero, tr_iso, & + tr_pond_topo, & first_ice, nbtrcr, & dfzsal, dflux_bio ) @@ -1038,6 +1063,7 @@ subroutine zap_small_areas (dt, ntrcr, & nslyr , & ! number of snow layers ntrcr , & ! number of tracers in use n_aero , & ! number of aerosol tracers + n_iso , & ! number of isotope tracers nbtrcr ! number of biology tracers real (kind=dbl_kind), intent(in) :: & @@ -1065,11 +1091,15 @@ subroutine zap_small_areas (dt, ntrcr, & real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) + real (kind=dbl_kind), dimension (:), intent(inout) :: & + dfiso_ocn ! zapped isotope flux (kg/m^2/s) + real (kind=dbl_kind), dimension (:), intent(inout) :: & dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & tr_aero, & ! aerosol flag + tr_iso, & ! isotope flag tr_pond_topo ! pond flag logical (kind=log_kind), dimension (:), intent(inout) :: & @@ -1123,6 +1153,14 @@ subroutine zap_small_areas (dt, ntrcr, & enddo endif + if (tr_iso) then + do it = 1, n_iso + xtmp = (vicen(n)*(trcrn(nt_iso+2+4*(it-1),n) & + + trcrn(nt_iso+3+4*(it-1),n)))/dt + dfiso_ocn(it) = dfiso_ocn(it) + xtmp + enddo + endif + if (solve_zsal) then do it = 1, nblyr xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& @@ -1187,8 +1225,9 @@ subroutine zap_small_areas (dt, ntrcr, & trcrn(:,n), vsnon(n), & dfresh, dfhocn, & dfaero_ocn, tr_aero, & + dfiso_ocn, tr_iso, & dflux_bio, nbtrcr, & - n_aero, & + n_aero, n_iso, & aicen(n), nblyr) if (icepack_warnings_aborted(subname)) return @@ -1245,6 +1284,17 @@ subroutine zap_small_areas (dt, ntrcr, & enddo ! it endif + if (tr_iso) then + do it = 1, n_iso + xtmp = (vsnon(n)*(trcrn(nt_iso +4*(it-1),n) & + + trcrn(nt_iso+1+4*(it-1),n)) & + + vicen(n)*(trcrn(nt_iso+2+4*(it-1),n) & + + trcrn(nt_iso+3+4*(it-1),n))) & + * (aice-c1)/aice / dt + dfiso_ocn(it) = dfiso_ocn(it) + xtmp + enddo ! it + endif + !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- @@ -1319,13 +1369,15 @@ subroutine zap_snow(dt, nslyr, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & + dfiso_ocn, tr_iso, & dflux_bio, nbtrcr, & - n_aero, & + n_aero, n_iso, & aicen, nblyr) integer (kind=int_kind), intent(in) :: & nslyr , & ! number of snow layers n_aero , & ! number of aerosol tracers + n_iso , & ! number of isotope tracers nblyr , & ! number of bio layers nbtrcr @@ -1346,11 +1398,15 @@ subroutine zap_snow(dt, nslyr, & real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) + real (kind=dbl_kind), dimension (:), intent(inout) :: & + dfiso_ocn ! zapped isotope flux (kg/m^2/s) + real (kind=dbl_kind), dimension (:), intent(inout) :: & dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & - tr_aero ! aerosol flag + tr_aero, & ! aerosol flag + tr_iso ! isotope flag ! local variables @@ -1370,6 +1426,15 @@ subroutine zap_snow(dt, nslyr, & enddo ! it endif ! tr_aero + ! isotopes + if (tr_iso) then + do it = 1, n_iso + xtmp = (vsnon*(trcrn(nt_iso +4*(it-1)) & + + trcrn(nt_iso+1+4*(it-1))))/dt + dfiso_ocn(it) = dfiso_ocn(it) + xtmp + enddo ! it + endif ! tr_iso + if (z_tracers) then dvssl = min(p5*vsnon, hs_ssl*aicen) !snow surface layer dvint = vsnon- dvssl !snow interior @@ -1406,13 +1471,15 @@ subroutine zap_snow_temperature(dt, ncat, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & + dfiso_ocn, tr_iso, & dflux_bio, nbtrcr, & - n_aero) + n_aero, n_iso) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nslyr , & ! number of snow layers n_aero, & ! number of aerosol tracers + n_iso, & ! number of isotope tracers nbtrcr, & ! number of z-tracers in use nblyr ! number of bio layers in ice @@ -1439,10 +1506,14 @@ subroutine zap_snow_temperature(dt, ncat, & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout) :: & + dfiso_ocn ! zapped isotope flux (kg/m^2/s) + + real (kind=dbl_kind), dimension (:),intent(inout) :: & dflux_bio ! zapped biology flux (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & - tr_aero ! aerosol flag + tr_aero, & ! aerosol flag + tr_iso ! isotope flag ! local variables @@ -1521,8 +1592,9 @@ subroutine zap_snow_temperature(dt, ncat, & trcrn(:,n), vsnon(n), & dfresh, dfhocn, & dfaero_ocn, tr_aero, & + dfiso_ocn, tr_iso, & dflux_bio, nbtrcr, & - n_aero, & + n_aero, n_iso, & aicen(n), nblyr) if (icepack_warnings_aborted(subname)) return diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 215b3790f..03232a821 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -41,9 +41,9 @@ module icepack_mechred use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg use icepack_parameters, only: heat_capacity - use icepack_tracers, only: tr_pond_topo, tr_aero, tr_brine, ntrcr, nbtrcr + use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice - use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero + use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_iso use icepack_tracers, only: nt_apnd, nt_hpnd use icepack_tracers, only: icepack_compute_tracers @@ -95,6 +95,7 @@ module icepack_mechred subroutine ridge_ice (dt, ndtd, & ncat, n_aero, & + n_iso, & nilyr, nslyr, & ntrcr, hin_max, & rdg_conv, rdg_shear, & @@ -110,7 +111,7 @@ subroutine ridge_ice (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -124,6 +125,7 @@ subroutine ridge_ice (dt, ndtd, & nilyr , & ! number of ice layers nslyr , & ! number of snow layers n_aero, & ! number of aerosol tracers + n_iso, & ! number of isotope tracers ntrcr ! number of tracers in use real (kind=dbl_kind), intent(in) :: & @@ -191,6 +193,9 @@ subroutine ridge_ice (dt, ndtd, & real (kind=dbl_kind), dimension(:), intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + ! local variables real (kind=dbl_kind), dimension (ncat) :: & @@ -221,6 +226,9 @@ subroutine ridge_ice (dt, ndtd, & real (kind=dbl_kind), dimension (n_aero) :: & maero ! aerosol mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension (n_iso) :: & + miso ! isotope mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension (0:ncat) :: & apartic ! participation function; fraction of ridging ! and closing associated w/ category n @@ -270,6 +278,7 @@ subroutine ridge_ice (dt, ndtd, & msnow_mlt = c0 esnow_mlt = c0 maero (:) = c0 + miso (:) = c0 mpond = c0 ardg1 = c0 ardg2 = c0 @@ -394,8 +403,10 @@ subroutine ridge_ice (dt, ndtd, & ardg1n, ardg2n, & virdgn, & nslyr, n_aero, & + n_iso, & msnow_mlt, esnow_mlt, & - maero, mpond, & + maero, miso, & + mpond, & aredistn, vredistn) if (icepack_warnings_aborted(subname)) return @@ -589,6 +600,13 @@ subroutine ridge_ice (dt, ndtd, & faero_ocn(it) = faero_ocn(it) + maero(it)*dti enddo endif + if (present(fiso_ocn)) then + if (tr_iso) then + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti + enddo + endif + endif if (present(fpond)) then fpond = fpond - mpond ! units change later endif @@ -1074,8 +1092,10 @@ subroutine ridge_shift (ntrcr, dt, & ardg1nn, ardg2nn, & virdgnn, & nslyr, n_aero, & + n_iso, & msnow_mlt, esnow_mlt, & - maero, mpond, & + maero, miso, & + mpond, & aredistn, vredistn) integer (kind=int_kind), intent(in) :: & @@ -1083,6 +1103,7 @@ subroutine ridge_shift (ntrcr, dt, & nslyr , & ! number of snow layers ntrcr , & ! number of tracers in use n_aero, & ! number of aerosol tracers + n_iso, & ! number of isotope tracers krdg_redist ! selects redistribution function real (kind=dbl_kind), intent(in) :: & @@ -1147,6 +1168,9 @@ subroutine ridge_shift (ntrcr, dt, & real (kind=dbl_kind), dimension(:), intent(inout) :: & maero ! aerosol mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension(:), intent(inout) :: & + miso ! isotope mass added to ocean (kg m-2) + real (kind=dbl_kind), dimension (:), intent(inout), optional :: & aredistn , & ! redistribution function: fraction of new ridge area vredistn ! redistribution function: fraction of new ridge volume @@ -1362,6 +1386,15 @@ subroutine ridge_shift (ntrcr, dt, & enddo endif + if (tr_iso) then + do it = 1, n_iso + miso(it) = miso(it) & + + vsrdgn*(c1-fsnowrdg) & + *(trcrn(nt_iso +4*(it-1),n) & + + trcrn(nt_iso+1+4*(it-1),n)) + enddo + endif + if (tr_pond_topo) then mpond = mpond + ardg1n * trcrn(nt_apnd,n) & * trcrn(nt_hpnd,n) @@ -1711,8 +1744,8 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - n_aero, & - faero_ocn, & + n_aero, n_iso, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1731,7 +1764,8 @@ subroutine icepack_step_ridge (dt, ndtd, & nblyr , & ! number of bio layers nilyr , & ! number of ice layers nslyr , & ! number of snow layers - n_aero ! number of aerosol tracers + n_aero, & ! number of aerosol tracers + n_iso ! number of isotope tracers real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: & hin_max ! category limits (m) @@ -1781,6 +1815,9 @@ subroutine icepack_step_ridge (dt, ndtd, & faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean + real (kind=dbl_kind), dimension(:), intent(inout) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & trcrn ! tracers @@ -1813,6 +1850,7 @@ subroutine icepack_step_ridge (dt, ndtd, & call ridge_ice (dt, ndtd, & ncat, n_aero, & + n_iso, & nilyr, nslyr, & ntrcr, hin_max, & rdg_conv, rdg_shear, & @@ -1830,7 +1868,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1842,6 +1880,7 @@ subroutine icepack_step_ridge (dt, ndtd, & call ridge_ice (dt, ndtd, & ncat, n_aero, & + n_iso, & nilyr, nslyr, & ntrcr, hin_max, & rdg_conv, rdg_shear, & @@ -1859,7 +1898,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1882,16 +1921,17 @@ subroutine icepack_step_ridge (dt, ndtd, & aicen, trcrn, & vicen, vsnon, & aice0, aice, & - n_aero, & + n_aero, n_iso, & nbtrcr, nblyr, & - tr_aero, & + tr_aero, tr_iso, & tr_pond_topo, heat_capacity, & first_ice, & trcr_depend, trcr_base, & n_trcr_strata, nt_strata, & fpond, fresh, & fsalt, fhocn, & - faero_ocn, fzsal, & + faero_ocn, fiso_ocn, & + fzsal, & flux_bio) if (icepack_warnings_aborted(subname)) return diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 0952e08ca..aa2852c16 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -30,12 +30,12 @@ module icepack_therm_itd use icepack_tracers, only: ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice - use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero + use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_iso use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd use icepack_tracers, only: nt_alvl, nt_vlvl use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo - use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, tr_fsd - use icepack_tracers, only: n_aero + use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd + use icepack_tracers, only: n_aero, n_iso use icepack_tracers, only: bio_index use icepack_warnings, only: warnstr, icepack_warnings_add @@ -45,6 +45,7 @@ module icepack_therm_itd use icepack_itd, only: reduce_area, cleanup_itd use icepack_itd, only: aggregate_area, shift_ice use icepack_itd, only: column_sum, column_conservation_check + use icepack_isotope, only: isoice_alpha, frac use icepack_mushy_physics, only: liquidus_temperature_mush, enthalpy_mush use icepack_therm_shared, only: hfrazilmin use icepack_therm_shared, only: hi_min @@ -853,9 +854,11 @@ end subroutine update_vertical_tracers ! subroutine lateral_melt (dt, ncat, & nilyr, nslyr, & - n_aero, fpond, & + n_aero, n_iso, & + fpond, & fresh, fsalt, & fhocn, faero_ocn, & + fiso_ocn, & rside, meltl, & fside, sss, & aicen, vicen, & @@ -875,6 +878,7 @@ subroutine lateral_melt (dt, ncat, & nblyr , & ! number of bio layers nslyr , & ! number of snow layers n_aero , & ! number of aerosol tracers + n_iso , & ! number of isotope tracers nbtrcr ! number of bio tracers real (kind=dbl_kind), dimension (:), intent(inout) :: & @@ -911,6 +915,9 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), dimension(:), intent(inout) :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), dimension(:), intent(inout) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + ! local variables integer (kind=int_kind) :: & @@ -1167,6 +1174,18 @@ subroutine lateral_melt (dt, ncat, & enddo ! k endif ! tr_aero + if (tr_iso) then + do k = 1, n_iso + fiso_ocn(k) = fiso_ocn(k) + (vsnon(n) & + *(trcrn(nt_iso +4*(k-1),n) & + + trcrn(nt_iso+1+4*(k-1),n)) & + + vicen(n) & + *(trcrn(nt_iso+2+4*(k-1),n) & + + trcrn(nt_iso+3+4*(k-1),n))) & + * rside / dt + enddo ! k + endif ! tr_iso + !----------------------------------------------------------------- ! Biogeochemistry !----------------------------------------------------------------- @@ -1234,7 +1253,7 @@ end subroutine lateral_melt ! subroutine add_new_ice (ncat, nilyr, & nfsd, nblyr, & - n_aero, dt, & + n_aero, n_iso, dt, & ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & @@ -1251,6 +1270,9 @@ subroutine add_new_ice (ncat, nilyr, & nbtrcr, flux_bio, & ocean_bio, fzsal, & frazil_diag, & + fiso_ocn, & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn, & wave_sig_ht, & wave_spectrum, & wavefreq, & @@ -1269,6 +1291,7 @@ subroutine add_new_ice (ncat, nilyr, & ntrcr , & ! number of tracers nltrcr, & ! number of zbgc tracers n_aero, & ! number of aerosol tracers + n_iso, & ! number of isotope tracers ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy) real (kind=dbl_kind), dimension(0:ncat), intent(in) :: & @@ -1336,6 +1359,16 @@ subroutine add_new_ice (ncat, nilyr, & real (kind=dbl_kind), intent(inout) :: & fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) + ! water isotopes + + real (kind=dbl_kind), dimension(:), intent(inout) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + real (kind=dbl_kind), intent(in) :: & + HDO_ocn , & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) + ! floe size distribution real (kind=dbl_kind), intent(in) :: & wave_sig_ht ! significant height of waves globally (m) @@ -1392,6 +1425,8 @@ subroutine add_new_ice (ncat, nilyr, & vice_init, vice_final, & ! ice volume summed over categories eice_init, eice_final ! ice energy summed over categories + real (kind=dbl_kind) :: frazil_conc + real (kind=dbl_kind), dimension (nilyr) :: & Sprofile ! salinity profile used for new ice additions @@ -1666,6 +1701,33 @@ subroutine add_new_ice (ncat, nilyr, & enddo endif + if (tr_iso .and. vtmp > puny) then + do it=1,n_iso + if (it==1) & + frazil_conc = isoice_alpha(c0,'HDO',frac) & + *HDO_ocn + if (it==2) & + frazil_conc = isoice_alpha(c0,'H2_16O',frac) & + *H2_16O_ocn + if (it==3) & + frazil_conc = isoice_alpha(c0,'H2_18O',frac) & + *H2_18O_ocn + + ! dilution in the ssl + trcrn(nt_iso+2+4*(it-1),n) & + = (trcrn(nt_iso+2+4*(it-1),n)*vicen(n)) & + / vtmp + ! dilution and uptake in the int + trcrn(nt_iso+3+4*(it-1),n) & + = (trcrn(nt_iso+3+4*(it-1),n)*vicen(n) & + + frazil_conc*rhoi*vsurp) & + / vtmp + + fiso_ocn(it) = fiso_ocn(it) & + - frazil_conc*rhoi*vsurp/dt + enddo + endif + ! update category volumes vicen(n) = vtmp @@ -1767,6 +1829,29 @@ subroutine add_new_ice (ncat, nilyr, & enddo endif + if (tr_iso) then + do it=1,n_iso + if (it==1) & + frazil_conc = isoice_alpha(c0,'HDO',frac) & + *HDO_ocn + if (it==2) & + frazil_conc = isoice_alpha(c0,'H2_16O',frac) & + *H2_16O_ocn + if (it==3) & + frazil_conc = isoice_alpha(c0,'H2_18O',frac) & + *H2_18O_ocn + + trcrn(nt_iso+2+4*(it-1),1) = & + (trcrn(nt_iso+2+4*(it-1),1)*vice1)/vicen(1) + trcrn(nt_iso+3+4*(it-1),1) = & + (trcrn(nt_iso+3+4*(it-1),1)*vice1) & + +frazil_conc*rhoi*vi0new/vicen(1) + + fiso_ocn(it) = fiso_ocn(it) & + - frazil_conc*rhoi*vi0new/dt + enddo + endif + if (tr_lvl) then alvl = trcrn(nt_alvl,n) trcrn(nt_alvl,n) = & @@ -1880,6 +1965,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & flux_bio, ocean_bio, & frazil_diag, & frz_onset, yday, & + fiso_ocn, HDO_ocn, & + H2_16O_ocn, H2_18O_ocn, & nfsd, wave_sig_ht, & wave_spectrum, & wavefreq, & @@ -1987,6 +2074,14 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & real (kind=dbl_kind), intent(in), optional :: & yday ! day of year + ! water isotopes + real (kind=dbl_kind), dimension(:), intent(inout) :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + + real (kind=dbl_kind), intent(in) :: & + HDO_ocn , & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) !autodocument_end character(len=*),parameter :: subname='(icepack_step_therm2)' @@ -2051,7 +2146,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & call add_new_ice (ncat, nilyr, & nfsd, nblyr, & - n_aero, dt, & + n_aero, n_iso, dt, & ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & @@ -2067,12 +2162,15 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag, & + frazil_diag, fiso_ocn, & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn, & wave_sig_ht, & wave_spectrum, & wavefreq, dwavefreq, & d_afsd_latg, d_afsd_newi, & floe_rad_c, floe_binwidth) + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -2081,9 +2179,11 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & call lateral_melt (dt, ncat, & nilyr, nslyr, & - n_aero, fpond, & + n_aero, n_iso, & + fpond, & fresh, fsalt, & fhocn, faero_ocn, & + fiso_ocn, & rside, meltl, & fside, sss, & aicen, vicen, & @@ -2126,17 +2226,17 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & aicen, trcrn(1:ntrcr,:), & vicen, vsnon, & aice0, aice, & - n_aero, & + n_aero, n_iso, & nbtrcr, nblyr, & - tr_aero, & + tr_aero, tr_iso, & tr_pond_topo, heat_capacity, & first_ice, & trcr_depend, trcr_base, & n_trcr_strata, nt_strata, & fpond, fresh, & fsalt, fhocn, & - faero_ocn, fzsal, & - flux_bio) + faero_ocn, fiso_ocn, & + fzsal, flux_bio) if (icepack_warnings_aborted(subname)) return end subroutine icepack_step_therm2 diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 907ef565a..2cb99f915 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -28,9 +28,9 @@ module icepack_therm_vertical use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd use icepack_parameters, only: phi_i_mushy - use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd + use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_iso, tr_pond, tr_fsd use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo - use icepack_tracers, only: n_aero + use icepack_tracers, only: n_aero, n_iso use icepack_therm_shared, only: ferrmax, l_brine use icepack_therm_shared, only: calculate_tin_from_qin, Tmin @@ -47,6 +47,7 @@ module icepack_therm_vertical use icepack_mushy_physics, only: enthalpy_mush, enthalpy_of_melting use icepack_aerosol, only: update_aerosol + use icepack_isotope, only: update_isotope use icepack_atmo, only: neutral_drag_coeffs, icepack_atm_boundary use icepack_age, only: increment_age use icepack_firstyear, only: update_FYarea @@ -2040,11 +2041,14 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & ipnd , & iage , FY , & aerosno , aeroice , & + isosno , isoice , & uatm , vatm , & wind , zlvl , & Qa , rhoa , & + Qa_iso , & Tair , Tref , & Qref , Uref , & + Qref_iso , & Cdn_atm_ratio, & Cdn_ocn , Cdn_ocn_skin, & Cdn_ocn_floe, Cdn_ocn_keel, & @@ -2083,6 +2087,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & flatn_f , fsensn_f , & fsurfn_f , fcondtopn_f , & faero_atm , faero_ocn , & + fiso_atm , fiso_ocn , & + fiso_evap , & + HDO_ocn , H2_16O_ocn , & + H2_18O_ocn , & dhsn , ffracn , & meltt , melttn , & meltb , meltbn , & @@ -2186,6 +2194,18 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & mlt_onset , & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) + real (kind=dbl_kind), dimension(n_iso), intent(inout) :: & + Qa_iso , & ! isotope specific humidity (kg/kg) + Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg) + fiso_atm , & ! isotope deposition rate (kg/m^2 s) + fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) + fiso_evap ! isotope evaporation (kg/m^2/s) + + real (kind=dbl_kind), intent(in) :: & + HDO_ocn , & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) + real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice vicen_init , & ! volume per unit area of ice (m) @@ -2235,6 +2255,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & aerosno , & ! snow aerosol tracer (kg/m^2) aeroice ! ice aerosol tracer (kg/m^2) + real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & + isosno , & ! snow isotope tracer (kg/m^2) + isoice ! ice isotope tracer (kg/m^2) !autodocument_end ! local variables @@ -2265,6 +2288,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & lhcoef , & ! transfer coefficient for latent heat rfrac ! water fraction retained for melt ponds + real (kind=dbl_kind), dimension(n_iso) :: & + Qrefn_iso , & ! isotope air sp hum reference level (kg/kg) + fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s) + fiso_evapn ! isotope evaporation (kg/m^2/s) + real (kind=dbl_kind) :: & pond ! water retained in ponds (m) @@ -2326,6 +2354,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & Trefn = c0 Qrefn = c0 + Qrefn_iso(:) = c0 + fiso_ocnn(:) = c0 + fiso_evapn(:) = c0 Urefn = c0 lhcoef = c0 shcoef = c0 @@ -2356,7 +2387,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - uvel, vvel, & + n_iso=n_iso, tr_iso=tr_iso, & + Qa_iso=Qa_iso(:), Qref_iso=Qrefn_iso(:), & + uvel=uvel, vvel=vvel, & Uref=Urefn) if (icepack_warnings_aborted(subname)) return @@ -2469,6 +2502,26 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (icepack_warnings_aborted(subname)) return endif + if (tr_iso) then + call update_isotope (dt = dt, & + nilyr = nilyr, nslyr = nslyr, n_iso = n_iso, & + meltt = melttn(n),melts = meltsn(n), & + meltb = meltbn(n),congel=congeln(n), & + snoice=snoicen(n),evap=evapn, & + fsnow=fsnow, Tsfc=Tsfc(n), & + Qref_iso=Qrefn_iso(:), & + isosno=isosno(:,:,n),isoice=isoice(:,:,n), & + aice_old=aicen_init(n),vice_old=vicen_init(n), & + vsno_old=vsnon_init(n), & + vicen=vicen(n),vsnon=vsnon(n), & + aicen=aicen(n), & + fiso_atm=fiso_atm(:), & + fiso_evapn=fiso_evapn(:), & + fiso_ocnn=fiso_ocnn(:), & + HDO_ocn=HDO_ocn,H2_16O_ocn=H2_16O_ocn, & + H2_18O_ocn=H2_18O_ocn) + + endif endif ! aicen_init !----------------------------------------------------------------- @@ -2570,7 +2623,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & meltt, melts, & meltb, congel, & snoice, & - Uref, Urefn) + Uref=Uref, Urefn=Urefn, & + Qref_iso=Qref_iso(:),Qrefn_iso=Qrefn_iso(:), & + fiso_ocn=fiso_ocn(:),fiso_ocnn=fiso_ocnn(:), & + fiso_evap=fiso_evap(:),fiso_evapn=fiso_evapn(:)) + if (icepack_warnings_aborted(subname)) return enddo ! ncat diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 04adc96d1..f100b9462 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -37,6 +37,7 @@ module icepack_tracers nmodal1 = 10 , & ! dimension for modal aerosol radiation parameters nmodal2 = 8 , & ! dimension for modal aerosol radiation parameters max_aero = 6 , & ! maximum number of aerosols + max_iso = 3 , & ! maximum number of isotopes max_nbtrcr = max_algae*2 & ! algal nitrogen and chlorophyll + max_dic & ! dissolved inorganic carbon + max_doc & ! dissolved organic carbon @@ -55,6 +56,7 @@ module icepack_tracers nblyr = 0, & ! number of bio/brine layers per category nfsd = 0, & ! number of fsd layers n_aero = 0, & ! number of aerosols in use + n_iso = 0, & ! number of isotopes in use n_zaero = 0, & ! number of z aerosols in use n_algae = 0, & ! number of algae in use n_doc = 0, & ! number of DOC pools in use @@ -78,6 +80,7 @@ module icepack_tracers nt_ipnd = 0, & ! melt pond refrozen lid thickness nt_fsd = 0, & ! floe size distribution nt_aero = 0, & ! starting index for aerosols in ice + nt_iso = 0, & ! starting index for isotopes in ice nt_bgc_Nit = 0, & ! nutrients nt_bgc_Am = 0, & ! nt_bgc_Sil = 0, & ! @@ -98,6 +101,7 @@ module icepack_tracers tr_pond_lvl = .false., & ! if .true., use level-ice pond tracer tr_pond_topo = .false., & ! if .true., use explicit topography-based ponds tr_aero = .false., & ! if .true., use aerosol tracers + tr_iso = .false., & ! if .true., use isotope tracers tr_brine = .false., & ! if .true., brine height differs from ice thickness tr_fsd = .false. ! if .true., use floe size distribution @@ -197,7 +201,7 @@ module icepack_tracers subroutine icepack_init_tracer_flags(& tr_iage_in, tr_FY_in, tr_lvl_in, & tr_pond_in, tr_pond_cesm_in, tr_pond_lvl_in, tr_pond_topo_in, & - tr_fsd_in, tr_aero_in, tr_brine_in, tr_zaero_in, & + tr_fsd_in, tr_aero_in, tr_iso_in, tr_brine_in, tr_zaero_in, & tr_bgc_Nit_in, tr_bgc_N_in, tr_bgc_DON_in, tr_bgc_C_in, tr_bgc_chl_in, & tr_bgc_Am_in, tr_bgc_Sil_in, tr_bgc_DMS_in, tr_bgc_Fe_in, tr_bgc_hum_in, & tr_bgc_PON_in) @@ -212,6 +216,7 @@ subroutine icepack_init_tracer_flags(& tr_pond_topo_in , & ! if .true., use explicit topography-based ponds tr_fsd_in , & ! if .true., use floe size distribution tracers tr_aero_in , & ! if .true., use aerosol tracers + tr_iso_in , & ! if .true., use isotope tracers tr_brine_in , & ! if .true., brine height differs from ice thickness tr_zaero_in , & ! if .true., black carbon is tracers (n_zaero) tr_bgc_Nit_in , & ! if .true., Nitrate tracer in ice @@ -239,6 +244,7 @@ subroutine icepack_init_tracer_flags(& if (present(tr_pond_topo_in)) tr_pond_topo = tr_pond_topo_in if (present(tr_fsd_in) ) tr_fsd = tr_fsd_in if (present(tr_aero_in) ) tr_aero = tr_aero_in + if (present(tr_iso_in) ) tr_iso = tr_iso_in if (present(tr_brine_in) ) tr_brine = tr_brine_in if (present(tr_zaero_in) ) tr_zaero = tr_zaero_in if (present(tr_bgc_Nit_in)) tr_bgc_Nit = tr_bgc_Nit_in @@ -262,7 +268,7 @@ end subroutine icepack_init_tracer_flags subroutine icepack_query_tracer_flags(& tr_iage_out, tr_FY_out, tr_lvl_out, & tr_pond_out, tr_pond_cesm_out, tr_pond_lvl_out, tr_pond_topo_out, & - tr_fsd_out, tr_aero_out, tr_brine_out, tr_zaero_out, & + tr_fsd_out, tr_aero_out, tr_iso_out, tr_brine_out, tr_zaero_out, & tr_bgc_Nit_out, tr_bgc_N_out, tr_bgc_DON_out, tr_bgc_C_out, tr_bgc_chl_out, & tr_bgc_Am_out, tr_bgc_Sil_out, tr_bgc_DMS_out, tr_bgc_Fe_out, tr_bgc_hum_out, & tr_bgc_PON_out) @@ -277,6 +283,7 @@ subroutine icepack_query_tracer_flags(& tr_pond_topo_out , & ! if .true., use explicit topography-based ponds tr_fsd_out , & ! if .true., use floe size distribution tr_aero_out , & ! if .true., use aerosol tracers + tr_iso_out , & ! if .true., use isotope tracers tr_brine_out , & ! if .true., brine height differs from ice thickness tr_zaero_out , & ! if .true., black carbon is tracers (n_zaero) tr_bgc_Nit_out , & ! if .true., Nitrate tracer in ice @@ -304,6 +311,7 @@ subroutine icepack_query_tracer_flags(& if (present(tr_pond_topo_out)) tr_pond_topo_out = tr_pond_topo if (present(tr_fsd_out) ) tr_fsd_out = tr_fsd if (present(tr_aero_out) ) tr_aero_out = tr_aero + if (present(tr_iso_out) ) tr_iso_out = tr_iso if (present(tr_brine_out) ) tr_brine_out = tr_brine if (present(tr_zaero_out) ) tr_zaero_out = tr_zaero if (present(tr_bgc_Nit_out)) tr_bgc_Nit_out = tr_bgc_Nit @@ -342,6 +350,7 @@ subroutine icepack_write_tracer_flags(iounit) write(iounit,*) " tr_pond_topo = ",tr_pond_topo write(iounit,*) " tr_fsd = ",tr_fsd write(iounit,*) " tr_aero = ",tr_aero + write(iounit,*) " tr_iso = ",tr_iso write(iounit,*) " tr_brine = ",tr_brine write(iounit,*) " tr_zaero = ",tr_zaero write(iounit,*) " tr_bgc_Nit = ",tr_bgc_Nit @@ -366,7 +375,7 @@ subroutine icepack_init_tracer_indices(& nt_Tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in, & nt_fbri_in, nt_iage_in, nt_FY_in, & nt_alvl_in, nt_vlvl_in, nt_apnd_in, nt_hpnd_in, nt_ipnd_in, & - nt_fsd_in, nt_aero_in, nt_zaero_in, nt_bgc_C_in, & + nt_fsd_in, nt_aero_in, nt_iso_in, nt_zaero_in, nt_bgc_C_in, & nt_bgc_N_in, nt_bgc_chl_in, nt_bgc_DOC_in, nt_bgc_DON_in, & nt_bgc_DIC_in, nt_bgc_Fed_in, nt_bgc_Fep_in, nt_bgc_Nit_in, nt_bgc_Am_in, & nt_bgc_Sil_in, nt_bgc_DMSPp_in, nt_bgc_DMSPd_in, nt_bgc_DMS_in, nt_bgc_hum_in, & @@ -393,6 +402,7 @@ subroutine icepack_init_tracer_indices(& nt_ipnd_in, & ! melt pond refrozen lid thickness nt_fsd_in, & ! floe size distribution nt_aero_in, & ! starting index for aerosols in ice + nt_iso_in, & ! starting index for isotopes in ice nt_bgc_Nit_in, & ! nutrients nt_bgc_Am_in, & ! nt_bgc_Sil_in, & ! @@ -468,6 +478,7 @@ subroutine icepack_init_tracer_indices(& if (present(nt_ipnd_in)) nt_ipnd = nt_ipnd_in if (present(nt_fsd_in) ) nt_fsd = nt_fsd_in if (present(nt_aero_in)) nt_aero = nt_aero_in + if (present(nt_iso_in)) nt_iso = nt_iso_in if (present(nt_bgc_Nit_in) ) nt_bgc_Nit = nt_bgc_Nit_in if (present(nt_bgc_Am_in) ) nt_bgc_Am = nt_bgc_Am_in if (present(nt_bgc_Sil_in) ) nt_bgc_Sil = nt_bgc_Sil_in @@ -714,7 +725,7 @@ subroutine icepack_query_tracer_indices(& nt_Tsfc_out, nt_qice_out, nt_qsno_out, nt_sice_out, & nt_fbri_out, nt_iage_out, nt_FY_out, & nt_alvl_out, nt_vlvl_out, nt_apnd_out, nt_hpnd_out, nt_ipnd_out, & - nt_fsd_out, nt_aero_out, nt_zaero_out, nt_bgc_C_out, & + nt_fsd_out, nt_aero_out, nt_iso_out, nt_zaero_out, nt_bgc_C_out, & nt_bgc_N_out, nt_bgc_chl_out, nt_bgc_DOC_out, nt_bgc_DON_out, & nt_bgc_DIC_out, nt_bgc_Fed_out, nt_bgc_Fep_out, nt_bgc_Nit_out, nt_bgc_Am_out, & nt_bgc_Sil_out, nt_bgc_DMSPp_out, nt_bgc_DMSPd_out, nt_bgc_DMS_out, nt_bgc_hum_out, & @@ -741,6 +752,7 @@ subroutine icepack_query_tracer_indices(& nt_ipnd_out, & ! melt pond refrozen lid thickness nt_fsd_out, & ! floe size distribution nt_aero_out, & ! starting index for aerosols in ice + nt_iso_out, & ! starting index for isotopes in ice nt_bgc_Nit_out, & ! nutrients nt_bgc_Am_out, & ! nt_bgc_Sil_out, & ! @@ -814,6 +826,7 @@ subroutine icepack_query_tracer_indices(& if (present(nt_ipnd_out)) nt_ipnd_out = nt_ipnd if (present(nt_fsd_out) ) nt_fsd_out = nt_fsd if (present(nt_aero_out)) nt_aero_out = nt_aero + if (present(nt_iso_out)) nt_iso_out = nt_iso if (present(nt_bgc_Nit_out) ) nt_bgc_Nit_out = nt_bgc_Nit if (present(nt_bgc_Am_out) ) nt_bgc_Am_out = nt_bgc_Am if (present(nt_bgc_Sil_out) ) nt_bgc_Sil_out = nt_bgc_Sil @@ -887,6 +900,7 @@ subroutine icepack_write_tracer_indices(iounit) write(iounit,*) " nt_ipnd = ",nt_ipnd write(iounit,*) " nt_fsd = ",nt_fsd write(iounit,*) " nt_aero = ",nt_aero + write(iounit,*) " nt_iso = ",nt_iso write(iounit,*) " nt_bgc_Nit = ",nt_bgc_Nit write(iounit,*) " nt_bgc_Am = ",nt_bgc_Am write(iounit,*) " nt_bgc_Sil = ",nt_bgc_Sil @@ -964,7 +978,7 @@ end subroutine icepack_write_tracer_indices subroutine icepack_init_tracer_sizes(& ncat_in, nilyr_in, nslyr_in, nblyr_in, nfsd_in , & - n_algae_in, n_DOC_in, n_aero_in, & + n_algae_in, n_DOC_in, n_aero_in, n_iso_in, & n_DON_in, n_DIC_in, n_fed_in, n_fep_in, n_zaero_in, & ntrcr_in, ntrcr_o_in, nbtrcr_in, nbtrcr_sw_in) @@ -982,6 +996,7 @@ subroutine icepack_init_tracer_sizes(& n_fep_in , & ! n_zaero_in, & ! n_aero_in , & ! + n_iso_in , & ! ntrcr_in , & ! number of tracers in use ntrcr_o_in, & ! number of non-bio tracers in use nbtrcr_in , & ! number of bio tracers in use @@ -1005,6 +1020,7 @@ subroutine icepack_init_tracer_sizes(& if (present(n_fep_in) ) n_fep = n_fep_in if (present(n_zaero_in) ) n_zaero = n_zaero_in if (present(n_aero_in) ) n_aero = n_aero_in + if (present(n_iso_in) ) n_iso = n_iso_in if (present(ntrcr_in) ) ntrcr = ntrcr_in if (present(ntrcr_o_in) ) ntrcr_o = ntrcr_o_in @@ -1022,7 +1038,7 @@ subroutine icepack_query_tracer_sizes(& max_don_out , max_fe_out , nmodal1_out , & nmodal2_out , max_aero_out , max_nbtrcr_out , & ncat_out, nilyr_out, nslyr_out, nblyr_out, nfsd_out, & - n_algae_out, n_DOC_out, n_aero_out, & + n_algae_out, n_DOC_out, n_aero_out, n_iso_out, & n_DON_out, n_DIC_out, n_fed_out, n_fep_out, n_zaero_out, & ntrcr_out, ntrcr_o_out, nbtrcr_out, nbtrcr_sw_out) @@ -1051,6 +1067,7 @@ subroutine icepack_query_tracer_sizes(& n_fep_out , & ! n_zaero_out, & ! n_aero_out , & ! + n_iso_out , & ! ntrcr_out , & ! number of tracers in use ntrcr_o_out, & ! number of non-bio tracers in use nbtrcr_out , & ! number of bio tracers in use @@ -1084,6 +1101,7 @@ subroutine icepack_query_tracer_sizes(& if (present(n_fep_out) ) n_fep_out = n_fep if (present(n_zaero_out) ) n_zaero_out = n_zaero if (present(n_aero_out) ) n_aero_out = n_aero + if (present(n_iso_out) ) n_iso_out = n_iso if (present(ntrcr_out) ) ntrcr_out = ntrcr if (present(ntrcr_o_out) ) ntrcr_o_out = ntrcr_o @@ -1114,6 +1132,7 @@ subroutine icepack_write_tracer_sizes(iounit) write(iounit,*) " nmodal1_out =", nmodal1 write(iounit,*) " nmodal2_out =", nmodal2 write(iounit,*) " max_aero_out =", max_aero + write(iounit,*) " max_iso_out =", max_iso write(iounit,*) " max_nbtrcr_out=", max_nbtrcr write(iounit,*) " model defined parameters: " @@ -1130,6 +1149,7 @@ subroutine icepack_write_tracer_sizes(iounit) write(iounit,*) " n_fep = ",n_fep write(iounit,*) " n_zaero = ",n_zaero write(iounit,*) " n_aero = ",n_aero + write(iounit,*) " n_iso = ",n_iso write(iounit,*) " ntrcr = ",ntrcr write(iounit,*) " ntrcr_o = ",ntrcr_o write(iounit,*) " nbtrcr = ",nbtrcr diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index 70d9c50b6..ab4e177f1 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -43,7 +43,7 @@ subroutine icedrv_initialize use icedrv_flux, only: init_coupler_flux, init_history_therm, & init_flux_atm_ocn use icedrv_forcing, only: init_forcing, get_forcing, get_wave_spec - use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, init_forcing_bgc + use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, fiso_default, init_forcing_bgc use icedrv_restart_shared, only: restart use icedrv_init, only: input_data, init_state, init_grid2, init_fsd use icedrv_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc @@ -53,6 +53,7 @@ subroutine icedrv_initialize skl_bgc, & ! from icepack z_tracers, & ! from icepack tr_aero, & ! from icepack + tr_iso, & ! from icepack tr_zaero, & ! from icepack tr_fsd, wave_spec @@ -127,6 +128,7 @@ subroutine icedrv_initialize call icepack_query_parameters(z_tracers_out=z_tracers) call icepack_query_parameters(wave_spec_out=wave_spec) call icepack_query_tracer_flags(tr_aero_out=tr_aero) + call icepack_query_tracer_flags(tr_iso_out=tr_iso) call icepack_query_tracer_flags(tr_zaero_out=tr_zaero) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -141,6 +143,7 @@ subroutine icedrv_initialize ! if (tr_aero) call faero_data ! data file ! if (tr_zaero) call fzaero_data ! data file (gx1) if (tr_aero .or. tr_zaero) call faero_default ! default values + if (tr_iso) call fiso_default ! default values if (skl_bgc .or. z_tracers) call get_forcing_bgc ! biogeochemistry if (.not. restart) & diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index 3c03f7faa..83021df95 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -33,10 +33,10 @@ subroutine icedrv_run use icedrv_calendar, only: istep, istep1, time, dt, stop_now, calendar use icedrv_forcing, only: get_forcing, get_wave_spec - use icedrv_forcing_bgc, only: faero_default, get_forcing_bgc + use icedrv_forcing_bgc, only: faero_default, fiso_default, get_forcing_bgc use icedrv_flux, only: init_flux_atm_ocn - logical (kind=log_kind) :: skl_bgc, z_tracers, tr_aero, tr_zaero, & + logical (kind=log_kind) :: skl_bgc, z_tracers, tr_aero, tr_iso, tr_zaero, & wave_spec, tr_fsd character(len=*), parameter :: subname='(icedrv_run)' @@ -46,7 +46,7 @@ subroutine icedrv_run !-------------------------------------------------------------------- call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero, & - tr_fsd_out=tr_fsd) + tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -75,6 +75,8 @@ subroutine icedrv_run ! aerosols if (tr_aero .or. tr_zaero) call faero_default ! default values + if (tr_iso) call fiso_default ! default values + if (skl_bgc .or. z_tracers) call get_forcing_bgc ! biogeochemistry call init_flux_atm_ocn ! initialize atmosphere, ocean fluxes diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90 index 123ac8c56..e697dee35 100644 --- a/configuration/driver/icedrv_diagnostics.F90 +++ b/configuration/driver/icedrv_diagnostics.F90 @@ -9,6 +9,7 @@ module icedrv_diagnostics use icedrv_kinds use icedrv_constants, only: nu_diag, nu_diag_out use icedrv_domain_size, only: nx + use icedrv_domain_size, only: ncat, nfsd, n_iso use icepack_intfc, only: c0 use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_parameters @@ -41,6 +42,9 @@ module icedrv_diagnostics pdhs , & ! change in mean snow thickness (m) pde ! change in ice and snow energy (W m-2) + real (kind=dbl_kind), dimension(nx,n_iso) :: & + pdiso ! change in mean isotope concentration + !======================================================================= contains @@ -56,10 +60,9 @@ module icedrv_diagnostics subroutine runtime_diags (dt) use icedrv_arrays_column, only: floe_rad_c - use icedrv_domain_size, only: ncat, nfsd use icedrv_flux, only: evap, fsnow, frazil use icedrv_flux, only: fswabs, flw, flwout, fsens, fsurf, flat - use icedrv_flux, only: frain + use icedrv_flux, only: frain, fiso_evap, fiso_ocn, fiso_atm use icedrv_flux, only: Tair, Qa, fsw, fcondtop use icedrv_flux, only: meltt, meltb, meltl, snoice use icedrv_flux, only: dsnow, congel, sst, sss, Tf, fhocn @@ -74,7 +77,7 @@ subroutine runtime_diags (dt) n, nc, k logical (kind=log_kind) :: & - calc_Tsfc, tr_fsd + calc_Tsfc, tr_fsd, tr_iso ! fields at diagnostic points real (kind=dbl_kind) :: & @@ -83,13 +86,13 @@ subroutine runtime_diags (dt) pevap, pfhocn, fsdavg real (kind=dbl_kind), dimension (nx) :: & - work1, work2 + work1, work2, work3 real (kind=dbl_kind) :: & Tffresh, rhos, rhow, rhoi logical (kind=log_kind) :: tr_brine - integer (kind=int_kind) :: nt_fbri, nt_Tsfc, nt_fsd + integer (kind=int_kind) :: nt_fbri, nt_Tsfc, nt_fsd, nt_iso character(len=*), parameter :: subname='(runtime_diags)' @@ -98,9 +101,9 @@ subroutine runtime_diags (dt) !----------------------------------------------------------------- call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc) - call icepack_query_tracer_flags(tr_brine_out=tr_brine,tr_fsd_out=tr_fsd) + call icepack_query_tracer_flags(tr_brine_out=tr_brine,tr_fsd_out=tr_fsd,tr_iso_out=tr_iso) call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc,& - nt_fsd_out=nt_fsd) + nt_fsd_out=nt_fsd,nt_iso_out=nt_iso) call icepack_query_parameters(Tffresh_out=Tffresh, rhos_out=rhos, & rhow_out=rhow, rhoi_out=rhoi) call icepack_warnings_flush(nu_diag) @@ -148,6 +151,16 @@ subroutine runtime_diags (dt) pde(n) =-(work1(n)- pde(n))/dt ! ice/snow energy change pfhocn = -fhocn(n) ! ocean heat used by ice + work3(:) = c0 + + do k = 1, n_iso + work3 (n) = (trcr(n,nt_iso +4*(k-1))*vsno(n) & + +trcr(n,nt_iso+1+4*(k-1))*vsno(n) & + +trcr(n,nt_iso+2+4*(k-1))*vice(n) & + +trcr(n,nt_iso+3+4*(k-1))*vice(n)) + pdiso(n,k) = work3(n) - pdiso(n,k) + enddo + !----------------------------------------------------------------- ! start spewing !----------------------------------------------------------------- @@ -202,6 +215,15 @@ subroutine runtime_diags (dt) write(nu_diag_out+n-1,900) 'freezing temp (C) = ',Tf(n) ! freezing temperature write(nu_diag_out+n-1,900) 'heat used (W/m^2) = ',pfhocn ! ocean heat used by ice + if (tr_iso) then + do k = 1,n_iso + write(nu_diag_out+n-1,900) 'Isotopic precip = ',fiso_atm(n,k)*dt + write(nu_diag_out+n-1,900) 'Isotopic evap/cond = ',fiso_evap(n,k)*dt + write(nu_diag_out+n-1,900) 'Isotopic loss to ocn = ',fiso_ocn(n,k)*dt + write(nu_diag_out+n-1,900) 'Isotopic gain/loss = ',(fiso_atm(n,k)-fiso_ocn(n,k)+fiso_evap(n,k))*dt + write(nu_diag_out+n-1,900) 'Isotopic conc chg = ',pdiso(n,k) + enddo + endif end do 899 format (43x,a24) 900 format (a25,2x,f24.17) @@ -216,20 +238,27 @@ end subroutine runtime_diags subroutine init_mass_diags - use icedrv_domain_size, only: nx - use icedrv_state, only: vice, vsno + use icedrv_state, only: vice, vsno, trcr - integer (kind=int_kind) :: i + integer (kind=int_kind) :: i, k, nt_iso real (kind=dbl_kind), dimension (nx) :: work1 character(len=*), parameter :: subname='(init_mass_diags)' + call icepack_query_tracer_indices(nt_iso_out=nt_iso) + call total_energy (work1) do i = 1, nx pdhi(i) = vice (i) pdhs(i) = vsno (i) pde (i) = work1(i) + do k = 1, n_iso + pdiso(i,k) = (trcr(i,nt_iso +4*(k-1))*vsno(i) & + +trcr(i,nt_iso+1+4*(k-1))*vsno(i) & + +trcr(i,nt_iso+2+4*(k-1))*vice(i) & + +trcr(i,nt_iso+3+4*(k-1))*vice(i)) + enddo enddo end subroutine init_mass_diags @@ -242,7 +271,7 @@ end subroutine init_mass_diags subroutine total_energy (work) - use icedrv_domain_size, only: ncat, nilyr, nslyr, nx + use icedrv_domain_size, only: ncat, nilyr, nslyr use icedrv_state, only: vicen, vsnon, trcrn real (kind=dbl_kind), dimension (nx), intent(out) :: & @@ -303,7 +332,7 @@ end subroutine total_energy subroutine total_salt (work) - use icedrv_domain_size, only: ncat, nilyr, nx + use icedrv_domain_size, only: ncat, nilyr use icedrv_state, only: vicen, trcrn real (kind=dbl_kind), dimension (nx), & @@ -390,7 +419,7 @@ end subroutine icedrv_diagnostics_debug subroutine print_state(plabel,i) use icedrv_calendar, only: istep1, time - use icedrv_domain_size, only: ncat, nilyr, nslyr, nfsd + use icedrv_domain_size, only: nilyr, nslyr use icedrv_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, trcrn use icedrv_flux, only: uatm, vatm, potT, Tair, Qa, flw, frain, fsnow use icedrv_flux, only: fsens, flat, evap, flwout diff --git a/configuration/driver/icedrv_domain_size.F90 b/configuration/driver/icedrv_domain_size.F90 index 37b825e27..fc5c83d79 100644 --- a/configuration/driver/icedrv_domain_size.F90 +++ b/configuration/driver/icedrv_domain_size.F90 @@ -20,6 +20,7 @@ module icedrv_domain_size nilyr = NICELYR , & ! number of ice layers per category nslyr = NSNWLYR , & ! number of snow layers per category n_aero = NTRAERO , & ! number of aerosols in use + n_iso = NTRISO , & ! number of isotopes in use n_zaero = TRZAERO , & ! number of z aerosols in use n_algae = TRALG , & ! number of algae in use n_doc = TRDOC , & ! number of DOC pools in use @@ -48,6 +49,7 @@ module icedrv_domain_size + TRLVL*2 & ! level/deformed ice + TRPND*3 & ! ponds + n_aero*4 & ! number of aerosols * 4 aero layers + + n_iso*4 & ! number of isotopes * 4 isotope layers + TRBRI & ! brine height + TRBGCS*n_bgc & ! skeletal layer BGC + TRZS *TRBRI* nblyr & ! zsalinity (off if TRBRI=0) diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 0abecd448..d1bff5316 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -12,7 +12,7 @@ module icedrv_flux use icedrv_domain_size, only: ncat, nilyr, nx use icedrv_constants, only: c0, c1, c5, c10, c20, c180 use icedrv_constants, only: nu_diag - use icepack_intfc, only: icepack_max_aero, icepack_max_nbtrcr, icepack_max_fe + use icepack_intfc, only: icepack_max_aero, icepack_max_iso, icepack_max_nbtrcr, icepack_max_fe use icepack_intfc, only: icepack_max_algae, icepack_max_doc, icepack_max_don use icepack_intfc, only: icepack_max_dic use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -129,8 +129,11 @@ module icedrv_flux qdp , & ! deep ocean heat flux (W/m^2), negative upward hmix ! mixed layer depth (m) - ! out to atmosphere (if calc_Tsfc) - ! note Tsfc is in ice_state.F + ! water isotopes + real (kind=dbl_kind), dimension (nx), public :: & + HDO_ocn , & ! seawater concentration of HDO (kg/kg) + H2_16O_ocn , & ! seawater concentration of H2_16O (kg/kg) + H2_18O_ocn ! seawater concentration of H2_18O (kg/kg) real (kind=dbl_kind), dimension (nx), public :: & fsens , & ! sensible heat flux (W/m^2) @@ -176,6 +179,10 @@ module icedrv_flux fhocn , & ! net heat flux to ocean (W/m^2) fswthru ! shortwave penetrating to ocean (W/m^2) + real (kind=dbl_kind), dimension (nx,icepack_max_iso), public :: & + Qa_iso , & ! isotope specific humidity (kg/kg) + Qref_iso ! 2m atm reference isotope spec humidity (kg/kg) + ! internal real (kind=dbl_kind), & @@ -286,6 +293,11 @@ module icedrv_flux dimension (nx,icepack_max_aero), public :: & faero_atm ! aerosol deposition rate (kg/m^2 s) + real (kind=dbl_kind), & !coupling variable for tr_iso + dimension (nx,icepack_max_iso), public :: & + fiso_atm , & ! isotope deposition rate (kg/m^2 s) + fiso_evap ! isotope evaporation rate (kg/m^2 s) + real (kind=dbl_kind), & dimension (nx,icepack_max_nbtrcr), public :: & flux_bio_atm ! all bio fluxes to ice from atmosphere @@ -296,6 +308,10 @@ module icedrv_flux dimension (nx,icepack_max_aero), public :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) + real (kind=dbl_kind), & + dimension (nx,icepack_max_iso), public :: & + fiso_ocn ! isotope flux to ocean (kg/m^2/s) + ! out to ocean real (kind=dbl_kind), & @@ -460,6 +476,8 @@ subroutine init_coupler_flux endif ! l_winter faero_atm (:,:) = c0 ! aerosol deposition rate (kg/m2/s) + fiso_atm (:,:) = c0 ! isotope deposition rate (kg/m2/s) + fiso_evap (:,:) = c0 ! isotope evaporation rate (kg/m2/s) flux_bio_atm (:,:) = c0 ! zaero and bio deposition rate (kg/m2/s) !----------------------------------------------------------------- @@ -473,6 +491,11 @@ subroutine init_coupler_flux sst (:) = -1.8_dbl_kind ! sea surface temperature (C) sstdat (:) = sst(:) ! sea surface temperature (C) + ! water isotopes from ocean + HDO_ocn(:) = c0 + H2_16O_ocn(:) = c0 + H2_18O_ocn(:) = c0 + do i = 1, nx Tf (i) = icepack_liquidus_temperature(sss(i)) ! freezing temp (C) enddo @@ -565,6 +588,7 @@ subroutine init_flux_atm_ocn fswabs (:) = c0 flwout (:) = c0 evap (:) = c0 + fiso_evap(:,:) = c0 ! isotope evaporation rate (kg/m2/s) evaps (:) = c0 evapi (:) = c0 Tref (:) = c0 @@ -580,6 +604,7 @@ subroutine init_flux_atm_ocn fhocn (:) = c0 fswthru (:) = c0 faero_ocn(:,:) = c0 + fiso_ocn(:,:) = c0 end subroutine init_flux_atm_ocn diff --git a/configuration/driver/icedrv_forcing_bgc.F90 b/configuration/driver/icedrv_forcing_bgc.F90 index 650919417..6ad30e580 100644 --- a/configuration/driver/icedrv_forcing_bgc.F90 +++ b/configuration/driver/icedrv_forcing_bgc.F90 @@ -19,7 +19,7 @@ module icedrv_forcing_bgc implicit none private - public :: get_forcing_bgc, faero_default, init_forcing_bgc + public :: get_forcing_bgc, faero_default, fiso_default, init_forcing_bgc real (kind=dbl_kind), dimension(365) :: & ! hardwired for now sil_data, nit_data @@ -165,6 +165,23 @@ subroutine faero_default end subroutine faero_default +!======================================================================= + +! constant values for atmospheric water isotopes +! +! authors: Elizabeth Hunke, LANL + + subroutine fiso_default + + use icedrv_flux, only: fiso_atm + character(len=*), parameter :: subname='(fiso_default)' + + fiso_atm(:,1) = 1.e-12_dbl_kind ! kg/m^2 s + fiso_atm(:,2) = 1.e-13_dbl_kind + fiso_atm(:,3) = 1.e-14_dbl_kind + + end subroutine fiso_default + !======================================================================= end module icedrv_forcing_bgc diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 24f56ecfb..b760791ca 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -57,7 +57,7 @@ module icedrv_init subroutine input_data use icedrv_diagnostics, only: diag_file, nx_names - use icedrv_domain_size, only: nilyr, nslyr, nblyr, max_ntrcr, ncat, n_aero, nfsd + use icedrv_domain_size, only: nilyr, nslyr, nblyr, max_ntrcr, ncat, n_aero, n_iso, nfsd use icedrv_calendar, only: year_init, istep0 use icedrv_calendar, only: dumpfreq, diagfreq, dump_last use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years @@ -98,11 +98,11 @@ subroutine input_data logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_iso, tr_fsd logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo, wave_spec integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & - nt_aero, nt_fsd + nt_aero, nt_iso, nt_fsd real (kind=real_kind) :: rpcesm, rplvl, rptopo real (kind=dbl_kind) :: Cf, puny @@ -162,6 +162,7 @@ subroutine input_data tr_pond_lvl, & tr_pond_topo, & tr_aero, & + tr_iso, & tr_fsd !----------------------------------------------------------------- @@ -251,6 +252,7 @@ subroutine input_data tr_pond_lvl = .false. ! level-ice melt ponds tr_pond_topo = .false. ! explicit melt ponds (topographic) tr_aero = .false. ! aerosols + tr_iso = .false. ! isotopes tr_fsd = .false. ! floe size distribution !----------------------------------------------------------------- @@ -390,6 +392,13 @@ subroutine input_data call icedrv_system_abort(file=__FILE__,line=__LINE__) endif + if (tr_iso .and. n_iso==0) then + write (nu_diag,*) 'WARNING: isotopes activated but' + write (nu_diag,*) 'WARNING: not allocated in tracer array.' + write (nu_diag,*) 'WARNING: Activate in compilation script.' + call icedrv_system_abort(file=__FILE__,line=__LINE__) + endif + if (tr_aero .and. trim(shortwave) /= 'dEdd') then write (nu_diag,*) 'WARNING: aerosols activated but dEdd' write (nu_diag,*) 'WARNING: shortwave is not.' @@ -590,6 +599,7 @@ subroutine input_data write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo write(nu_diag,1010) ' tr_aero = ', tr_aero + write(nu_diag,1010) ' tr_iso = ', tr_iso write(nu_diag,1010) ' tr_fsd = ', tr_fsd nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1 @@ -655,6 +665,12 @@ subroutine input_data ntrcr = ntrcr + 4*n_aero ! 4 dEdd layers, n_aero species endif + nt_iso = max_ntrcr - 4*n_iso + if (tr_iso) then + nt_iso = ntrcr + 1 + ntrcr = ntrcr + 4*n_iso ! 4 dEdd layers, n_iso species + endif + if (ntrcr > max_ntrcr-1) then write(nu_diag,*) 'max_ntrcr-1 < number of namelist tracers' write(nu_diag,*) 'max_ntrcr-1 = ',max_ntrcr-1,' ntrcr = ',ntrcr @@ -737,9 +753,10 @@ subroutine input_data wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec) call icepack_init_tracer_sizes(ntrcr_in=ntrcr, & ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & - nfsd_in=nfsd, n_aero_in=n_aero) + nfsd_in=nfsd, n_aero_in=n_aero, n_iso_in=n_iso) call icepack_init_tracer_flags(tr_iage_in=tr_iage, & tr_FY_in=tr_FY, tr_lvl_in=tr_lvl, tr_aero_in=tr_aero, & + tr_iso_in=tr_iso, & tr_pond_in=tr_pond, tr_pond_cesm_in=tr_pond_cesm, & tr_pond_lvl_in=tr_pond_lvl, & tr_pond_topo_in=tr_pond_topo, tr_fsd_in=tr_fsd) @@ -748,7 +765,7 @@ subroutine input_data nt_qsno_in=nt_qsno, nt_iage_in=nt_iage, & nt_fy_in=nt_fy, nt_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl, & nt_apnd_in=nt_apnd, nt_hpnd_in=nt_hpnd, nt_ipnd_in=nt_ipnd, & - nt_aero_in=nt_aero, nt_fsd_in=nt_fsd) + nt_aero_in=nt_aero, nt_iso_in=nt_iso, nt_fsd_in=nt_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -815,7 +832,7 @@ end subroutine init_grid2 subroutine init_state use icepack_intfc, only: icepack_aggregate - use icedrv_domain_size, only: ncat, nilyr, nslyr, nblyr, max_ntrcr, n_aero, nfsd + use icedrv_domain_size, only: ncat, nilyr, nslyr, nblyr, max_ntrcr, n_aero, n_iso, nfsd use icedrv_flux, only: sst, Tf, Tair, salinz, Tmltz use icedrv_state, only: trcr_depend, aicen, trcrn, vicen, vsnon use icedrv_state, only: aice0, aice, vice, vsno, trcr, aice_init @@ -830,11 +847,11 @@ subroutine init_state heat_capacity ! from icepack integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_fsd logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_fy integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, & - nt_ipnd, nt_aero, nt_fsd + nt_ipnd, nt_aero, nt_iso, nt_fsd character(len=*), parameter :: subname='(init_state)' @@ -845,7 +862,7 @@ subroutine init_state call icepack_query_parameters(heat_capacity_out=heat_capacity) call icepack_query_tracer_sizes(ntrcr_out=ntrcr) call icepack_query_tracer_flags(tr_iage_out=tr_iage, & - tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, & + tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, & tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, & @@ -853,7 +870,7 @@ subroutine init_state nt_qsno_out=nt_qsno, nt_iage_out=nt_iage, nt_fy_out=nt_fy, & nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, & - nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd) + nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, nt_iso_out=nt_iso, nt_fsd_out=nt_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -936,6 +953,15 @@ subroutine init_state trcr_depend(nt_aero+(it-1)*4+3) = 1 ! ice enddo endif + if (tr_iso) then ! volume-weighted isosols + do it = 1, n_iso + trcr_depend(nt_iso+(it-1)*4 ) = 2 ! snow + trcr_depend(nt_iso+(it-1)*4+1) = 2 ! snow + trcr_depend(nt_iso+(it-1)*4+2) = 1 ! ice + trcr_depend(nt_iso+(it-1)*4+3) = 1 ! ice + enddo + endif + do it = 1, ntrcr ! mask for base quantity on which tracers are carried diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90 index 0110df0f9..5f2fcf17b 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -22,6 +22,7 @@ module icedrv_restart write_restart_lvl, read_restart_lvl, & write_restart_pond_cesm, read_restart_pond_cesm, & write_restart_pond_lvl, read_restart_pond_lvl, & + write_restart_iso, read_restart_iso, & write_restart_fsd, read_restart_fsd, & write_restart_aero, read_restart_aero @@ -61,7 +62,7 @@ subroutine dumpfile nt_Tsfc, nt_sice, nt_qice, nt_qsno logical (kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, & + tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, & tr_pond_topo, tr_pond_cesm, tr_pond_lvl, tr_fsd ! solve_zsal, skl_bgc, z_tracers @@ -79,7 +80,7 @@ subroutine dumpfile call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, tr_brine_out=tr_brine, & tr_pond_topo_out=tr_pond_topo, tr_pond_cesm_out=tr_pond_cesm, & tr_pond_lvl_out=tr_pond_lvl,tr_fsd_out=tr_fsd) ! call icepack_query_parameters(solve_zsal_out=solve_zsal, & @@ -136,7 +137,8 @@ subroutine dumpfile if (tr_pond_cesm) call write_restart_pond_cesm() ! CESM melt ponds if (tr_pond_lvl) call write_restart_pond_lvl() ! level-ice melt ponds if (tr_pond_topo) call write_restart_pond_topo() ! topographic melt ponds - if (tr_aero) call write_restart_aero() ! ice aerosol + if (tr_aero) call write_restart_aero() ! ice aerosols + if (tr_iso) call write_restart_iso() ! ice isotopes if (tr_brine) call write_restart_hbrine() ! brine height if (tr_fsd) call write_restart_fsd() ! floe size distribution ! called in icedrv_RunMod.F90 to prevent circular dependencies @@ -175,7 +177,7 @@ subroutine restartfile (ice_ic) nt_Tsfc, nt_sice, nt_qice, nt_qsno logical (kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, & + tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, & tr_pond_topo, tr_pond_cesm, tr_pond_lvl, tr_fsd character(len=char_len_long) :: filename @@ -196,7 +198,7 @@ subroutine restartfile (ice_ic) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, tr_brine_out=tr_brine, & tr_pond_topo_out=tr_pond_topo, tr_pond_cesm_out=tr_pond_cesm, & tr_pond_lvl_out=tr_pond_lvl,tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) @@ -261,7 +263,8 @@ subroutine restartfile (ice_ic) if (tr_pond_cesm) call read_restart_pond_cesm() ! CESM melt ponds if (tr_pond_lvl) call read_restart_pond_lvl() ! level-ice melt ponds if (tr_pond_topo) call read_restart_pond_topo() ! topographic melt ponds - if (tr_aero) call read_restart_aero() ! ice aerosol + if (tr_aero) call read_restart_aero() ! ice aerosols + if (tr_iso) call read_restart_iso() ! ice isotopes if (tr_brine) call read_restart_hbrine ! brine height if (tr_fsd) call read_restart_fsd() ! floe size distribution !----------------------------------------------------------------- @@ -826,6 +829,82 @@ subroutine read_restart_aero() end subroutine read_restart_aero +!======================================================================= + +! Dumps all values needed for restarting +! +! authors Elizabeth Hunke, LANL +! David Bailey, NCAR +! Marika Holland, NCAR + + subroutine write_restart_iso() + + use icedrv_domain_size, only: n_iso + use icedrv_state, only: trcrn + use icedrv_domain_size, only: ncat + + ! local variables + + integer (kind=int_kind) :: & + k ! loop indices + integer (kind=int_kind) :: nt_iso + character(len=*), parameter :: subname='(write_restart_iso)' + + call icepack_query_tracer_indices(nt_iso_out=nt_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + write(nu_diag,*) 'write_restart_iso (isotopes)' + + do k = 1, n_iso + call write_restart_field(nu_dump, trcrn(:,nt_iso +(k-1)*4,:), ncat) + call write_restart_field(nu_dump, trcrn(:,nt_iso+1+(k-1)*4,:), ncat) + call write_restart_field(nu_dump, trcrn(:,nt_iso+2+(k-1)*4,:), ncat) + call write_restart_field(nu_dump, trcrn(:,nt_iso+3+(k-1)*4,:), ncat) + enddo + + end subroutine write_restart_iso + +!======================================================================= + +! Reads all values needed for an ice isotope restart +! +! authors Elizabeth Hunke, LANL +! David Bailey, NCAR +! Marika Holland, NCAR + + subroutine read_restart_iso() + + use icedrv_domain_size, only: n_iso + use icedrv_state, only: trcrn + use icedrv_domain_size, only: ncat + + ! local variables + + integer (kind=int_kind) :: & + k ! loop indices + integer (kind=int_kind) :: nt_iso + character(len=*), parameter :: subname='(read_restart_iso)' + + !----------------------------------------------------------------- + + call icepack_query_tracer_indices(nt_iso_out=nt_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + write(nu_diag,*) 'read_restart_iso (isotopes)' + + do k = 1, n_iso + call read_restart_field(nu_restart, trcrn(:,nt_iso +(k-1)*4,:), ncat) + call read_restart_field(nu_restart, trcrn(:,nt_iso+1+(k-1)*4,:), ncat) + call read_restart_field(nu_restart, trcrn(:,nt_iso+2+(k-1)*4,:), ncat) + call read_restart_field(nu_restart, trcrn(:,nt_iso+3+(k-1)*4,:), ncat) + enddo + + end subroutine read_restart_iso + !======================================================================= subroutine write_restart_hbrine() diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 787e25c26..f9e1efc14 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -102,19 +102,20 @@ subroutine step_therm1 (dt) use icedrv_arrays_column, only: hkeel, dkeel, lfloe, dfloe use icedrv_arrays_column, only: fswsfcn, fswintn, fswthrun, Sswabsn, Iswabsn use icedrv_calendar, only: yday - use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nx + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, & fbot, Tbot, Tsnice use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm - use icedrv_flux, only: wind, rhoa, potT, Qa, zlvl, strax, stray, flatn + use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn use icedrv_flux, only: fsensn, fsurfn, fcondtopn, fcondbotn use icedrv_flux, only: flw, fsnow, fpond, sss, mlt_onset, frz_onset use icedrv_flux, only: frain, Tair, strairxT, strairyT, fsurf use icedrv_flux, only: fcondtop, fcondbot, fsens, fresh, fsalt, fhocn - use icedrv_flux, only: flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, Uref + use icedrv_flux, only: flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, Qref_iso, Uref use icedrv_flux, only: fswthru, meltt, melts, meltb, congel, snoice use icedrv_flux, only: flatn_f, fsensn_f, fsurfn_f, fcondtopn_f - use icedrv_flux, only: dsnown, faero_atm, faero_ocn + use icedrv_flux, only: dsnown, faero_atm, faero_ocn, fiso_atm, fiso_ocn, fiso_evap + use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn use icedrv_init, only: lmask_n, lmask_s use icedrv_state, only: aice, aicen, aice_init, aicen_init, vicen_init use icedrv_state, only: vice, vicen, vsno, vsnon, trcrn, uvel, vvel, vsnon_init @@ -137,15 +138,18 @@ subroutine step_therm1 (dt) integer (kind=int_kind) :: & ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, & - nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_qsno + nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_iso, nt_qsno logical (kind=log_kind) :: & - tr_iage, tr_FY, tr_aero, tr_pond, tr_pond_cesm, & + tr_iage, tr_FY, tr_aero, tr_iso, tr_pond, tr_pond_cesm, & tr_pond_lvl, tr_pond_topo, calc_Tsfc real (kind=dbl_kind), dimension(n_aero,2,ncat) :: & aerosno, aeroice ! kg/m^2 + real (kind=dbl_kind), dimension(n_iso,2,ncat) :: & + isosno, isoice ! kg/m^2 + real (kind=dbl_kind) :: & puny @@ -168,19 +172,19 @@ subroutine step_therm1 (dt) file=__FILE__,line= __LINE__) call icepack_query_tracer_flags( & - tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_aero_out=tr_aero, tr_pond_out=tr_pond, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo) + tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_aero_out=tr_aero, tr_iso_out=tr_iso, tr_pond_out=tr_pond, tr_pond_cesm_out=tr_pond_cesm, & + tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_indices( & - nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & - nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & - nt_iage_out=nt_iage, nt_FY_out=nt_FY, & - nt_qice_out=nt_qice, nt_sice_out=nt_sice, & - nt_aero_out=nt_aero, nt_qsno_out=nt_qsno) + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, & + nt_qice_out=nt_qice, nt_sice_out=nt_sice, & + nt_aero_out=nt_aero, nt_iso_out=nt_iso, nt_qsno_out=nt_qsno) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -190,6 +194,8 @@ subroutine step_therm1 (dt) prescribed_ice = .false. aerosno(:,:,:) = c0 aeroice(:,:,:) = c0 + isosno(:,:,:) = c0 + isoice(:,:,:) = c0 do i = 1, nx @@ -223,7 +229,21 @@ subroutine step_therm1 (dt) enddo enddo endif ! tr_aero - + + if (tr_iso) then + ! trcrn(nt_iso) has units kg/m^3 + do n=1,ncat + do k=1,n_iso + isosno (k,:,n) = & + trcrn(i,nt_iso+(k-1)*4 :nt_iso+(k-1)*4+1,n) & + * vsnon_init(i,n) + isoice (k,:,n) = & + trcrn(i,nt_iso+(k-1)*4+2:nt_iso+(k-1)*4+3,n) & + * vicen_init(i,n) + enddo + enddo + endif ! tr_iso + call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, & aicen_init = aicen_init(i,:), & vicen_init = vicen_init(i,:), & @@ -245,11 +265,15 @@ subroutine step_therm1 (dt) FY = trcrn(i,nt_FY,:), & aerosno = aerosno(:,:,:), & aeroice = aeroice(:,:,:), & + isosno = isosno(:,:,:), & + isoice = isoice(:,:,:), & uatm = uatm(i), vatm = vatm(i), & wind = wind(i), zlvl = zlvl(i), & Qa = Qa(i), rhoa = rhoa(i), & + Qa_iso = Qa_iso(i,:), & Tair = Tair(i), Tref = Tref(i), & Qref = Qref(i), Uref = Uref(i), & + Qref_iso = Qref_iso(i,:), & Cdn_atm_ratio = Cdn_atm_ratio(i),& Cdn_ocn = Cdn_ocn(i), & Cdn_ocn_skin = Cdn_ocn_skin(i), & @@ -289,6 +313,12 @@ subroutine step_therm1 (dt) fcondtopn_f = fcondtopn_f(i,:), & faero_atm = faero_atm(i,1:n_aero), & faero_ocn = faero_ocn(i,1:n_aero), & + fiso_atm = fiso_atm (i,:), & + fiso_ocn = fiso_ocn (i,:), & + fiso_evap = fiso_evap (i,:), & + HDO_ocn = HDO_ocn (i), & + H2_16O_ocn = H2_16O_ocn (i), & + H2_18O_ocn = H2_18O_ocn (i), & Sswabsn = Sswabsn(i,:,:),Iswabsn = Iswabsn(i,:,:), & evap = evap(i), evaps = evaps(i), evapi = evapi(i), & dhsn = dhsn(i,:), ffracn = ffracn(i,:), & @@ -317,6 +347,21 @@ subroutine step_therm1 (dt) enddo endif ! tr_aero + if (tr_iso) then + do n = 1, ncat + if (vicen(i,n) > puny) & + isoice(:,:,n) = isoice(:,:,n)/vicen(i,n) + if (vsnon(i,n) > puny) & + isosno(:,:,n) = isosno(:,:,n)/vsnon(i,n) + do k = 1, n_iso + do kk = 1, 2 + trcrn(i,nt_iso+(k-1)*4+kk-1,n)=isosno(k,kk,n) + trcrn(i,nt_iso+(k-1)*4+kk+1,n)=isoice(k,kk,n) + enddo + enddo + enddo + endif ! tr_iso + enddo ! i call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -340,11 +385,12 @@ subroutine step_therm2 (dt) d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld use icedrv_arrays_column, only: first_ice, bgrid, cgrid, igrid use icedrv_calendar, only: yday - use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, & + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nblyr, & nltrcr, nx, nfsd use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside - use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn + use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn + use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn use icedrv_init, only: tmask use icedrv_state, only: aice, aicen, aice0, trcr_depend use icedrv_state, only: aicen_init, vicen_init, trcrn, vicen, vsnon @@ -419,6 +465,10 @@ subroutine step_therm2 (dt) frazil_diag=frazil_diag(i), & frz_onset=frz_onset(i), & yday=yday, & + fiso_ocn=fiso_ocn(i,:), & + HDO_ocn=HDO_ocn(i), & + H2_16O_ocn=H2_16O_ocn(i), & + H2_18O_ocn=H2_18O_ocn(i), & nfsd=nfsd, wave_sig_ht=wave_sig_ht(i), & wave_spectrum=wave_spectrum(i,:), & wavefreq=wavefreq(:), & @@ -606,11 +656,11 @@ end subroutine step_dyn_wave subroutine step_dyn_ridge (dt, ndtd) use icedrv_arrays_column, only: hin_max, fzsal, first_ice - use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, nx + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nblyr, nx use icedrv_flux, only: rdg_conv, rdg_shear, dardg1dt, dardg2dt use icedrv_flux, only: dvirdgdt, opening, closing, fpond, fresh, fhocn use icedrv_flux, only: aparticn, krdgn, aredistn, vredistn, dardg1ndt, dardg2ndt - use icedrv_flux, only: dvirdgndt, araftn, vraftn, fsalt, flux_bio, faero_ocn + use icedrv_flux, only: dvirdgndt, araftn, vraftn, fsalt, flux_bio, faero_ocn, fiso_ocn use icedrv_init, only: tmask use icedrv_state, only: trcrn, vsnon, aicen, vicen use icedrv_state, only: aice, aice0, trcr_depend, n_trcr_strata @@ -674,8 +724,8 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgdt=dvirdgdt(i), opening=opening(i), & fpond=fpond(i), & fresh=fresh(i), fhocn=fhocn(i), & - n_aero=n_aero, & - faero_ocn=faero_ocn(i,:), & + n_aero=n_aero, n_iso=n_iso, & + faero_ocn=faero_ocn(i,:), fiso_ocn=fiso_ocn(i,:), & aparticn=aparticn(i,:), krdgn=krdgn(i,:), & aredistn=aredistn(i,:), vredistn=vredistn(i,:), & dardg1ndt=dardg1ndt(i,:), dardg2ndt=dardg2ndt(i,:), & @@ -717,8 +767,8 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgdt=dvirdgdt(i), opening=opening(i), & fpond=fpond(i), & fresh=fresh(i), fhocn=fhocn(i), & - n_aero=n_aero, & - faero_ocn=faero_ocn(i,:), & + n_aero=n_aero, n_iso=n_iso, & + faero_ocn=faero_ocn(i,:), fiso_ocn=fiso_ocn(i,:), & aparticn=aparticn(i,:), krdgn=krdgn(i,:), & aredistn=aredistn(i,:), vredistn=vredistn(i,:), & dardg1ndt=dardg1ndt(i,:), dardg2ndt=dardg2ndt(i,:), & @@ -1058,7 +1108,7 @@ subroutine biogeochemistry (dt) use icedrv_domain_size, only: n_doc, n_dic, n_don, n_fed, n_fep, nx use icedrv_flux, only: meltbn, melttn, congeln, snoicen use icedrv_flux, only: sst, sss, fsnow, meltsn - use icedrv_flux, only: hin_old, flux_bio, flux_bio_atm, faero_atm + use icedrv_flux, only: hin_old, flux_bio, flux_bio_atm, faero_atm use icedrv_flux, only: nit, amm, sil, dmsp, dms, algalN, doc, don, dic, fed, fep, zaeros, hum use icedrv_state, only: aicen_init, vicen_init, aicen, vicen, vsnon use icedrv_state, only: trcrn, vsnon_init, aice0 diff --git a/configuration/scripts/icepack.build b/configuration/scripts/icepack.build index 715b4a30a..d92456529 100755 --- a/configuration/scripts/icepack.build +++ b/configuration/scripts/icepack.build @@ -24,8 +24,8 @@ endif if !(-d ${ICE_OBJDIR}) mkdir -p ${ICE_OBJDIR} cd ${ICE_OBJDIR} -#setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNYGLOB=${ICE_NYGLOB} -DBLCKX=${ICE_BLCKX} -DBLCKY=${ICE_BLCKY} -DMXBLCKS=${ICE_MXBLCKS} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} -DNUMIN=${NUMIN} -DNUMAX=${NUMAX}" -setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " +#setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNYGLOB=${ICE_NYGLOB} -DBLCKX=${ICE_BLCKX} -DBLCKY=${ICE_BLCKY} -DMXBLCKS=${ICE_MXBLCKS} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DNTRISO=${NTRISO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} -DNUMIN=${NUMIN} -DNUMAX=${NUMAX}" +setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DNTRISO=${NTRISO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " if (${ICE_NTASKS} == 1) then setenv ICE_COMMDIR serial diff --git a/configuration/scripts/icepack.settings b/configuration/scripts/icepack.settings index e08ffdd55..6098144b2 100755 --- a/configuration/scripts/icepack.settings +++ b/configuration/scripts/icepack.settings @@ -51,6 +51,7 @@ setenv TRPND 1 # set to 1 for melt pond tracers setenv NTRAERO 1 # number of aerosol tracers # (up to max_aero in ice_domain_size.F90) # CESM uses 3 aerosol tracers +setenv NTRISO 1 # number of isotope tracers setenv TRBRI 0 # set to 1 for brine height tracer setenv TRZS 0 # set to 1 for zsalinity tracer # (needs TRBRI = 1) diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index 324267a6b..d1c0b75ba 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -28,6 +28,7 @@ tr_pond_topo = .false. tr_pond_lvl = .true. tr_aero = .false. + tr_iso = .false. tr_fsd = .false. / diff --git a/configuration/scripts/options/set_env.isotope b/configuration/scripts/options/set_env.isotope new file mode 100644 index 000000000..2ea2111a8 --- /dev/null +++ b/configuration/scripts/options/set_env.isotope @@ -0,0 +1 @@ +setenv NTRISO 3 # number of isotope tracers diff --git a/configuration/scripts/options/set_nml.isotope b/configuration/scripts/options/set_nml.isotope new file mode 100644 index 000000000..f53e29b83 --- /dev/null +++ b/configuration/scripts/options/set_nml.isotope @@ -0,0 +1 @@ +tr_iso = .true.