From 3ec02519bee3234b2b541166b9ea4cd7e1f00dec Mon Sep 17 00:00:00 2001 From: David Bailey Date: Mon, 8 Jul 2019 09:21:37 -0600 Subject: [PATCH 01/16] Water isotope branch --- columnphysics/icepack_aerosol.F90 | 19 +- columnphysics/icepack_atmo.F90 | 37 +- columnphysics/icepack_intfc.F90 | 1 + columnphysics/icepack_isotope.F90 | 729 ++++++++++++++++++++ columnphysics/icepack_itd.F90 | 105 ++- columnphysics/icepack_mechred.F90 | 62 +- columnphysics/icepack_therm_itd.F90 | 100 ++- columnphysics/icepack_therm_vertical.F90 | 34 +- columnphysics/icepack_tracers.F90 | 24 +- configuration/driver/icedrv_InitMod.F90 | 2 + configuration/driver/icedrv_domain_size.F90 | 2 + configuration/driver/icedrv_flux.F90 | 16 +- configuration/driver/icedrv_init.F90 | 45 +- configuration/driver/icedrv_restart.F90 | 93 ++- configuration/driver/icedrv_step.F90 | 78 ++- 15 files changed, 1245 insertions(+), 102 deletions(-) create mode 100644 columnphysics/icepack_isotope.F90 diff --git a/columnphysics/icepack_aerosol.F90 b/columnphysics/icepack_aerosol.F90 index af0156c60..39f461f9f 100644 --- a/columnphysics/icepack_aerosol.F90 +++ b/columnphysics/icepack_aerosol.F90 @@ -126,23 +126,18 @@ subroutine update_aerosol(dt, & ar = c1/aicen hs = vsnon*ar hi = vicen*ar - dhs_melts = -melts*ar - dhi_snoice = snoice*ar - dhs_snoice = dhi_snoice*rhoi/rhos - dhi_meltt = -meltt*ar - dhi_meltb = -meltb*ar - dhi_congel = congel*ar else ! ice disappeared during time step hs = vsnon/aice_old hi = vicen/aice_old - dhs_melts = -melts/aice_old - dhi_snoice = snoice/aice_old - dhs_snoice = dhi_snoice*rhoi/rhos - dhi_meltt = -meltt/aice_old - dhi_meltb = -meltb/aice_old - dhi_congel = congel/aice_old endif + dhs_melts = -melts + dhi_snoice = snoice + dhs_snoice = dhi_snoice*rhoi/rhos + dhi_meltt = -meltt + dhi_meltb = -meltb + dhi_congel = congel + dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice & + fsnow/rhos*dt) dhi_evap = hi - (hi_old + dhi_meltt + dhi_meltb & diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 792c018da..93fb59b8e 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -61,6 +61,8 @@ subroutine atmo_boundary_layer (sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, & + Qa_iso, Qref_iso, & uvel, vvel, & Uref ) @@ -75,6 +77,9 @@ subroutine atmo_boundary_layer (sfctype, & integer (kind=int_kind), intent(in) :: & natmiter ! number of iterations for boundary layer calculations + integer (kind=int_kind), intent(in), optional :: & + 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 +109,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(out), 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 +125,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 +147,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) @@ -173,6 +185,7 @@ subroutine atmo_boundary_layer (sfctype, & Tref = c0 Qref = c0 Uref = c0 + if (present(Qref_iso)) Qref_iso(:) = c0 delt = c0 delq = c0 shcoef = c0 @@ -356,6 +369,15 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif + if (present(Qref_iso)) then + do n = 1, n_iso + ratio = c1 + if (Qa_iso(2) > puny) & + ratio = Qa_iso(n)/Qa_iso(2) + Qref_iso(n) = Qa_iso(n) - ratio*delq*fac + enddo + endif + end subroutine atmo_boundary_layer !======================================================================= @@ -809,6 +831,8 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, & + Qa_iso, Qref_iso, & uvel, vvel, & Uref) @@ -825,6 +849,9 @@ subroutine icepack_atm_boundary(sfctype, & Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) + integer (kind=int_kind), intent(in), optional :: & + n_iso ! number of isotopes + real (kind=dbl_kind), intent(inout) :: & Cdn_atm , & ! neutral drag coefficient Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff @@ -842,6 +869,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(out), 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) @@ -890,6 +923,8 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, & + Qa_iso, Qref_iso, & worku, workv, & workr) if (icepack_warnings_aborted(subname)) return diff --git a/columnphysics/icepack_intfc.F90 b/columnphysics/icepack_intfc.F90 index 1ed3f8e0c..195ca303b 100644 --- a/columnphysics/icepack_intfc.F90 +++ b/columnphysics/icepack_intfc.F90 @@ -22,6 +22,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..88d448526 --- /dev/null +++ b/columnphysics/icepack_isotope.F90 @@ -0,0 +1,729 @@ +!======================================================================= +! +!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 +! +!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 + +!======================================================================= + + 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, & + meltt, melts, & + meltb, congel, & + snoice, evap, & + fsnow, trcrn, & + 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: c0, c1, c2, p001, p1, p5, puny + use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh + use icepack_tracers, only: nt_iso, nt_Tsfc +! +! !INPUT/OUTPUT PARAMETERS: +! + integer (kind=int_kind), intent(in) :: & + nilyr, nslyr + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + real (kind=dbl_kind), intent(in) :: & + 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) + + real (kind=dbl_kind), dimension(ntrcr), & + intent(inout) :: & + trcrn ! ice/snow tracer array + +! +! 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) + + do k=1,n_iso + isosno(k,:)=& ! isotope in snow + trcrn(nt_iso+(k-1)*4 :nt_iso+(k-1)*4+1)*vsno_old + isoice(k,:)=& ! isotope in ice + trcrn(nt_iso+(k-1)*4+2:nt_iso+(k-1)*4+3)*vice_old + isotot0(k)=isosno(k,2)+isosno(k,1)+isoice(k,2)+isoice(k,1) + enddo + +! condensation of vapor onto snow and ice + + TsfK = trcrn(nt_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 (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 (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 + + do k=1,n_iso + ! Update tracers only when vsnon/vicen is large enough. + ! Otherwise, they are unchanged from last timestep. + if (vsnon > puny) then + trcrn(nt_iso+(k-1)*4 ) = isosno(k,1)/vsnon + trcrn(nt_iso+(k-1)*4+1) = isosno(k,2)/vsnon + endif + if (vicen > puny) then + trcrn(nt_iso+(k-1)*4+2) = isoice(k,1)/vicen + trcrn(nt_iso+(k-1)*4+3) = isoice(k,2)/vicen + endif + + !do n = 1,2 + ! limit the trcrn to be positive + ! if (trcrn(nt_iso+(k-1)*4+n-1) < puny) then + ! trcrn(nt_iso+(k-1)*4+n-1) = c0 + ! fiso_ocnn(k) = fiso_ocnn(k) + & + ! trcrn(nt_iso+(k-1)*4+n-1)*vsnon/dt + ! endif + ! if (trcrn(nt_iso+(k-1)*4+n+1) < puny) then + ! trcrn(nt_iso+(k-1)*4+n+1) = c0 + ! fiso_ocnn(k) = fiso_ocnn(k) + & + ! trcrn(nt_iso+(k-1)*4+n+1)*vicen/dt + ! endif + !enddo + + enddo ! n_iso + +! 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 + +! if (trcrn(nt_iso) < -puny .or. trcrn(nt_iso+1) < -puny & +! .or. trcrn(nt_iso+2) < -puny .or. trcrn(nt_iso+3) < -puny) then +! write(nu_diag,*) 'isotope negative in isotope code' +! write(nu_diag,*) 'INT neg in isotope my_task = ',& +! my_task & +! ,' printing point i and j = ',i,j +! write(nu_diag,*) 'Int Neg iso new snowssl= ',isosno(1,1) +! write(nu_diag,*) 'Int Neg iso new snowint= ',isosno(1,2) +! write(nu_diag,*) 'Int Neg iso new ice ssl= ',isoice(1,1) +! write(nu_diag,*) 'Int Neg iso new ice int= ',isoice(1,2) +! write(nu_diag,*) 'Int Neg iso vicen = ',vice_old +! write(nu_diag,*) 'Int Neg iso vsnon = ',vsno_old +! write(nu_diag,*) 'Int Neg iso aicen = ',aicen +! write(nu_diag,*) 'Int Neg iso new vicen = ',vicen +! write(nu_diag,*) 'Int Neg iso new vsnon = ',vsnon +! write(nu_diag,*) 'Int Neg iso melts = ',melts +! write(nu_diag,*) 'Int Neg iso meltt = ',meltt +! write(nu_diag,*) 'Int Neg iso meltb = ',meltb +! write(nu_diag,*) 'Int Neg iso congel = ',congel +! write(nu_diag,*) 'Int Neg iso snoice = ',snoice +! write(nu_diag,*) 'Int Neg iso evap sno = ',evaps +! write(nu_diag,*) 'Int Neg iso evap ice = ',evapi +! write(nu_diag,*) 'Int Neg iso fsnow = ',fsnow +! write(nu_diag,*) 'Int Neg iso fiso_atm = ',fiso_atm(1) +! write(nu_diag,*) 'Int Neg iso fiso_ocnn = ',fiso_ocnn(1) +! 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) + + use icepack_parameters, only: c1 + +!----------------------------------------------------------------------- +! 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(dbl_kind), intent(in) :: tk ! temperature (k) + real(dbl_kind) :: wiso_alpi ! return fractionation +!----------------------------------------------------------------------- + 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 347858606..502e79b57 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 @@ -833,7 +836,8 @@ subroutine cleanup_itd (dt, ntrcr, & real (kind=dbl_kind), dimension (:), & intent(inout), optional :: & - faero_ocn ! aerosol flux to ocean (kg/m^2/s) + faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) + 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 @@ -855,6 +859,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) @@ -878,6 +885,7 @@ subroutine cleanup_itd (dt, ntrcr, & dfsalt = c0 dfhocn = c0 dfaero_ocn(:) = c0 + dfiso_ocn(:) = c0 dflux_bio(:) = c0 dfzsal = c0 @@ -932,7 +940,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, & @@ -940,8 +949,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 ) @@ -967,8 +978,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 !------------------------------------------------------------------- @@ -988,6 +1000,11 @@ subroutine cleanup_itd (dt, ntrcr, & faero_ocn(it) = faero_ocn(it) + dfaero_ocn(it) enddo endif + if (present(fiso_ocn)) then + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it) + enddo + endif if (present(flux_bio)) then do it = 1, nbtrcr flux_bio (it) = flux_bio(it) + dflux_bio(it) @@ -1019,7 +1036,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, & @@ -1027,8 +1045,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 ) @@ -1039,6 +1059,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) :: & @@ -1066,11 +1087,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) :: & @@ -1124,6 +1149,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& @@ -1188,8 +1221,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 @@ -1246,6 +1280,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 !----------------------------------------------------------------- @@ -1320,13 +1365,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 @@ -1347,11 +1394,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 @@ -1371,6 +1422,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 @@ -1407,13 +1467,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,11 +1501,15 @@ subroutine zap_snow_temperature(dt, ncat, & 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 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 @@ -1522,8 +1588,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 2eab9a9d0..4cdb26227 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,11 @@ subroutine ridge_ice (dt, ndtd, & faero_ocn(it) = faero_ocn(it) + maero(it)*dti enddo endif + if (present(fiso_ocn)) then + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti + enddo + endif if (present(fpond)) then fpond = fpond - mpond ! units change later endif @@ -1074,8 +1090,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 +1101,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 +1166,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 +1384,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) @@ -1709,8 +1740,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, & @@ -1729,7 +1760,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) @@ -1777,6 +1809,7 @@ subroutine icepack_step_ridge (dt, ndtd, & aredistn , & ! redistribution function: fraction of new ridge area vredistn , & ! redistribution function: fraction of new ridge volume faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) + fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean real (kind=dbl_kind), dimension(:,:), intent(inout) :: & @@ -1809,6 +1842,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, & @@ -1826,7 +1860,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1838,6 +1872,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, & @@ -1855,7 +1890,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1878,16 +1913,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 c35c3a3c4..5d74aab7e 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -30,11 +30,11 @@ 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 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 + use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, tr_iso use icepack_tracers, only: bio_index use icepack_warnings, only: warnstr, icepack_warnings_add @@ -850,9 +850,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, & aicen, vicen, & vsnon, trcrn, & @@ -868,6 +870,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) :: & @@ -896,6 +899,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) :: & @@ -975,6 +981,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_aero + !----------------------------------------------------------------- ! Biogeochemistry !----------------------------------------------------------------- @@ -1025,8 +1043,8 @@ end subroutine lateral_melt ! Adrian Turner, LANL ! subroutine add_new_ice (ncat, nilyr, nblyr, & - n_aero, dt, & - ntrcr, nltrcr, & + n_aero, n_iso, dt, & + ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & vicen, vsnon1, & @@ -1050,6 +1068,7 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & 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) :: & @@ -1337,6 +1356,33 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & 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 + + ! 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 @@ -1413,6 +1459,29 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & 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(ij))/vicen(1) + trcrn(nt_iso+3+4*(it-1),1) = & + (trcrn(nt_iso+3+4*(it-1),1)*vice1(ij) & + +frazil_conc*rhoi*vi0new(m))/vicen(1) + + fiso_ocn(it) = fiso_ocn(it) & + - frazil_conc*rhoi*vi0new(m)/dt + enddo + endif + if (tr_lvl) then alvl = trcrn(nt_alvl,1) trcrn(nt_alvl,1) = & @@ -1499,7 +1568,7 @@ end subroutine add_new_ice ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL - subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & + subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, nltrcr, & nilyr, nslyr, & hin_max, nblyr, & aicen, & @@ -1519,6 +1588,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & fhocn, update_ocn_f, & bgrid, cgrid, & igrid, faero_ocn, & + fiso_ocn, & first_ice, fzsal, & flux_bio, ocean_bio, & frazil_diag, & @@ -1530,7 +1600,8 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & 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 logical (kind=log_kind), intent(in) :: & update_ocn_f ! if true, update fresh water and salt fluxes @@ -1589,6 +1660,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & vicen , & ! volume per unit area of ice (m) vsnon , & ! volume per unit area of snow (m) faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) + fiso_ocn, & ! isotope flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean real (kind=dbl_kind), dimension(:,:), intent(inout) :: & @@ -1664,7 +1736,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & call add_new_ice (ncat, nilyr, & nblyr, & - n_aero, dt, & + n_aero, n_iso, dt, & ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & @@ -1689,9 +1761,11 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, 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, & aicen, vicen, & vsnon, trcrn, & @@ -1723,17 +1797,17 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, 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 12f50b603..2a5e22203 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -28,7 +28,7 @@ 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 + use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_iso, tr_pond use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo use icepack_therm_shared, only: ferrmax, l_brine @@ -46,6 +46,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 @@ -2015,7 +2016,7 @@ end subroutine update_state_vthermo ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL - subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & + subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & aicen_init , & vicen_init , vsnon_init , & aice , aicen , & @@ -2029,6 +2030,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & ipnd , & iage , FY , & aerosno , aeroice , & + isosno , isoice , & uatm , vatm , & wind , zlvl , & Qa , rhoa , & @@ -2071,6 +2073,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & flatn_f , fsensn_f , & fsurfn_f , fcondtopn_f , & faero_atm , faero_ocn , & + fiso_atm , fiso_ocn , & dhsn , ffracn , & meltt , melttn , & meltb , meltbn , & @@ -2086,7 +2089,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nslyr , & ! number of snow layers - n_aero ! number of aerosol tracers in use + n_aero , & ! number of aerosol tracers in use + n_iso ! number of isotope tracers in use real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -2203,6 +2207,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & fswintn , & ! SW absorbed in ice interior, below surface (W m-2) faero_atm , & ! aerosol deposition rate (kg/m^2 s) faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) + fiso_atm , & ! isotope deposition rate (kg/m^2 s) + fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond meltsn , & ! snow melt (m) @@ -2223,6 +2229,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & 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) + ! local variables integer (kind=int_kind) :: & @@ -2453,6 +2463,24 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & if (icepack_warnings_aborted(subname)) return endif + if (tr_iso) then + call update_isotope (dt, & + melttn(n), meltsn(n), & + meltbn(n), congeln(n), & + snoicen(n), evapn, & + fsnow, & + Qrefn_iso(:), trcrn(:,n), & + aicen_init(n), vicen_init(n), & + vsnon_init, & + vicen(n), vsnon(n), & + aicen(n), & + fiso_atm(:), & + fiso_evapn(:), & + fiso_ocnn(:), & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn) + + endif endif ! aicen_init !----------------------------------------------------------------- diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 850399b2d..428e09f5d 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -52,6 +52,7 @@ module icepack_tracers ntrcr , & ! number of tracers in use ntrcr_o , & ! number of non-bio tracers in use n_aero , & ! number of aerosols in use + n_iso , & ! number of isotopes in use n_zaero , & ! number of z aerosols in use n_algae , & ! number of algae in use n_doc , & ! number of DOC pools in use @@ -74,6 +75,7 @@ module icepack_tracers nt_hpnd , & ! melt pond depth nt_ipnd , & ! melt pond refrozen lid thickness nt_aero , & ! starting index for aerosols in ice + nt_iso , & ! starting index for isotopes in ice nt_bgc_Nit, & ! nutrients nt_bgc_Am, & ! nt_bgc_Sil, & ! @@ -94,6 +96,7 @@ module icepack_tracers tr_pond_lvl , & ! if .true., use level-ice pond tracer tr_pond_topo, & ! if .true., use explicit topography-based ponds tr_aero , & ! if .true., use aerosol tracers + tr_iso , & ! if .true., use isotope tracers tr_brine ! if .true., brine height differs from ice thickness !----------------------------------------------------------------- @@ -191,7 +194,8 @@ module icepack_tracers subroutine icepack_query_tracer_sizes( & max_algae_out , max_dic_out , max_doc_out , & max_don_out , max_fe_out , nmodal1_out , & - nmodal2_out , max_aero_out , max_nbtrcr_out ) + nmodal2_out , max_aero_out , & + max_nbtrcr_out ) integer (kind=int_kind), intent(out), optional :: & max_algae_out , & ! maximum number of algal types @@ -242,7 +246,7 @@ end subroutine icepack_write_tracer_sizes 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_aero_in, tr_brine_in, tr_zaero_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) @@ -256,6 +260,7 @@ subroutine icepack_init_tracer_flags(& tr_pond_lvl_in , & ! if .true., use level-ice pond tracer tr_pond_topo_in , & ! if .true., use explicit topography-based ponds 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 @@ -278,6 +283,7 @@ subroutine icepack_init_tracer_flags(& if (present(tr_pond_lvl_in) ) tr_pond_lvl = tr_pond_lvl_in if (present(tr_pond_topo_in)) tr_pond_topo = tr_pond_topo_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 @@ -300,7 +306,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_aero_out, tr_brine_out, tr_zaero_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) @@ -314,6 +320,7 @@ subroutine icepack_query_tracer_flags(& tr_pond_lvl_out , & ! if .true., use level-ice pond tracer tr_pond_topo_out , & ! if .true., use explicit topography-based ponds 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 @@ -336,6 +343,7 @@ subroutine icepack_query_tracer_flags(& if (present(tr_pond_lvl_out) ) tr_pond_lvl_out = tr_pond_lvl if (present(tr_pond_topo_out)) tr_pond_topo_out = tr_pond_topo 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 @@ -368,6 +376,7 @@ subroutine icepack_write_tracer_flags(iounit) write(iounit,*) " tr_pond_lvl = ",tr_pond_lvl write(iounit,*) " tr_pond_topo = ",tr_pond_topo 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 @@ -391,7 +400,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_aero_in, nt_zaero_in, & + nt_aero_in, nt_iso_in, nt_zaero_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, & @@ -418,6 +427,7 @@ subroutine icepack_init_tracer_indices(& nt_hpnd_in, & ! melt pond depth nt_ipnd_in, & ! melt pond refrozen lid thickness 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, & ! @@ -499,6 +509,7 @@ subroutine icepack_init_tracer_indices(& if (present(nt_hpnd_in)) nt_hpnd = nt_hpnd_in if (present(nt_ipnd_in)) nt_ipnd = nt_ipnd_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 @@ -604,7 +615,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_aero_out, nt_zaero_out, & + nt_aero_out, nt_iso_out, nt_zaero_out, & nt_bgc_N_out, nt_bgc_C_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, & @@ -630,6 +641,7 @@ subroutine icepack_query_tracer_indices(& nt_hpnd_out, & ! melt pond depth nt_ipnd_out, & ! melt pond refrozen lid thickness 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, & ! @@ -698,6 +710,7 @@ subroutine icepack_query_tracer_indices(& if (present(nt_hpnd_out)) nt_hpnd_out = nt_hpnd if (present(nt_ipnd_out)) nt_ipnd_out = nt_ipnd 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 @@ -766,6 +779,7 @@ subroutine icepack_write_tracer_indices(iounit) write(iounit,*) " nt_hpnd = ",nt_hpnd write(iounit,*) " nt_ipnd = ",nt_ipnd 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 diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index 98eff083e..b8d900b46 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -48,6 +48,7 @@ subroutine icedrv_initialize skl_bgc, & ! from icepack z_tracers, & ! from icepack tr_aero, & ! from icepack + tr_iso, & ! from icepack tr_zaero ! from icepack character(len=*), parameter :: subname='(icedrv_initialize)' @@ -97,6 +98,7 @@ subroutine icedrv_initialize call icepack_query_parameters(skl_bgc_out=skl_bgc) call icepack_query_parameters(z_tracers_out=z_tracers) 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, & diff --git a/configuration/driver/icedrv_domain_size.F90 b/configuration/driver/icedrv_domain_size.F90 index f13db9997..a97406502 100644 --- a/configuration/driver/icedrv_domain_size.F90 +++ b/configuration/driver/icedrv_domain_size.F90 @@ -19,6 +19,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 @@ -45,6 +46,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 4cbd7388e..5a3127eeb 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 @@ -176,6 +176,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), & @@ -285,6 +289,10 @@ 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 ! aerosol deposition 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 @@ -295,6 +303,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), & @@ -459,6 +471,7 @@ 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) flux_bio_atm (:,:) = c0 ! zaero and bio deposition rate (kg/m2/s) !----------------------------------------------------------------- @@ -579,6 +592,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_init.F90 b/configuration/driver/icedrv_init.F90 index 3052f0d26..da20b2c1e 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -56,7 +56,7 @@ module icedrv_init subroutine input_data use icedrv_diagnostics, only: diag_file, nx_names - use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, n_aero + use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, n_aero, n_iso 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 @@ -97,10 +97,10 @@ 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 + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_iso 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 + integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nt_iso real (kind=real_kind) :: rpcesm, rplvl, rptopo real (kind=dbl_kind) :: Cf, puny @@ -159,7 +159,8 @@ subroutine input_data tr_pond_cesm, & tr_pond_lvl, & tr_pond_topo, & - tr_aero + tr_aero, & + tr_iso !----------------------------------------------------------------- ! query Icepack values @@ -246,6 +247,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 !----------------------------------------------------------------- ! read from input file @@ -384,6 +386,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.' @@ -578,6 +587,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 nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1 ntrcr = 1 ! count tracers, starting with Tsfc = 1 @@ -636,6 +646,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 @@ -715,6 +731,7 @@ subroutine input_data call icepack_init_tracer_numbers(ntrcr_in=ntrcr) 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) @@ -723,7 +740,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_aero_in=nt_aero, nt_iso_in=nt_iso) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -790,7 +807,7 @@ end subroutine init_grid2 subroutine init_state use icepack_intfc, only: icepack_aggregate - use icedrv_domain_size, only: ncat, nilyr, nslyr, max_ntrcr, n_aero + use icedrv_domain_size, only: ncat, nilyr, nslyr, max_ntrcr, n_aero, n_iso 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 @@ -805,10 +822,10 @@ subroutine init_state heat_capacity ! from icepack integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso 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 + integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nt_iso character(len=*), parameter :: subname='(init_state)' @@ -820,6 +837,7 @@ subroutine init_state call icepack_query_tracer_numbers(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_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) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, & @@ -827,7 +845,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_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, 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__) @@ -905,6 +923,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 74a42b63e..8d53567e4 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -22,7 +22,8 @@ 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_aero, read_restart_aero + write_restart_aero, read_restart_aero, & + write_restart_iso, read_restart_iso public :: dumpfile, restartfile, final_restart, & write_restart_field, read_restart_field @@ -60,7 +61,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 ! solve_zsal, skl_bgc, z_tracers @@ -78,7 +79,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) ! call icepack_query_parameters(solve_zsal_out=solve_zsal, & @@ -135,7 +136,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 ! called in icedrv_RunMod.F90 to prevent circular dependencies ! if (solve_zsal .or. skl_bgc .or. z_tracers) & @@ -173,7 +175,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 character(len=char_len_long) :: filename @@ -194,7 +196,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) call icepack_warnings_flush(nu_diag) @@ -259,7 +261,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 !----------------------------------------------------------------- @@ -776,6 +779,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 2b7c6d7b7..adbee8c7c 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -103,7 +103,7 @@ 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, 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 @@ -114,7 +114,7 @@ subroutine step_therm1 (dt) use icedrv_flux, only: flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, 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 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 +137,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 @@ -169,7 +172,7 @@ subroutine step_therm1 (dt) 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_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, & @@ -180,7 +183,7 @@ subroutine step_therm1 (dt) 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_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 +193,8 @@ subroutine step_therm1 (dt) prescribed_ice = .false. aerosno(:,:,:) = c0 aeroice(:,:,:) = c0 + isosno(:,:,:) = c0 + isoice(:,:,:) = c0 do i = 1, nx @@ -224,6 +229,20 @@ subroutine step_therm1 (dt) 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, ncat, nilyr, nslyr, n_aero, & aicen_init (i,:), & vicen_init (i,:), vsnon_init (i,:), & @@ -243,6 +262,7 @@ subroutine step_therm1 (dt) trcrn (i,nt_iage,:), & trcrn (i,nt_FY ,:), & aerosno (:,:,:), aeroice (:,:,:), & + isosno (:,:,:), isoice (:,:,:), & uatm (i), vatm (i), & wind (i), zlvl (i), & Qa (i), rhoa (i), & @@ -286,6 +306,8 @@ subroutine step_therm1 (dt) fsurfn_f (i,:), fcondtopn_f (i,:), & faero_atm (i,1:n_aero), & faero_ocn (i,1:n_aero), & + fiso_atm (i,1:n_iso), & + fiso_ocn (i,1:n_iso), & dhsn (i,:), ffracn (i,:), & meltt (i), melttn (i,:), & meltb (i), meltbn (i,:), & @@ -312,6 +334,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, & @@ -331,10 +368,10 @@ subroutine step_therm2 (dt) use icedrv_arrays_column, only: hin_max, fzsal, ocean_bio 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, nltrcr, nx + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nblyr, nltrcr, nx use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside - 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_init, only: tmask use icedrv_state, only: aice, aicen, aice0, trcr_depend use icedrv_state, only: aicen_init, vicen_init, trcrn, vicen, vsnon @@ -371,7 +408,7 @@ subroutine step_therm2 (dt) if (tmask(i)) then - call icepack_step_therm2(dt, ncat, n_aero, nltrcr, & + call icepack_step_therm2(dt, ncat, n_aero, n_iso, nltrcr, & nilyr, nslyr, & hin_max (:), nblyr, & aicen (i,:), & @@ -390,6 +427,7 @@ subroutine step_therm2 (dt) fhocn (i), update_ocn_f, & bgrid, cgrid, & igrid, faero_ocn (i,:), & + fiso_ocn(i,:), & first_ice (i,:), fzsal (i), & flux_bio (i,1:nbtrcr), & ocean_bio (i,1:nbtrcr), & @@ -517,11 +555,11 @@ end subroutine update_state 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 @@ -583,8 +621,8 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgdt (i), opening (i), & fpond (i), & fresh (i), fhocn (i), & - n_aero, & - faero_ocn(i,:), & + n_aero, n_iso, & + faero_ocn(i,:), fiso_ocn(i,:), & aparticn (i,:), krdgn (i,:), & aredistn (i,:), vredistn (i,:), & dardg1ndt(i,:), dardg2ndt(i,:), & @@ -624,8 +662,8 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgdt (i), opening (i), & fpond (i), & fresh (i), fhocn (i), & - n_aero, & - faero_ocn(i,:), & + n_aero, n_iso, & + faero_ocn(i,:), fiso_ocn(i,:), & aparticn (i,:), krdgn (i,:), & aredistn (i,:), vredistn (i,:), & dardg1ndt(i,:), dardg2ndt(i,:), & @@ -663,7 +701,7 @@ subroutine step_radiation (dt) use icedrv_arrays_column, only: kaer_tab, waer_tab, gaer_tab, kaer_bc_tab, waer_bc_tab use icedrv_arrays_column, only: gaer_bc_tab, bcenh, swgrid, igrid use icedrv_calendar, only: calendar_type, days_per_year, nextsw_cday, yday, sec - use icedrv_domain_size, only: ncat, n_aero, nilyr, nslyr, n_zaero, n_algae, nblyr, nx + use icedrv_domain_size, only: ncat, n_aero, n_iso, nilyr, nslyr, n_zaero, n_algae, nblyr, nx use icedrv_flux, only: swvdr, swvdf, swidr, swidf, coszen, fsnow use icedrv_init, only: TLAT, TLON, tmask use icedrv_state, only: aicen, vicen, vsnon, trcrn @@ -681,7 +719,7 @@ subroutine step_radiation (dt) integer (kind=int_kind) :: & max_aero, nt_Tsfc, nt_alvl, & - nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nlt_chl_sw, & + nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nt_iso, nlt_chl_sw, & ntrcr, nbtrcr_sw, nt_fbri integer (kind=int_kind), dimension(:), allocatable :: & @@ -728,7 +766,7 @@ subroutine step_radiation (dt) call icepack_query_tracer_indices( & nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, & - nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, & + nt_iso_out=nt_iso, nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, & nt_fbri_out=nt_fbri, nt_zaero_out=nt_zaero) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -764,6 +802,7 @@ subroutine step_radiation (dt) nblyr, ntrcr, & nbtrcr_sw, nilyr, & nslyr, n_aero, & + n_iso, & n_zaero, dEdd_algae, & nlt_chl_sw, nlt_zaero_sw(:), & swgrid(:), igrid(:), & @@ -776,6 +815,7 @@ subroutine step_radiation (dt) trcrn(i,nt_hpnd,:), & trcrn(i,nt_ipnd,:), & trcrn(i,nt_aero:nt_aero+4*n_aero-1,:), & + trcrn(i,nt_iso:nt_iso+4*n_iso-1,:), & ztrcr_sw, ztrcr, & TLAT(i), TLON(i), & calendar_type, days_per_year, & @@ -971,7 +1011,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, fiso_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 From b5da658c800277baf7695fea678b20918116e3a2 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Fri, 19 Jul 2019 11:20:23 -0600 Subject: [PATCH 02/16] Water isotopes --- columnphysics/icepack_flux.F90 | 27 ++++++++++- columnphysics/icepack_isotope.F90 | 27 ++++++++--- columnphysics/icepack_therm_itd.F90 | 46 ++++++++++++++---- columnphysics/icepack_therm_vertical.F90 | 59 +++++++++++++++++++----- columnphysics/icepack_tracers.F90 | 1 + configuration/driver/icedrv_flux.F90 | 10 ++-- configuration/driver/icedrv_forcing.F90 | 6 +++ configuration/driver/icedrv_step.F90 | 40 +++++++++------- configuration/scripts/icepack.build | 4 +- configuration/scripts/icepack.settings | 1 + 10 files changed, 173 insertions(+), 48 deletions(-) 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_isotope.F90 b/columnphysics/icepack_isotope.F90 index 88d448526..cb88f09be 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -23,6 +23,7 @@ module icepack_isotope ! !USES: ! use icepack_kinds + use icepack_parameters, only: c0, c1, c2, p001, p1, p5, puny ! !EOP ! @@ -34,6 +35,15 @@ module icepack_isotope ! 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 @@ -54,6 +64,7 @@ module icepack_isotope ! subroutine update_isotope (dt, & nilyr, nslyr, & + n_iso, ntrcr, & meltt, melts, & meltb, congel, & snoice, evap, & @@ -70,14 +81,13 @@ subroutine update_isotope (dt, & ! !USES: ! ! use water_isotopes, only: wiso_alpi - use icepack_parameters, only: c0, c1, c2, p001, p1, p5, puny use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh use icepack_tracers, only: nt_iso, nt_Tsfc ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & - nilyr, nslyr + nilyr, nslyr, n_iso, ntrcr real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -697,15 +707,20 @@ end function isoice_alpha !======================================================================= function wiso_alpi(isp,tk) - use icepack_parameters, only: c1 - !----------------------------------------------------------------------- ! 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(dbl_kind), intent(in) :: tk ! temperature (k) - real(dbl_kind) :: wiso_alpi ! return fractionation + 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 diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 5d74aab7e..4808a9e4d 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -43,6 +43,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 @@ -1059,7 +1060,10 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & bgrid, cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag ) + frazil_diag, & + fiso_ocn, & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1136,6 +1140,16 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & 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) + ! local variables integer (kind=int_kind) :: & @@ -1165,6 +1179,8 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & 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 @@ -1472,13 +1488,13 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & *H2_18O_ocn trcrn(nt_iso+2+4*(it-1),1) = & - (trcrn(nt_iso+2+4*(it-1),1)*vice1(ij))/vicen(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(ij) & - +frazil_conc*rhoi*vi0new(m))/vicen(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(m)/dt + - frazil_conc*rhoi*vi0new/dt enddo endif @@ -1588,11 +1604,12 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, nltrcr, & fhocn, update_ocn_f, & bgrid, cgrid, & igrid, faero_ocn, & - fiso_ocn, & first_ice, fzsal, & flux_bio, ocean_bio, & frazil_diag, & - frz_onset, yday) + frz_onset, yday, & + fiso_ocn, HDO_ocn, & + H2_16O_ocn, H2_18O_ocn) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1660,7 +1677,6 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, nltrcr, & vicen , & ! volume per unit area of ice (m) vsnon , & ! volume per unit area of snow (m) faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) - fiso_ocn, & ! isotope flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean real (kind=dbl_kind), dimension(:,:), intent(inout) :: & @@ -1675,6 +1691,15 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, 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) + character(len=*),parameter :: subname='(icepack_step_therm2)' !----------------------------------------------------------------- @@ -1752,7 +1777,10 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, nltrcr, & cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag ) + frazil_diag, fiso_ocn, & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn) + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 2a5e22203..b8a6307f0 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -74,7 +74,7 @@ module icepack_therm_vertical ! authors: William H. Lipscomb, LANL ! C. M. Bitz, UW - subroutine thermo_vertical (nilyr, nslyr, & + subroutine thermo_vertical (nilyr, nslyr, ntrcr, & dt, aicen, & vicen, vsnon, & Tsf, zSin, & @@ -104,7 +104,8 @@ subroutine thermo_vertical (nilyr, nslyr, & integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers - nslyr ! number of snow layers + nslyr , & ! number of snow layers + ntrcr ! number of tracers in use real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2016,12 +2017,13 @@ end subroutine update_state_vthermo ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL - subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & + subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & aicen_init , & vicen_init , vsnon_init , & aice , aicen , & vice , vicen , & vsno , vsnon , & + trcrn , & uvel , vvel , & Tsfc , zqsn , & zqin , zSin , & @@ -2034,8 +2036,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & 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, & @@ -2073,7 +2077,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & flatn_f , fsensn_f , & fsurfn_f , fcondtopn_f , & faero_atm , faero_ocn , & - fiso_atm , fiso_ocn , & + fiso_atm , fiso_ocn , & + fiso_evap , & + HDO_ocn , H2_16O_ocn , & + H2_18O_ocn , & dhsn , ffracn , & meltt , melttn , & meltb , meltbn , & @@ -2089,6 +2096,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nslyr , & ! number of snow layers + ntrcr , & ! number of tracers in use n_aero , & ! number of aerosol tracers in use n_iso ! number of isotope tracers in use @@ -2178,6 +2186,18 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & 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) @@ -2207,8 +2227,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & fswintn , & ! SW absorbed in ice interior, below surface (W m-2) faero_atm , & ! aerosol deposition rate (kg/m^2 s) faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) - fiso_atm , & ! isotope deposition rate (kg/m^2 s) - fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond meltsn , & ! snow melt (m) @@ -2218,6 +2236,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & + trcrn + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zqin , & ! ice layer enthalpy (J m-3) @@ -2261,6 +2282,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & 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) @@ -2320,6 +2346,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & Trefn = c0 Qrefn = c0 + Qrefn_iso(:) = c0 + fiso_ocnn(:) = c0 + fiso_evapn(:) = c0 Urefn = c0 lhcoef = c0 shcoef = c0 @@ -2350,6 +2379,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & + n_iso, & + Qa_iso(:), Qrefn_iso(:), & uvel, vvel, & Uref=Urefn) if (icepack_warnings_aborted(subname)) return @@ -2404,7 +2435,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & if (icepack_warnings_aborted(subname)) return endif - call thermo_vertical(nilyr, nslyr, & + call thermo_vertical(nilyr, nslyr, ntrcr, & dt, aicen (n), & vicen (n), vsnon (n), & Tsfc (n), zSin (:,n), & @@ -2465,13 +2496,15 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & if (tr_iso) then call update_isotope (dt, & + nilyr, nslyr, n_iso, ntrcr, & melttn(n), meltsn(n), & meltbn(n), congeln(n), & snoicen(n), evapn, & - fsnow, & - Qrefn_iso(:), trcrn(:,n), & + fsnow, trcrn(:,n), & + Qrefn_iso(:), & + isosno(:,:,n), isoice(:,:,n), & aicen_init(n), vicen_init(n), & - vsnon_init, & + vsnon_init(n), & vicen(n), vsnon(n), & aicen(n), & fiso_atm(:), & @@ -2582,7 +2615,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & meltt, melts, & meltb, congel, & snoice, & - Uref, Urefn) + Uref, Urefn, & + Qref_iso(:),Qrefn_iso(:), & + fiso_ocn(:),fiso_ocnn(:), & + fiso_evap(:),fiso_evapn(:)) + if (icepack_warnings_aborted(subname)) return enddo ! ncat diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 428e09f5d..bd2271695 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -39,6 +39,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 diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 5a3127eeb..3d03e7d3e 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -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) @@ -291,7 +294,8 @@ module icedrv_flux real (kind=dbl_kind), & !coupling variable for tr_iso dimension (nx,icepack_max_iso), public :: & - fiso_atm ! aerosol deposition rate (kg/m^2 s) + 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 :: & diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90 index 126a78235..8b38d718a 100644 --- a/configuration/driver/icedrv_forcing.F90 +++ b/configuration/driver/icedrv_forcing.F90 @@ -19,6 +19,7 @@ module icedrv_forcing use icedrv_flux, only: zlvl, Tair, potT, rhoa, uatm, vatm, wind, & strax, stray, fsw, swvdr, swvdf, swidr, swidf, Qa, flw, frain, & fsnow, sst, sss, uocn, vocn, qdp, hmix, Tf, opening, closing, sstdat + use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn implicit none private @@ -455,6 +456,11 @@ subroutine get_forcing(timestep) endif + ! water isotopes + HDO_ocn(:) = c0 + H2_16O_ocn(:) = c0 + H2_18O_ocn(:) = c0 + end subroutine get_forcing !======================================================================= diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index adbee8c7c..0115365e7 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -106,15 +106,16 @@ subroutine step_therm1 (dt) use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, 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, fiso_atm, fiso_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 @@ -243,12 +244,13 @@ subroutine step_therm1 (dt) enddo endif ! tr_iso - call icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & + call icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & aicen_init (i,:), & vicen_init (i,:), vsnon_init (i,:), & aice (i), aicen (i,:), & vice (i), vicen (i,:), & vsno (i), vsnon (i,:), & + trcrn (i,:,:), & uvel (i), vvel (i), & trcrn (i,nt_Tsfc,:), & trcrn (i,nt_qsno:nt_qsno+nslyr-1,:), & @@ -266,8 +268,10 @@ subroutine step_therm1 (dt) uatm (i), vatm (i), & wind (i), zlvl (i), & Qa (i), rhoa (i), & + Qa_iso (i,:), & Tair (i), Tref (i), & Qref (i), Uref (i), & + Qref_iso (i,:), & Cdn_atm_ratio(i), & Cdn_ocn (i), Cdn_ocn_skin(i), & Cdn_ocn_floe(i), Cdn_ocn_keel(i), & @@ -304,10 +308,14 @@ subroutine step_therm1 (dt) fhocn (i), fswthru (i), & flatn_f (i,:), fsensn_f (i,:), & fsurfn_f (i,:), fcondtopn_f (i,:), & - faero_atm (i,1:n_aero), & - faero_ocn (i,1:n_aero), & - fiso_atm (i,1:n_iso), & - fiso_ocn (i,1:n_iso), & + faero_atm (i,1:n_aero), & + faero_ocn (i,1:n_aero), & + fiso_atm (i,1:n_iso), & + fiso_ocn (i,1:n_iso), & + fiso_evap (i,1:n_iso), & + HDO_ocn (i), & + H2_16O_ocn (i), & + H2_18O_ocn (i), & dhsn (i,:), ffracn (i,:), & meltt (i), melttn (i,:), & meltb (i), meltbn (i,:), & @@ -372,6 +380,7 @@ subroutine step_therm2 (dt) use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside 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 @@ -427,12 +436,13 @@ subroutine step_therm2 (dt) fhocn (i), update_ocn_f, & bgrid, cgrid, & igrid, faero_ocn (i,:), & - fiso_ocn(i,:), & first_ice (i,:), fzsal (i), & flux_bio (i,1:nbtrcr), & ocean_bio (i,1:nbtrcr), & frazil_diag(i), & - frz_onset (i), yday) + frz_onset (i), yday, & + fiso_ocn(i,1:n_iso), HDO_ocn(i), & + H2_16O_ocn(i), H2_18O_ocn(i)) endif ! tmask @@ -701,7 +711,7 @@ subroutine step_radiation (dt) use icedrv_arrays_column, only: kaer_tab, waer_tab, gaer_tab, kaer_bc_tab, waer_bc_tab use icedrv_arrays_column, only: gaer_bc_tab, bcenh, swgrid, igrid use icedrv_calendar, only: calendar_type, days_per_year, nextsw_cday, yday, sec - use icedrv_domain_size, only: ncat, n_aero, n_iso, nilyr, nslyr, n_zaero, n_algae, nblyr, nx + use icedrv_domain_size, only: ncat, n_aero, nilyr, nslyr, n_zaero, n_algae, nblyr, nx use icedrv_flux, only: swvdr, swvdf, swidr, swidf, coszen, fsnow use icedrv_init, only: TLAT, TLON, tmask use icedrv_state, only: aicen, vicen, vsnon, trcrn @@ -719,7 +729,7 @@ subroutine step_radiation (dt) integer (kind=int_kind) :: & max_aero, nt_Tsfc, nt_alvl, & - nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nt_iso, nlt_chl_sw, & + nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nlt_chl_sw, & ntrcr, nbtrcr_sw, nt_fbri integer (kind=int_kind), dimension(:), allocatable :: & @@ -766,7 +776,7 @@ subroutine step_radiation (dt) call icepack_query_tracer_indices( & nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, & - nt_iso_out=nt_iso, nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, & + nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, & nt_fbri_out=nt_fbri, nt_zaero_out=nt_zaero) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -802,7 +812,6 @@ subroutine step_radiation (dt) nblyr, ntrcr, & nbtrcr_sw, nilyr, & nslyr, n_aero, & - n_iso, & n_zaero, dEdd_algae, & nlt_chl_sw, nlt_zaero_sw(:), & swgrid(:), igrid(:), & @@ -815,7 +824,6 @@ subroutine step_radiation (dt) trcrn(i,nt_hpnd,:), & trcrn(i,nt_ipnd,:), & trcrn(i,nt_aero:nt_aero+4*n_aero-1,:), & - trcrn(i,nt_iso:nt_iso+4*n_iso-1,:), & ztrcr_sw, ztrcr, & TLAT(i), TLON(i), & calendar_type, days_per_year, & @@ -1011,7 +1019,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, fiso_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 cb92e941c..61e53202d 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} -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} -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} -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} -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 b1de66970..718cd138f 100755 --- a/configuration/scripts/icepack.settings +++ b/configuration/scripts/icepack.settings @@ -50,6 +50,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 3 # 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) From b812402b5644d359bdc144853a9e24609e0a9212 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Fri, 19 Jul 2019 13:52:24 -0600 Subject: [PATCH 03/16] Revert aerosol bug fix for another commit --- columnphysics/icepack_aerosol.F90 | 19 ++++++++++++------- configuration/driver/icedrv_flux.F90 | 1 + 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/columnphysics/icepack_aerosol.F90 b/columnphysics/icepack_aerosol.F90 index 9dc07b25f..45fe24f6f 100644 --- a/columnphysics/icepack_aerosol.F90 +++ b/columnphysics/icepack_aerosol.F90 @@ -126,18 +126,23 @@ subroutine update_aerosol(dt, & ar = c1/aicen hs = vsnon*ar hi = vicen*ar + dhs_melts = -melts*ar + dhi_snoice = snoice*ar + dhs_snoice = dhi_snoice*rhoi/rhos + dhi_meltt = -meltt*ar + dhi_meltb = -meltb*ar + dhi_congel = congel*ar else ! ice disappeared during time step hs = vsnon/aice_old hi = vicen/aice_old + dhs_melts = -melts/aice_old + dhi_snoice = snoice/aice_old + dhs_snoice = dhi_snoice*rhoi/rhos + dhi_meltt = -meltt/aice_old + dhi_meltb = -meltb/aice_old + dhi_congel = congel/aice_old endif - dhs_melts = -melts - dhi_snoice = snoice - dhs_snoice = dhi_snoice*rhoi/rhos - dhi_meltt = -meltt - dhi_meltb = -meltb - dhi_congel = congel - dhs_evap = hs - (hs_old + dhs_melts - dhs_snoice & + fsnow/rhos*dt) dhi_evap = hi - (hi_old + dhi_meltt + dhi_meltb & diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 3d03e7d3e..bb2de666f 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -476,6 +476,7 @@ subroutine init_coupler_flux 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) !----------------------------------------------------------------- From 63984aec00671be33a49a4e459a651112be886c6 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Fri, 19 Jul 2019 15:25:35 -0600 Subject: [PATCH 04/16] water isotopes --- configuration/driver/icedrv_flux.F90 | 5 +++++ configuration/driver/icedrv_forcing.F90 | 6 ------ configuration/scripts/icepack.settings | 2 +- configuration/scripts/icepack_in | 1 + 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index bb2de666f..09f4c45a3 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -490,6 +490,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 diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90 index 8b38d718a..126a78235 100644 --- a/configuration/driver/icedrv_forcing.F90 +++ b/configuration/driver/icedrv_forcing.F90 @@ -19,7 +19,6 @@ module icedrv_forcing use icedrv_flux, only: zlvl, Tair, potT, rhoa, uatm, vatm, wind, & strax, stray, fsw, swvdr, swvdf, swidr, swidf, Qa, flw, frain, & fsnow, sst, sss, uocn, vocn, qdp, hmix, Tf, opening, closing, sstdat - use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn implicit none private @@ -456,11 +455,6 @@ subroutine get_forcing(timestep) endif - ! water isotopes - HDO_ocn(:) = c0 - H2_16O_ocn(:) = c0 - H2_18O_ocn(:) = c0 - end subroutine get_forcing !======================================================================= diff --git a/configuration/scripts/icepack.settings b/configuration/scripts/icepack.settings index 718cd138f..611660ac0 100755 --- a/configuration/scripts/icepack.settings +++ b/configuration/scripts/icepack.settings @@ -50,7 +50,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 3 # number of isotope 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 81072bf9e..3d1647fcf 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. / &thermo_nml From 04c87e55bfc5de9db8477cf043b48a21911a67f8 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Fri, 26 Jul 2019 14:08:31 -0600 Subject: [PATCH 05/16] isotopes --- columnphysics/icepack_mechred.F90 | 8 ++++++-- columnphysics/icepack_therm_vertical.F90 | 7 +++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 4cdb26227..19a5fa8d2 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -1808,10 +1808,14 @@ subroutine icepack_step_ridge (dt, ndtd, & vraftn , & ! rafting ice volume aredistn , & ! redistribution function: fraction of new ridge area vredistn , & ! redistribution function: fraction of new ridge volume - faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) - fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean + 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) + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & trcrn ! tracers diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index b8a6307f0..c13df1021 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -74,7 +74,7 @@ module icepack_therm_vertical ! authors: William H. Lipscomb, LANL ! C. M. Bitz, UW - subroutine thermo_vertical (nilyr, nslyr, ntrcr, & + subroutine thermo_vertical (nilyr, nslyr, & dt, aicen, & vicen, vsnon, & Tsf, zSin, & @@ -104,8 +104,7 @@ subroutine thermo_vertical (nilyr, nslyr, ntrcr, & integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers - nslyr , & ! number of snow layers - ntrcr ! number of tracers in use + nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2435,7 +2434,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & if (icepack_warnings_aborted(subname)) return endif - call thermo_vertical(nilyr, nslyr, ntrcr, & + call thermo_vertical(nilyr, nslyr, & dt, aicen (n), & vicen (n), vsnon (n), & Tsfc (n), zSin (:,n), & From 4825eff9f2875c803c5d77663f81de2be66e78b0 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Wed, 28 Aug 2019 11:12:08 -0600 Subject: [PATCH 06/16] Updates to the isotopes --- columnphysics/icepack_atmo.F90 | 52 ++++++++++++---- columnphysics/icepack_isotope.F90 | 78 ++---------------------- columnphysics/icepack_therm_vertical.F90 | 19 +++--- configuration/driver/icedrv_step.F90 | 7 +-- 4 files changed, 55 insertions(+), 101 deletions(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 93fb59b8e..2075104a2 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -61,7 +61,7 @@ subroutine atmo_boundary_layer (sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - n_iso, & + n_iso, tr_iso, & Qa_iso, Qref_iso, & uvel, vvel, & Uref ) @@ -69,6 +69,9 @@ subroutine atmo_boundary_layer (sfctype, & 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 @@ -77,7 +80,7 @@ subroutine atmo_boundary_layer (sfctype, & integer (kind=int_kind), intent(in) :: & natmiter ! number of iterations for boundary layer calculations - integer (kind=int_kind), intent(in), optional :: & + integer (kind=int_kind), optional, intent(in) :: & n_iso ! number of isotopes real (kind=dbl_kind), intent(in) :: & @@ -112,7 +115,7 @@ subroutine atmo_boundary_layer (sfctype, & real (kind=dbl_kind), intent(in), optional, dimension(:) :: & Qa_iso ! specific isotopic humidity (kg/kg) - real (kind=dbl_kind), intent(out), optional, dimension(:) :: & + real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & Qref_iso ! reference specific isotopic humidity (kg/kg) real (kind=dbl_kind), intent(in) :: & @@ -185,7 +188,6 @@ subroutine atmo_boundary_layer (sfctype, & Tref = c0 Qref = c0 Uref = c0 - if (present(Qref_iso)) Qref_iso(:) = c0 delt = c0 delq = c0 shcoef = c0 @@ -369,7 +371,8 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif - if (present(Qref_iso)) then + if (present(tr_iso) .and. tr_iso) then + Qref_iso(:) = c0 do n = 1, n_iso ratio = c1 if (Qa_iso(2) > puny) & @@ -820,7 +823,7 @@ end subroutine neutral_drag_coeffs !======================================================================= - subroutine icepack_atm_boundary(sfctype, & + subroutine icepack_atm_boundary(sfctype, & Tsf, potT, & uatm, vatm, & wind, zlvl, & @@ -831,7 +834,7 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - n_iso, & + n_iso, tr_iso, & Qa_iso, Qref_iso, & uvel, vvel, & Uref) @@ -849,8 +852,11 @@ subroutine icepack_atm_boundary(sfctype, & Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) - integer (kind=int_kind), intent(in), optional :: & + 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 @@ -872,7 +878,7 @@ subroutine icepack_atm_boundary(sfctype, & real (kind=dbl_kind), intent(in), optional, dimension(:) :: & Qa_iso ! specific isotopic humidity (kg/kg) - real (kind=dbl_kind), intent(out), optional, dimension(:) :: & + real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & Qref_iso ! reference specific isotopic humidity (kg/kg) real (kind=dbl_kind), optional, intent(in) :: & @@ -885,6 +891,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 @@ -896,6 +911,13 @@ subroutine icepack_atm_boundary(sfctype, & if (present(vvel)) then workv = vvel endif + if (present(n_iso)) then + allocate(Qaiso(n_iso),Qriso(n_iso)) + Qaiso = Qa_iso + Qriso = Qref_iso + niso = n_iso + triso = tr_iso + endif Cdn_atm_ratio_n = c1 @@ -923,10 +945,10 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - n_iso, & - Qa_iso, Qref_iso, & - 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 @@ -934,6 +956,10 @@ subroutine icepack_atm_boundary(sfctype, & Uref = workr endif + if (present(Qref_iso)) then + Qref_iso = Qriso + endif + end subroutine icepack_atm_boundary !------------------------------------------------------------ diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 index cb88f09be..dab88528c 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -64,11 +64,11 @@ module icepack_isotope ! subroutine update_isotope (dt, & nilyr, nslyr, & - n_iso, ntrcr, & + n_iso, & meltt, melts, & meltb, congel, & snoice, evap, & - fsnow, trcrn, & + fsnow, Tsfc, & Qref_iso, & isosno, isoice, & aice_old, & @@ -82,17 +82,18 @@ subroutine update_isotope (dt, & ! ! use water_isotopes, only: wiso_alpi use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh - use icepack_tracers, only: nt_iso, nt_Tsfc + use icepack_tracers, only: nt_iso ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & - nilyr, nslyr, n_iso, ntrcr + nilyr, nslyr, n_iso real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), intent(in) :: & + Tsfc, & meltt, & melts, & meltb, & @@ -121,10 +122,6 @@ subroutine update_isotope (dt, & real (kind=dbl_kind), dimension(:,:), intent(inout) :: & isosno, isoice ! mass of isotopes (kg) - real (kind=dbl_kind), dimension(ntrcr), & - intent(inout) :: & - trcrn ! ice/snow tracer array - ! ! local variables ! @@ -200,17 +197,9 @@ subroutine update_isotope (dt, & fsnow/rhos*dt) evapi = hi-(hi_old-meltt-meltb+congel+snoice) - do k=1,n_iso - isosno(k,:)=& ! isotope in snow - trcrn(nt_iso+(k-1)*4 :nt_iso+(k-1)*4+1)*vsno_old - isoice(k,:)=& ! isotope in ice - trcrn(nt_iso+(k-1)*4+2:nt_iso+(k-1)*4+3)*vice_old - isotot0(k)=isosno(k,2)+isosno(k,1)+isoice(k,2)+isoice(k,1) - enddo - ! condensation of vapor onto snow and ice - TsfK = trcrn(nt_Tsfc) + Tffresh + TsfK = Tsfc + Tffresh if (evaps > c0) then ! condensation to snow do k = 1, n_iso @@ -592,67 +581,12 @@ subroutine update_isotope (dt, & ! reload tracers - do k=1,n_iso - ! Update tracers only when vsnon/vicen is large enough. - ! Otherwise, they are unchanged from last timestep. - if (vsnon > puny) then - trcrn(nt_iso+(k-1)*4 ) = isosno(k,1)/vsnon - trcrn(nt_iso+(k-1)*4+1) = isosno(k,2)/vsnon - endif - if (vicen > puny) then - trcrn(nt_iso+(k-1)*4+2) = isoice(k,1)/vicen - trcrn(nt_iso+(k-1)*4+3) = isoice(k,2)/vicen - endif - - !do n = 1,2 - ! limit the trcrn to be positive - ! if (trcrn(nt_iso+(k-1)*4+n-1) < puny) then - ! trcrn(nt_iso+(k-1)*4+n-1) = c0 - ! fiso_ocnn(k) = fiso_ocnn(k) + & - ! trcrn(nt_iso+(k-1)*4+n-1)*vsnon/dt - ! endif - ! if (trcrn(nt_iso+(k-1)*4+n+1) < puny) then - ! trcrn(nt_iso+(k-1)*4+n+1) = c0 - ! fiso_ocnn(k) = fiso_ocnn(k) + & - ! trcrn(nt_iso+(k-1)*4+n+1)*vicen/dt - ! endif - !enddo - - enddo ! n_iso - ! 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 -! if (trcrn(nt_iso) < -puny .or. trcrn(nt_iso+1) < -puny & -! .or. trcrn(nt_iso+2) < -puny .or. trcrn(nt_iso+3) < -puny) then -! write(nu_diag,*) 'isotope negative in isotope code' -! write(nu_diag,*) 'INT neg in isotope my_task = ',& -! my_task & -! ,' printing point i and j = ',i,j -! write(nu_diag,*) 'Int Neg iso new snowssl= ',isosno(1,1) -! write(nu_diag,*) 'Int Neg iso new snowint= ',isosno(1,2) -! write(nu_diag,*) 'Int Neg iso new ice ssl= ',isoice(1,1) -! write(nu_diag,*) 'Int Neg iso new ice int= ',isoice(1,2) -! write(nu_diag,*) 'Int Neg iso vicen = ',vice_old -! write(nu_diag,*) 'Int Neg iso vsnon = ',vsno_old -! write(nu_diag,*) 'Int Neg iso aicen = ',aicen -! write(nu_diag,*) 'Int Neg iso new vicen = ',vicen -! write(nu_diag,*) 'Int Neg iso new vsnon = ',vsnon -! write(nu_diag,*) 'Int Neg iso melts = ',melts -! write(nu_diag,*) 'Int Neg iso meltt = ',meltt -! write(nu_diag,*) 'Int Neg iso meltb = ',meltb -! write(nu_diag,*) 'Int Neg iso congel = ',congel -! write(nu_diag,*) 'Int Neg iso snoice = ',snoice -! write(nu_diag,*) 'Int Neg iso evap sno = ',evaps -! write(nu_diag,*) 'Int Neg iso evap ice = ',evapi -! write(nu_diag,*) 'Int Neg iso fsnow = ',fsnow -! write(nu_diag,*) 'Int Neg iso fiso_atm = ',fiso_atm(1) -! write(nu_diag,*) 'Int Neg iso fiso_ocnn = ',fiso_ocnn(1) -! endif - end subroutine update_isotope !======================================================================= diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index c13df1021..c5bdf457d 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2016,13 +2016,12 @@ end subroutine update_state_vthermo ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL - subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & + subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & aicen_init , & vicen_init , vsnon_init , & aice , aicen , & vice , vicen , & vsno , vsnon , & - trcrn , & uvel , vvel , & Tsfc , zqsn , & zqin , zSin , & @@ -2031,7 +2030,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & ipnd , & iage , FY , & aerosno , aeroice , & - isosno , isoice , & + isosno , isoice , & uatm , vatm , & wind , zlvl , & Qa , rhoa , & @@ -2095,7 +2094,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nslyr , & ! number of snow layers - ntrcr , & ! number of tracers in use n_aero , & ! number of aerosol tracers in use n_iso ! number of isotope tracers in use @@ -2235,9 +2233,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) - real (kind=dbl_kind), dimension(:,:), intent(inout) :: & - trcrn - real (kind=dbl_kind), dimension(:,:), intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zqin , & ! ice layer enthalpy (J m-3) @@ -2378,9 +2373,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - n_iso, & - Qa_iso(:), Qrefn_iso(:), & - 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 @@ -2495,11 +2490,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & if (tr_iso) then call update_isotope (dt, & - nilyr, nslyr, n_iso, ntrcr, & + nilyr, nslyr, n_iso, & melttn(n), meltsn(n), & meltbn(n), congeln(n), & snoicen(n), evapn, & - fsnow, trcrn(:,n), & + fsnow, Tsfc(n), & Qrefn_iso(:), & isosno(:,:,n), isoice(:,:,n), & aicen_init(n), vicen_init(n), & diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 0115365e7..bee2f348f 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -244,13 +244,12 @@ subroutine step_therm1 (dt) enddo endif ! tr_iso - call icepack_step_therm1(dt, ncat, nilyr, nslyr, ntrcr, n_aero, n_iso, & + call icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & aicen_init (i,:), & vicen_init (i,:), vsnon_init (i,:), & aice (i), aicen (i,:), & vice (i), vicen (i,:), & vsno (i), vsnon (i,:), & - trcrn (i,:,:), & uvel (i), vvel (i), & trcrn (i,nt_Tsfc,:), & trcrn (i,nt_qsno:nt_qsno+nslyr-1,:), & @@ -263,8 +262,8 @@ subroutine step_therm1 (dt) trcrn (i,nt_ipnd,:), & trcrn (i,nt_iage,:), & trcrn (i,nt_FY ,:), & - aerosno (:,:,:), aeroice (:,:,:), & - isosno (:,:,:), isoice (:,:,:), & + aerosno (:,:,:), aeroice (:,:,:), & + isosno (:,:,:), isoice (:,:,:), & uatm (i), vatm (i), & wind (i), zlvl (i), & Qa (i), rhoa (i), & From 3d0b962a92294a30dd504ddd97b5e584a4f543d8 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Wed, 4 Sep 2019 09:59:09 -0600 Subject: [PATCH 07/16] Latest isotope changes --- columnphysics/icepack_atmo.F90 | 18 +++++---- columnphysics/icepack_itd.F90 | 13 ++++-- columnphysics/icepack_mechred.F90 | 12 +++--- columnphysics/icepack_therm_itd.F90 | 4 +- columnphysics/icepack_therm_vertical.F90 | 8 ++-- configuration/driver/icedrv_step.F90 | 50 ++++++++++++------------ 6 files changed, 56 insertions(+), 49 deletions(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 2075104a2..fcce8b231 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -371,14 +371,16 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif - if (present(tr_iso) .and. tr_iso) then - Qref_iso(:) = c0 - do n = 1, n_iso - ratio = c1 - if (Qa_iso(2) > puny) & - ratio = Qa_iso(n)/Qa_iso(2) - Qref_iso(n) = Qa_iso(n) - ratio*delq*fac - enddo + if (present(tr_iso)) then + if (tr_iso) then + Qref_iso(:) = c0 + do n = 1, n_iso + ratio = c1 + 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 diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 502e79b57..4cacae335 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -836,7 +836,10 @@ subroutine cleanup_itd (dt, ntrcr, & real (kind=dbl_kind), dimension (:), & intent(inout), optional :: & - faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) + 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 :: & @@ -1001,9 +1004,11 @@ subroutine cleanup_itd (dt, ntrcr, & enddo endif if (present(fiso_ocn)) then - do it = 1, n_iso - fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it) - enddo + 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 diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 19a5fa8d2..36cffaaed 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -601,9 +601,11 @@ subroutine ridge_ice (dt, ndtd, & enddo endif if (present(fiso_ocn)) then - do it = 1, n_iso - fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti - enddo + 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 @@ -1808,11 +1810,9 @@ subroutine icepack_step_ridge (dt, ndtd, & vraftn , & ! rafting ice volume aredistn , & ! redistribution function: fraction of new ridge area vredistn , & ! redistribution function: fraction of new ridge volume + faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean - 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) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 4808a9e4d..7790c0586 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -992,7 +992,7 @@ subroutine lateral_melt (dt, ncat, & + trcrn(nt_iso+3+4*(k-1),n))) & * rside / dt enddo ! k - endif ! tr_aero + endif ! tr_iso !----------------------------------------------------------------- ! Biogeochemistry @@ -1372,7 +1372,7 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & enddo endif - if (tr_iso) then + if (tr_iso .and. vtmp > puny) then do it=1,n_iso if (it==1) & frazil_conc = isoice_alpha(c0,'HDO',frac) & diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index c5bdf457d..706f9e16c 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2609,10 +2609,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & meltt, melts, & meltb, congel, & snoice, & - Uref, Urefn, & - Qref_iso(:),Qrefn_iso(:), & - fiso_ocn(:),fiso_ocnn(:), & - fiso_evap(:),fiso_evapn(:)) + 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 diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index bee2f348f..270814c85 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -263,15 +263,15 @@ subroutine step_therm1 (dt) trcrn (i,nt_iage,:), & trcrn (i,nt_FY ,:), & aerosno (:,:,:), aeroice (:,:,:), & - isosno (:,:,:), isoice (:,:,:), & - uatm (i), vatm (i), & - wind (i), zlvl (i), & - Qa (i), rhoa (i), & - Qa_iso (i,:), & - Tair (i), Tref (i), & - Qref (i), Uref (i), & - Qref_iso (i,:), & - Cdn_atm_ratio(i), & + isosno (:,:,:), isoice (:,:,:), & + uatm (i), vatm (i), & + wind (i), zlvl (i), & + Qa (i), rhoa (i), & + Qa_iso (i, :), & + Tair (i), Tref (i), & + Qref (i), Uref (i), & + Qref_iso (i, :), & + Cdn_atm_ratio(i), & Cdn_ocn (i), Cdn_ocn_skin(i), & Cdn_ocn_floe(i), Cdn_ocn_keel(i), & Cdn_atm (i), Cdn_atm_skin(i), & @@ -307,14 +307,14 @@ subroutine step_therm1 (dt) fhocn (i), fswthru (i), & flatn_f (i,:), fsensn_f (i,:), & fsurfn_f (i,:), fcondtopn_f (i,:), & - faero_atm (i,1:n_aero), & - faero_ocn (i,1:n_aero), & - fiso_atm (i,1:n_iso), & - fiso_ocn (i,1:n_iso), & - fiso_evap (i,1:n_iso), & - HDO_ocn (i), & - H2_16O_ocn (i), & - H2_18O_ocn (i), & + faero_atm (i,:), & + faero_ocn (i,:), & + fiso_atm (i,:), & + fiso_ocn (i,:), & + fiso_evap (i,:), & + HDO_ocn (i), & + H2_16O_ocn (i), & + H2_18O_ocn (i), & dhsn (i,:), ffracn (i,:), & meltt (i), melttn (i,:), & meltb (i), meltbn (i,:), & @@ -416,7 +416,7 @@ subroutine step_therm2 (dt) if (tmask(i)) then - call icepack_step_therm2(dt, ncat, n_aero, n_iso, nltrcr, & + call icepack_step_therm2(dt, ncat, n_aero, n_iso, nltrcr, & nilyr, nslyr, & hin_max (:), nblyr, & aicen (i,:), & @@ -440,7 +440,7 @@ subroutine step_therm2 (dt) ocean_bio (i,1:nbtrcr), & frazil_diag(i), & frz_onset (i), yday, & - fiso_ocn(i,1:n_iso), HDO_ocn(i), & + fiso_ocn(i,:), HDO_ocn(i), & H2_16O_ocn(i), H2_18O_ocn(i)) endif ! tmask @@ -628,14 +628,14 @@ subroutine step_dyn_ridge (dt, ndtd) nt_strata(1:ntrcr,:), & dardg1dt (i), dardg2dt (i), & dvirdgdt (i), opening (i), & - fpond (i), & + fpond (i), & fresh (i), fhocn (i), & - n_aero, n_iso, & - faero_ocn(i,:), fiso_ocn(i,:), & + n_aero, n_iso, & + faero_ocn(i,:), fiso_ocn(i,:), & aparticn (i,:), krdgn (i,:), & aredistn (i,:), vredistn (i,:), & dardg1ndt(i,:), dardg2ndt(i,:), & - dvirdgndt(i,:), & + dvirdgndt(i,:), & araftn (i,:), vraftn (i,:), & aice (i), fsalt (i), & first_ice(i,:), fzsal (i), & @@ -671,8 +671,8 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgdt (i), opening (i), & fpond (i), & fresh (i), fhocn (i), & - n_aero, n_iso, & - faero_ocn(i,:), fiso_ocn(i,:), & + n_aero, n_iso, & + faero_ocn(i,:), fiso_ocn(i,:), & aparticn (i,:), krdgn (i,:), & aredistn (i,:), vredistn (i,:), & dardg1ndt(i,:), dardg2ndt(i,:), & From a5449f52072494483b7d69cdcb70618f5731937e Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 5 Sep 2019 08:31:45 -0600 Subject: [PATCH 08/16] Another isotope fix --- columnphysics/icepack_atmo.F90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 2075104a2..0bb6b1eec 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -371,14 +371,17 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif - if (present(tr_iso) .and. tr_iso) then - Qref_iso(:) = c0 - do n = 1, n_iso - ratio = c1 - if (Qa_iso(2) > puny) & - ratio = Qa_iso(n)/Qa_iso(2) - Qref_iso(n) = Qa_iso(n) - ratio*delq*fac - enddo + 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 = c1 + 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 From b2925cbc229a93a453613be5e58d4d840a75c70b Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 5 Sep 2019 14:25:19 -0600 Subject: [PATCH 09/16] One more isotope fix --- columnphysics/icepack_atmo.F90 | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 0bb6b1eec..55a3645a5 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -371,8 +371,8 @@ 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 (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 @@ -914,12 +914,15 @@ subroutine icepack_atm_boundary(sfctype, & if (present(vvel)) then workv = vvel endif - if (present(n_iso)) then + 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 @@ -963,6 +966,8 @@ subroutine icepack_atm_boundary(sfctype, & Qref_iso = Qriso endif + deallocate(Qaiso,Qriso) + end subroutine icepack_atm_boundary !------------------------------------------------------------ From 0e8b66d4f88cf731cc41252c858634e8ea0ff263 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Fri, 27 Sep 2019 15:23:27 -0600 Subject: [PATCH 10/16] Fix some bugs when Qref_iso is zero --- columnphysics/icepack_atmo.F90 | 2 +- columnphysics/icepack_isotope.F90 | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 55a3645a5..75655ed3a 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -376,7 +376,7 @@ subroutine atmo_boundary_layer (sfctype, & if (tr_iso) then Qref_iso(:) = c0 do n = 1, n_iso - ratio = c1 + ratio = c0 if (Qa_iso(2) > puny) & ratio = Qa_iso(n)/Qa_iso(2) Qref_iso(n) = Qa_iso(n) - ratio*delq*fac diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 index dab88528c..c2b30fa5e 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -207,6 +207,7 @@ subroutine update_isotope (dt, & 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) @@ -223,6 +224,7 @@ subroutine update_isotope (dt, & 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) From d95397d674425b3dd4e6543bd0c06ecf43c3af35 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Tue, 11 Feb 2020 11:50:26 -0700 Subject: [PATCH 11/16] Update from master --- columnphysics/icepack_therm_itd.F90 | 62 ++++-------------------- columnphysics/icepack_therm_vertical.F90 | 21 +------- configuration/driver/icedrv_init.F90 | 59 +++------------------- 3 files changed, 18 insertions(+), 124 deletions(-) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 2653736a1..b879d4878 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -30,20 +30,12 @@ module icepack_therm_itd use icepack_tracers, only: ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice -<<<<<<< HEAD use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_iso - use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY - 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_iso -======= - use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero 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 ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + 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 @@ -1247,16 +1239,10 @@ end subroutine lateral_melt ! Adrian Turner, LANL ! Lettie Roach, NIWA/VUW ! -<<<<<<< HEAD - subroutine add_new_ice (ncat, nilyr, nblyr, & - n_aero, n_iso, dt, & - ntrcr, nltrcr, & -======= subroutine add_new_ice (ncat, nilyr, & nfsd, nblyr, & - n_aero, dt, & + n_aero, n_iso, dt, & ntrcr, nltrcr, & ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 hin_max, ktherm, & aicen, trcrn, & vicen, vsnon1, & @@ -1272,11 +1258,9 @@ subroutine add_new_ice (ncat, nilyr, & nbtrcr, flux_bio, & ocean_bio, fzsal, & frazil_diag, & -<<<<<<< HEAD fiso_ocn, & HDO_ocn, H2_16O_ocn, & - H2_18O_ocn ) -======= + H2_18O_ocn, & wave_sig_ht, & wave_spectrum, & wavefreq, & @@ -1286,7 +1270,6 @@ subroutine add_new_ice (ncat, nilyr, & floe_rad_c, floe_binwidth) use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1364,7 +1347,6 @@ subroutine add_new_ice (ncat, nilyr, & real (kind=dbl_kind), intent(inout) :: & fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) -<<<<<<< HEAD ! water isotopes real (kind=dbl_kind), dimension(:), intent(inout) :: & @@ -1374,7 +1356,7 @@ subroutine add_new_ice (ncat, nilyr, & 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) @@ -1399,7 +1381,6 @@ subroutine add_new_ice (ncat, nilyr, & ! change in thickness distribution (area) d_afsd_latg , & ! due to fsd lateral growth d_afsd_newi ! new ice formation ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 ! local variables @@ -1943,11 +1924,7 @@ end subroutine add_new_ice ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL -<<<<<<< HEAD - subroutine icepack_step_therm2 (dt, ncat, n_aero, n_iso, nltrcr, & -======= subroutine icepack_step_therm2 (dt, ncat, nltrcr, & ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 nilyr, nslyr, & hin_max, nblyr, & aicen, & @@ -1972,10 +1949,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & flux_bio, ocean_bio, & frazil_diag, & frz_onset, yday, & -<<<<<<< HEAD fiso_ocn, HDO_ocn, & - H2_16O_ocn, H2_18O_ocn) -======= + H2_16O_ocn, H2_18O_ocn, & nfsd, wave_sig_ht, & wave_spectrum, & wavefreq, & @@ -1983,7 +1958,6 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & d_afsd_latg, d_afsd_newi, & d_afsd_latm, d_afsd_weld, & floe_rad_c, floe_binwidth) ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1991,13 +1965,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & nltrcr , & ! number of zbgc tracers nblyr , & ! number of bio layers nilyr , & ! number of ice layers -<<<<<<< HEAD - nslyr , & ! number of snow layers - n_aero , & ! number of aerosol tracers - n_iso ! number of isotope tracers -======= nslyr ! number of snow layers ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 logical (kind=log_kind), intent(in) :: & update_ocn_f ! if true, update fresh water and salt fluxes @@ -2090,7 +2058,6 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & real (kind=dbl_kind), intent(in), optional :: & yday ! day of year -<<<<<<< HEAD ! water isotopes real (kind=dbl_kind), dimension(:), intent(inout) :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) @@ -2099,9 +2066,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & 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 ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 character(len=*),parameter :: subname='(icepack_step_therm2)' @@ -2164,13 +2129,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & ! identify ice-ocean cells call add_new_ice (ncat, nilyr, & -<<<<<<< HEAD - nblyr, & - n_aero, n_iso, dt, & -======= nfsd, nblyr, & - n_aero, dt, & ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + n_aero, n_iso, dt, & ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & @@ -2186,19 +2146,15 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & -<<<<<<< HEAD frazil_diag, fiso_ocn, & HDO_ocn, H2_16O_ocn, & - H2_18O_ocn) - -======= - frazil_diag, & + H2_18O_ocn, & wave_sig_ht, & wave_spectrum, & wavefreq, dwavefreq, & d_afsd_latg, d_afsd_newi, & floe_rad_c, floe_binwidth) ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index e18f4d88b..2dae0e5d7 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -28,13 +28,9 @@ module icepack_therm_vertical use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd use icepack_parameters, only: phi_i_mushy -<<<<<<< HEAD - use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_iso, tr_pond -======= - use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + 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 @@ -2031,11 +2027,7 @@ end subroutine update_state_vthermo ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL -<<<<<<< HEAD - subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, n_iso, & -======= subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 aicen_init , & vicen_init , vsnon_init , & aice , aicen , & @@ -2113,13 +2105,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nilyr , & ! number of ice layers -<<<<<<< HEAD - nslyr , & ! number of snow layers - n_aero , & ! number of aerosol tracers in use - n_iso ! number of isotope tracers in use -======= nslyr ! number of snow layers ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -2269,13 +2255,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & aerosno , & ! snow aerosol tracer (kg/m^2) aeroice ! ice aerosol tracer (kg/m^2) -<<<<<<< HEAD real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & isosno , & ! snow isotope tracer (kg/m^2) isoice ! ice isotope tracer (kg/m^2) -======= !autodocument_end ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 ! local variables diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 116b93424..4cef4d8f8 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -754,20 +754,8 @@ subroutine input_data ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & nfsd_in=nfsd) call icepack_init_tracer_flags(tr_iage_in=tr_iage, & -<<<<<<< HEAD - 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) - call icepack_init_tracer_indices(nt_Tsfc_in=nt_Tsfc, & - nt_sice_in=nt_sice, nt_qice_in=nt_qice, & - 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_iso_in=nt_iso) -======= 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) @@ -776,12 +764,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, & -<<<<<<< HEAD - nt_aero_in=nt_aero) ->>>>>>> 68bc397d7b879e70a60a964415016e5e86cac1f4 -======= - nt_aero_in=nt_aero, nt_fsd_in=nt_fsd) ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + 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, & @@ -848,11 +831,7 @@ end subroutine init_grid2 subroutine init_state use icepack_intfc, only: icepack_aggregate -<<<<<<< HEAD - use icedrv_domain_size, only: ncat, nilyr, nslyr, max_ntrcr, n_aero, n_iso -======= - use icedrv_domain_size, only: ncat, nilyr, nslyr, nblyr, max_ntrcr, n_aero, nfsd ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + 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 @@ -867,18 +846,11 @@ subroutine init_state heat_capacity ! from icepack integer (kind=int_kind) :: ntrcr -<<<<<<< HEAD - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso - 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_iso -======= - 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 ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + nt_ipnd, nt_aero, nt_iso, nt_fsd character(len=*), parameter :: subname='(init_state)' @@ -889,19 +861,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, & -<<<<<<< HEAD - 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) - 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, 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_iso_out=nt_iso) -======= - 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, & @@ -909,12 +869,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, & -<<<<<<< HEAD - nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero) ->>>>>>> 68bc397d7b879e70a60a964415016e5e86cac1f4 -======= - nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd) ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 + 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__) From a091a6265c0ab1e7b4a685cdb5eb8903770723e2 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Tue, 18 Feb 2020 13:33:15 -0700 Subject: [PATCH 12/16] Update from master --- configuration/driver/icedrv_restart.F90 | 3 +- configuration/driver/icedrv_step.F90 | 39 ++++--------------------- 2 files changed, 7 insertions(+), 35 deletions(-) diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90 index 92099454e..5f2fcf17b 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -22,8 +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_aero, read_restart_aero, & - write_restart_iso, read_restart_iso + write_restart_iso, read_restart_iso, & write_restart_fsd, read_restart_fsd, & write_restart_aero, read_restart_aero diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 3170ee3ae..f9e1efc14 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -245,7 +245,6 @@ subroutine step_therm1 (dt) endif ! tr_iso call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, & - n_aero=n_aero, n_iso=n_iso, & aicen_init = aicen_init(i,:), & vicen_init = vicen_init(i,:), & vsnon_init = vsnon_init(i,:), & @@ -275,7 +274,6 @@ subroutine step_therm1 (dt) Tair = Tair(i), Tref = Tref(i), & Qref = Qref(i), Uref = Uref(i), & Qref_iso = Qref_iso(i,:), & - Tair = Tair(i), Tref = Tref(i), & Cdn_atm_ratio = Cdn_atm_ratio(i),& Cdn_ocn = Cdn_ocn(i), & Cdn_ocn_skin = Cdn_ocn_skin(i), & @@ -392,6 +390,7 @@ subroutine step_therm2 (dt) 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, 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 @@ -435,37 +434,7 @@ subroutine step_therm2 (dt) if (tr_fsd) & wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:))) - call icepack_step_therm2(dt, ncat, n_aero, n_iso, nltrcr, & - nilyr, nslyr, & - hin_max (:), nblyr, & - aicen (i,:), & - vicen (i,:), vsnon (i,:), & - aicen_init(i,:), vicen_init(i,:), & - trcrn (i,1:ntrcr,:), & - aice0 (i), aice (i), & - trcr_depend(1:ntrcr), trcr_base(1:ntrcr,:), & - n_trcr_strata(1:ntrcr), nt_strata(1:ntrcr,:), & - Tf (i), sss (i), & - salinz (i,:), & - rside (i), meltl (i), & - frzmlt (i), frazil (i), & - frain (i), fpond (i), & - fresh (i), fsalt (i), & - fhocn (i), update_ocn_f, & - bgrid, cgrid, & - igrid, faero_ocn (i,:), & - first_ice (i,:), fzsal (i), & - flux_bio (i,1:nbtrcr), & - ocean_bio (i,1:nbtrcr), & - frazil_diag(i), & - frz_onset (i), yday, & - fiso_ocn(i,:), HDO_ocn(i), & - H2_16O_ocn(i), H2_18O_ocn(i)) -======= - call icepack_step_therm2(dt=dt, ncat=ncat, n_aero=n_aero, & -======= call icepack_step_therm2(dt=dt, ncat=ncat, & ->>>>>>> 74839d21286a7b02bc2b29ea49a06cdc013b8683 nltrcr=nltrcr, nilyr=nilyr, nslyr=nslyr, & hin_max=hin_max(:), nblyr=nblyr, & aicen=aicen(i,:), & @@ -496,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(:), & @@ -795,7 +768,7 @@ subroutine step_dyn_ridge (dt, ndtd) fpond=fpond(i), & fresh=fresh(i), fhocn=fhocn(i), & n_aero=n_aero, n_iso=n_iso, & - faero_ocn=faero_ocn(i,:), fiso_ocn=fiso_ocn(i,:) & + 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,:), & From 48925594bda74e7b5b0321d2dc8c555ca0f791be Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 20 Feb 2020 09:07:54 -0700 Subject: [PATCH 13/16] Update from master and isotope bug fixes --- columnphysics/icepack_isotope.F90 | 3 ++ columnphysics/icepack_therm_vertical.F90 | 36 +++++++------- columnphysics/icepack_tracers.F90 | 10 +++- configuration/driver/icedrv_InitMod.F90 | 3 +- configuration/driver/icedrv_RunMod.F90 | 8 +-- configuration/driver/icedrv_diagnostics.F90 | 55 ++++++++++++++++----- configuration/driver/icedrv_flux.F90 | 1 + configuration/driver/icedrv_forcing_bgc.F90 | 19 ++++++- configuration/driver/icedrv_init.F90 | 2 +- 9 files changed, 99 insertions(+), 38 deletions(-) diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 index c2b30fa5e..64c4ec532 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -196,6 +196,7 @@ subroutine update_isotope (dt, & evaps = hs-(hs_old-melts-dhs_snoice+& fsnow/rhos*dt) evapi = hi-(hi_old-meltt-meltb+congel+snoice) + print *,'evaps,evapi',evaps,evapi,n_iso ! condensation of vapor onto snow and ice @@ -290,6 +291,7 @@ subroutine update_isotope (dt, & ! if (isosno(k,2) < c0) then ! write(nu_diag,*) 'Neg isosno(k,2)',isosno(k,2),sloss2 ! endif + print *,'sloss1,sloss2,dzssl,dzint',sloss1,sloss2,dzssl,dzint fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt enddo @@ -336,6 +338,7 @@ subroutine update_isotope (dt, & sloss2 = isoice(k,2) isoice(k,2) = c0 endif + print *,'sloss1,sloss2,dzssli,dzinti',sloss1,sloss2,dzssli,dzinti fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt enddo diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 2dae0e5d7..c05c10f64 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2503,23 +2503,25 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & endif if (tr_iso) then - call update_isotope (dt, & - nilyr, nslyr, n_iso, & - melttn(n), meltsn(n), & - meltbn(n), congeln(n), & - snoicen(n), evapn, & - fsnow, Tsfc(n), & - Qrefn_iso(:), & - isosno(:,:,n), isoice(:,:,n), & - aicen_init(n), vicen_init(n), & - vsnon_init(n), & - vicen(n), vsnon(n), & - aicen(n), & - fiso_atm(:), & - fiso_evapn(:), & - fiso_ocnn(:), & - HDO_ocn, H2_16O_ocn, & - H2_18O_ocn) + 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) + + print *,'fiso_evap',fiso_evapn(1),fiso_evapn(2),fiso_evapn(3) endif endif ! aicen_init diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 8d2bf14a2..f100b9462 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -978,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) @@ -996,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 @@ -1019,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 @@ -1036,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) @@ -1065,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 @@ -1098,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 @@ -1128,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: " @@ -1144,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 ee6922a6e..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 @@ -143,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_flux.F90 b/configuration/driver/icedrv_flux.F90 index 841be3e06..d1bff5316 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -588,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 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 4cef4d8f8..f4c340762 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -752,7 +752,7 @@ 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) + 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, & From b9b803e2eb79a16250d216ea6dedd9836d559eeb Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 20 Feb 2020 13:40:14 -0700 Subject: [PATCH 14/16] Remove print statements --- columnphysics/icepack_isotope.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 index 64c4ec532..c2b30fa5e 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -196,7 +196,6 @@ subroutine update_isotope (dt, & evaps = hs-(hs_old-melts-dhs_snoice+& fsnow/rhos*dt) evapi = hi-(hi_old-meltt-meltb+congel+snoice) - print *,'evaps,evapi',evaps,evapi,n_iso ! condensation of vapor onto snow and ice @@ -291,7 +290,6 @@ subroutine update_isotope (dt, & ! if (isosno(k,2) < c0) then ! write(nu_diag,*) 'Neg isosno(k,2)',isosno(k,2),sloss2 ! endif - print *,'sloss1,sloss2,dzssl,dzint',sloss1,sloss2,dzssl,dzint fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt enddo @@ -338,7 +336,6 @@ subroutine update_isotope (dt, & sloss2 = isoice(k,2) isoice(k,2) = c0 endif - print *,'sloss1,sloss2,dzssli,dzinti',sloss1,sloss2,dzssli,dzinti fiso_evapn(k) = fiso_evapn(k)-(sloss1+sloss2)/dt enddo From dbc1499a6ee4bc621ddeed3cb479b6cdbcc22875 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 20 Feb 2020 13:41:58 -0700 Subject: [PATCH 15/16] Remove print statements --- columnphysics/icepack_therm_vertical.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index c05c10f64..2cb99f915 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2521,8 +2521,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & HDO_ocn=HDO_ocn,H2_16O_ocn=H2_16O_ocn, & H2_18O_ocn=H2_18O_ocn) - print *,'fiso_evap',fiso_evapn(1),fiso_evapn(2),fiso_evapn(3) - endif endif ! aicen_init From 70a4dfdf027c2c15c153390dd2c47d575654410d Mon Sep 17 00:00:00 2001 From: David Bailey Date: Wed, 26 Feb 2020 13:38:30 -0700 Subject: [PATCH 16/16] Add test settings for isotopes --- configuration/scripts/options/set_env.isotope | 1 + configuration/scripts/options/set_nml.isotope | 1 + 2 files changed, 2 insertions(+) create mode 100644 configuration/scripts/options/set_env.isotope create mode 100644 configuration/scripts/options/set_nml.isotope 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.