module water_isotopes !----------------------------------------------------------------------- ! ! Provides the functions and constants needed to calculate the isotopc flux from ! water on the ocean surface into the atmosphere. ! ! All interface routine are identified by wiso_*, etc. ! ! This code works over species indices, rather than the constituent indices ! used in the water_tracers module. As such, MAKE SURE you call these ! routines with species indicies! The tracer variable names do not need to ! match the species names, which are privided just for diagnostic output. ! ! * This module MUST be includable by CAM and CLM * (be careful with uses) ! ! ! This routine has a bunch of "qtiny" - which could be standardized. ! ! Original Code Author: David Noone <dcn@colorado.edu> - March 2003 ! ! Module added to CESM's csm_share by: Jesse Nusbaumer <nusbaume@colorado.edu> - March 2011 ! !----------------------------------------------------------------------- #undef NOFRAC /* all fractionation factors = 1 */ #undef NOKIN /* all kinetic effects off */ !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 ! use abortutils, only: endrun use shr_const_mod, only: SHR_CONST_TKTRIP ! SHR_CONST_RSTD_H2ODEV, & ! SHR_CONST_VSMOW_16O, & ! SHR_CONST_VSMOW_18O, & ! SHR_CONST_VSMOW_D , & ! SHR_CONST_VSMOW_H implicit none private save ! Public interfaces !Initialization routines: public :: wiso_init ! initilize water isotopes/tracers. public :: wiso_get_ispec ! lookup a species index by name !Fractionation routines: public :: wiso_alpl ! look-up liquid/vapor equil. fractn. public :: wiso_alpi ! look-up ice/vapor equil. fractn. public :: wiso_kmol ! kinetic effects for ocean evap (Brutsaert) public :: wiso_kmolv10 ! kmol (as above) from 10 meter wind (M&J) public :: wiso_akel ! kinetic fractionation at liq. evaporation public :: wiso_akci ! kinetic fractnation at ice condensation !Calculation routines: public :: wiso_get_roce ! retrive ocean isotope ratio. public :: wiso_flxoce ! calculate isotopic ocean evaporation. public :: wiso_ssatf ! supersaturation function public :: wiso_heff ! effective humidity function !Data checking routines: public :: wiso_get_rstd !retrive standard isotope ratio public :: wiso_get_fisub !retrive isotope subsitutions !aka number of iso. atoms per molec. public :: wiso_ratio !calculate mass ratio of isotope . public :: wiso_delta !calculate the delta value for isotopes. !configuration pointers/indices (Added from water_tracers - JN) ! integer, public :: iwspec(pcnst+pnats) ! flag for water (isotope) species ! integer, public :: ixwti, ixwtx ! lowest and highest index to search ! 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 ! Module parameters integer , parameter, public :: pwtspec = 4 ! number of water species (h2o,hdo,h218o,h216o) ! Tunable prameters for fractionation scheme real(r8), parameter :: dkfac = 0.58_r8 ! diffusive evap. kinetic power law ! real(r8), parameter :: tkini = SHR_CONST_TKTRIP ! min temp. for kinetic effects as ice appears ! real(r8), parameter :: tkini = 258.15_r8 !From Bony et. al., 2008 real(r8), parameter :: tkini = 253.15_r8 !From Jouzel and Merlivat, 1984 real(r8), parameter :: recrit = 1.0_r8 ! critical raynolds number for kmol real(r8), parameter :: fsata = 1.000_r8 ! supersaturation peramater s = a + !bTdegC (Hoffman) ! real(r8), parameter :: fsatb = -0.003_r8 ! supersaturation parameter s = a + real(r8), parameter :: fsatb = -0.002_r8 !tuned to match Antarctic d-excess in precip. - JN !bTdegC (Hoffman) real(r8), parameter :: ssatmx = 2.00_r8 ! maximum supersaturation real(r8), parameter :: fkhum = 0.25_r8 ! effective humidity factor real(r8), parameter :: tzero = SHR_CONST_TKTRIP ! supercooled water in stratiform character(len=8), dimension(pwtspec), parameter, public :: & ! species names spnam = (/ 'H2O ', 'H216O ', 'HD16O ', 'H218O ' /) ! Private isotopic constants ! ! ! Physical constants for isotopic molecules ! real(r8), dimension(pwtspec), parameter :: & ! isotopic subs. fisub = (/ 1._r8, 1._r8, 2._r8, 1._r8 /) ! TBD: Ideally this should be controlled by something like a namelist parameter, ! but it needs to be something that can be made consistent between models. real(r8), dimension(pwtspec), parameter :: & ! diffusivity ratio (note D/H, not HDO/H2O) ! difrm = (/ 1._r8, 1._r8, 0.9836504_r8, 0.9686999_r8 /) ! kinetic theory ! difrm = (/ 1._r8, 1._r8, 1._r8, 1._r8 /) ! no kinetic fractination ! difrm = (/ 1._r8, 1._r8, 0.9836504_r8, 0.9686999_r8 /) ! this with expk ! difrm = (/ 1._r8, 1._r8, 0.9755_r8, 0.9723_r8 /) ! Merlivat 1978 (tuned for isoCAM3) difrm = (/ 1._r8, 1._r8, 0.9757_r8, 0.9727_r8 /) ! Merlivat 1978 (direct from paper) ! difrm = (/ 1._r8, 1._r8, 0.9839_r8, 0.9691_r8 /) ! Cappa etal 2003 ! Prescribed isotopic ratios (largely arbitrary and tunable) real(r8), dimension(pwtspec), parameter :: & ! model standard isotope ratio !suggested by D. Noone: rstd = (/ 1._r8, 1._r8, 1._r8, 1._r8 /) ! best numerics ! rstd = (/ 1._r8, 0.5_r8, 0.25_r8, 0.2_r8, 0.1_r8 /) ! test numerics ! rstd = (/ 1._r8, 0.9976_r8, 155.76e-6_r8, 2005.20e-6_r8 /) ! natural abundance ! rstd = (/ SHR_CONST_RSTD_H2ODEV, SHR_CONST_VSMOW_16O, SHR_CONST_VSMOW_D, SHR_CONST_VSMOW_18O /) ! natural abundance ! rstd = (/ SHR_CONST_RSTD_H2ODEV, SHR_CONST_RSTD_H2ODEV, SHR_CONST_RSTD_H2ODEV, SHR_CONST_RSTD_H2ODEV /) !all 1.0 ! Isotope enrichment at ocean surface (better to be computed or read from file) real(r8), dimension(pwtspec), parameter :: & ! mean ocean surface enrichent ! boce = (/ 1._r8, 1._r8, 1.004_r8, 1.0005_r8 /) ! boce = (/ 1._r8, 1._r8, 1.0128_r8, 1.0016_r8, 1.0008_r8, 1.00671_r8 /) ! LGM boce = (/ 1._r8, 1._r8, 1._r8, 1._r8 /) ! Ocean surface kinetic fractionation parameters for M&J method: ! TBD: Check to make sure that the entries for h216o are correct. real(r8), parameter, dimension(pwtspec) :: & ! surface kinetic exchange aksmc = (/ 0._r8, 0._r8, 0.00528_r8, 0.006_r8 /), & akrfa = (/ 0._r8, 0._r8, 0.2508e-3_r8, 0.285e-3_r8 /), & akrfb = (/ 0._r8, 0._r8, 0.7216e-3_r8, 0.82e-3_r8 /) ! Coefficients for fractionation ! TBD: Check to make sure that the entries for h216o are correct. !From Majoube, 1971a: ! real(r8), parameter, dimension(pwtspec) :: & ! liquid/vapour ! alpal = (/ 0._r8, 0._r8, 24.844e+3_r8, 1.137e+3_r8 /) , & ! alpbl = (/ 0._r8, 0._r8, -76.248_r8, -0.4156_r8 /) , & ! alpcl = (/ 0._r8, 0._r8, 52.612e-3_r8, -2.0667e-3_r8 /) !From Horita and Wesolowski, 1994: real(r8), parameter, dimension(pwtspec) :: & ! liquid/vapour alpal = (/ 0._r8, 0._r8, 1158.8e-12_r8, 0.35041e+6_r8 /), & alpbl = (/ 0._r8, 0._r8, -1620.1e-9_r8, -1.6664e+3_r8 /), & alpcl = (/ 0._r8, 0._r8, 794.84e-6_r8, 6.7123_r8 /), & alpdl = (/ 0._r8, 0._r8, -161.04e-3_r8, -7.685e-3_r8 /), & alpel = (/ 0._r8, 0._r8, 2.9992e+6_r8, 0._r8 /) !isoCAM3 values: ! real(r8), parameter, dimension(pwtspec) :: & ! ice/vapour ! alpai = (/ 0._r8, 0._r8, 16288._r8, 0._r8 /), & ! alpbi = (/ 0._r8, 0._r8, 0._r8, 11.839_r8 /), & ! alpci = (/ 0._r8, 0._r8, -9.34e-2_r8, -28.224e-3_r8 /) !From Merlivat & Nief,1967 for HDO, and Majoube, 1971b for H218O: real(r8), parameter, dimension(pwtspec) :: & ! ice/vapour alpai = (/ 0._r8, 0._r8, 16289._r8, 0._r8 /), & alpbi = (/ 0._r8, 0._r8, 0._r8, 11.839_r8 /), & alpci = (/ 0._r8, 0._r8, -9.45e-2_r8, -28.224e-3_r8 /) contains !----------------------- !Initialization routines: !----------------------- !======================================================================= subroutine wiso_init !----------------------------------------------------------------------- ! Purpose: Initialize module internal data arrays ! Author: David Noone <dcn@caltech.edu> - Sun Jun 29 20:29:26 MDT 2003 !----------------------------------------------------------------------- write(6,*) 'WISO_INIT: Initializing water isotopes.' return end subroutine wiso_init !---------------------- !Fractionation routines: !---------------------- !======================================================================= subroutine wiso_kmol(isp,rbot,zbot,ustar,alpkn) !----------------------------------------------------------------------- ! ! Purpose: compute kinetic modifier for drag coefficient (Merlivat & Jouzel) ! ! Method: ! Code solves Brutsaert equations for theturbulent layer using GCM computed ! quantities. Operates on a vector of points. ! ! Author: David Noone <dcn@caltech.edu> - Mon Jun 30 14:05:38 MDT 2003 ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use shr_const_mod, only: shr_const_g, shr_const_karman implicit none real(r8), parameter :: difair = 2.36e-5_r8 ! molecular diffusivity of air real(r8), parameter :: muair = 1.7e-5_r8 ! dynamic viscosity of air ! about 17 degC, 1.73 at STP (Salby) real(r8), parameter :: gravit = shr_const_g ! gravity real(r8), parameter :: karman = shr_const_karman ! Von Karman constant !---------------------------- Arguments -------------------------------- integer , intent(in) :: isp ! species flag real(r8), intent(in) :: rbot ! density of lowest layer (kg/m3) real(r8), intent(in) :: zbot ! height of lowest level (m) real(r8), intent(in) :: ustar ! Friction velocity (m/s) ! real(r8), intent(out) :: alpkn ! kinetic fractionation factor (1-kmol) !------------------------- Local Variables ----------------------------- real(r8) z0 ! roughness length (constant in cam 9.5e-5) real(r8) reno ! surface reynolds number real(r8) tmr ! ratio of turbulen to molecular resistance real(r8) enn ! diffusive power real(r8) sc ! Schmidt number (Prandtl number) real(r8) vmu ! kinematic viscocity of air real(r8) difn ! ratio of difusivities to the power of n real(r8) difrmj ! isotopic diffusion with substitutions real(r8) kmol ! Merlivals k_mol !----------------------------------------------------------------------- ! !! difrmj = difrm(isp)/fisub(isp) difrmj = difrm(isp) ! z0 = (ustar**2._r8)/(81.1_r8*gravit) ! Charnock's equation vmu = muair / rbot ! kinematic viscosity Sc = vmu/difair reno = ustar*z0 / vmu ! reynolds number ! if (reno < recrit) then ! Smooth (Re < 0.13) enn = 2._r8/3._r8 tmr = ( (1._r8/karman)*log(ustar*zbot / (30._r8 * vmu)) ) / (13.6_r8 * Sc**(2._r8/3._r8)) else ! Rough (Re > 2) enn = 1._r8/2._r8 tmr = ( (1._r8/karman)*log(zbot/z0) - 5._r8) / (7.3_r8 * reno**(1._r8/4._r8) * Sc**(1._r8/2._r8)) end if difn = (1._r8/difrmj)**enn ! use D/Di, not Di/D kmol = (difn - 1._r8) / (difn + tmr) alpkn = 1._r8 - kmol #ifdef NOKIN ! alpkn = 1._r8 #endif ! return end subroutine wiso_kmol !======================================================================= subroutine wiso_kmolv10(isp,ustar,alpkn) !----------------------------------------------------------------------- ! ! Purpose: compute kinetic modifier for drag coefficient (Merlivat & ! Jouzel, 1979) ! ! Method: ! Uses everyones favorite empirical relation to 10 meter windspeed ! ! Author: David Noone <dcn@caltech.edu> - Mon Jun 30 14:05:38 MDT 2003 ! ! Modified U(z=10 m) calculation: Jesse Nusbaumer <nusbaume@colorado.edu> - Sept. ! 2011 ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use shr_const_mod, only: shr_const_g, shr_const_karman implicit none real(r8), parameter :: gravit = shr_const_g ! gravity real(r8), parameter :: karman = shr_const_karman ! Von Karman constant !---------------------------- Arguments -------------------------------- integer , intent(in) :: isp ! species flag real(r8), intent(in) :: ustar !friction velocity ! real(r8), intent(out) :: alpkn ! kinetic fractionation fatcor !------------------------- Local Variables ----------------------------- real(r8) z0 ! roughness length real(r8) v10 ! 10 meter winds real(r8) kmol ! Merlivat's K_mol !----------------------------------------------------------------------- ! z0 = (ustar**2._r8)/(81.1_r8*gravit) ! Charnock's equation ! v10 = ustar*log(10._r8/z0)/karman !calculate U(z=10 m) wind speed. ! ! Compute the kinetic fractionation: ! if (v10 < 7.0_r8) then ! smooth regime kmol = aksmc(isp) else ! rough regime kmol = akrfa(isp)*v10 + akrfb(isp) end if ! alpkn = 1._r8 - kmol ! !#ifdef NOKIN ! alpkn = 1.0_r8 !#endif ! return end subroutine wiso_kmolv10 !======================================================================= function wiso_alpl(isp,tk) !----------------------------------------------------------------------- ! Purpose: return liquid/vapour fractionation from look-up tables ! Author: David Noone <dcn@caltech.edu> - Mon Jun 30 10:59:13 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species indes real(r8), intent(in) :: tk ! temperature (k) real(r8) :: wiso_alpl ! return fractionation !----------------------------------------------------------------------- ! if (isp == isph2o) then wiso_alpl = 1._r8 return end if !Majoube, 1971: ! wiso_alpl = exp(alpal(isp)/tk**2 + alpbl(isp)/tk + alpcl(isp)) !Horita and Wesolowski, 1994: if(isp == isphdo) then !HDO has different formulation: wiso_alpl = exp(alpal(isp)*tk**3 + alpbl(isp)*tk**2 + alpcl(isp)*tk + alpdl(isp) + alpel(isp)/tk**3) else wiso_alpl = exp(alpal(isp)/tk**3 + alpbl(isp)/tk**2 + alpcl(isp)/tk + alpdl(isp)) end if #ifdef NOFRAC wiso_alpl = 1._r8 #endif ! return end function wiso_alpl !======================================================================= function wiso_alpi(isp,tk) !----------------------------------------------------------------------- ! Purpose: return ice/vapour fractionation from loop-up tables ! Author: David Noone <dcn@caltech.edu> - Tue Jul 1 12:02:24 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species indes real(r8), intent(in) :: tk ! temperature (k) real(r8) :: wiso_alpi ! return fractionation !----------------------------------------------------------------------- if (isp == isph2o) then wiso_alpi = 1._r8 return end if wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp)) #ifdef NOFRAC wiso_alpi = 1._r8 #endif ! return end function wiso_alpi !======================================================================= function wiso_akel(isp,tk,hum0,alpeq) !----------------------------------------------------------------------- ! Purpose: return modified fractination for kinetic effects during ! liquid evaporation into unsaturated air. ! Author: David Noone <dcn@caltech.edu> - Tue Jul 1 12:02:24 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species indes real(r8), intent(in) :: tk ! Temperature (K) real(r8), intent(in) :: hum0 ! initial humidity () real(r8), intent(in) :: alpeq ! equilibrium fractionation factor real(r8) :: wiso_akel ! return effective fractionation real(r8) :: h0 ! humidity real(r8) :: heff ! effective humidity real(r8) :: difrmj ! diffusivity for iso. sub. hum. real(r8) :: dondi ! (D / Di)^fdif, (rather than Di/D) !----------------------------------------------------------------------- !! if (tk > tkinl) then ! also do it for supercooled water h0 = min(1.0_r8,hum0) !! difrmj = difrm(isp)/fisub(isp) difrmj = difrm(isp) heff = wiso_heff(h0) dondi = (1/difrmj)**dkfac wiso_akel = alpeq*heff / (alpeq*dondi*(heff-1._r8) + 1._r8) !! else !! wiso_akel = alpeq !! end if ! ! Modify for non-standard isotope ! !! wiso_akel = wiso_akel**expk(isp) #ifdef NOKIN wiso_akel = alpeq #endif return end function wiso_akel !======================================================================= function wiso_akci(isp,tk,alpeq) !----------------------------------------------------------------------- ! Purpose: return modified fractination for kinetic effects during ! condensation to ice. ! Make use of supersaturation function. ! Author: David Noone <dcn@caltech.edu> - Tue Jul 1 12:02:24 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species indes real(r8), intent(in) :: tk ! temperature (k) real(r8), intent(in) :: alpeq ! equilibrium fractionation factor real(r8) :: wiso_akci ! return effective fractionation real(r8) :: sat1 ! super sturation real(r8) :: difrmj ! isotopic diffusion for subs. molec. real(r8) :: dondi ! D / Di, (rather than Di/D) !----------------------------------------------------------------------- ! if (tk < tkini) then ! anytime below freezing sat1 = max(1._r8, wiso_ssatf(tk)) !! difrmj = difrm(isp)/fisub(isp) difrmj = difrm(isp) dondi = 1._r8/difrmj wiso_akci = alpeq*sat1 / (alpeq*dondi*(sat1-1._r8) + 1._r8) else wiso_akci = alpeq end if ! ! Modify for non-standard isotope ! !! wiso_akci = wiso_akci**expk(isp) #ifdef NOKIN wiso_akci = alpeq #endif ! return end function wiso_akci !-------------------- !Calculation routines !-------------------- !======================================================================= function wiso_get_roce(isp) !----------------------------------------------------------------------- ! Purpose: Retrieve internal Roce variable, based on species index ! Author: David Noone <dcn@caltech.edu> - Sun Jun 29 20:29:04 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species index real(r8) :: wiso_get_roce ! return isotope ratio !----------------------------------------------------------------------- wiso_get_roce = boce(isp)*rstd(isp) return end function wiso_get_roce !======================================================================= !======================================================================= subroutine wiso_flxoce( iso ,rbot ,zbot ,wtbot , & ts , rocn, ustar ,re , & ssq, qflx, qbot, qe ) !----------------------------------------------------------------------- ! ! Purpose: compute water tracer exchange from ocean ! ! Method: ! Used diagnostics output from (./dom/)flxoce to ensure ! quantities are exactly equal for constituent number 1. ! Isotopic fractionation (equilibrium and kinetci) is applied, ! when needed. ! ! E = fac (q - qs(ts)) ! ! where fac is some exchange efficiency and qs is the saturation ! vapour mixing rati at the surface temperature. These are needed ! from calling routine to solve isotopic equivilent. ! ! Ei = fac (1-kmol) (qi - qs(ts)*Rocn/alpha) ! ! To compute the kinetic drag modifneed also ! ! Author: ! David Noone <dcn@caltech.edu> - Mon Jun 30 10:24:49 MDT 2003 ! ! Ported to CAM5, and added Schmidt, 1999 scheme - Jesse Nusbaumer <nusbaume@colorado.edu> - April, 2012 ! !----------------------------------------------------------------------- ! use shr_kind_mod, only: r8 => shr_kind_r8 ! use water_tracers, only: trace_water, wtrc_is_vap, iwspec, ixwti, ixwtx ! use water_isotopes, only: wisotope, wiso_kmol, wiso_alpl,wiso_get_roce, & ! wiso_alpi implicit none !---------------------------- Arguments -------------------------------- ! integer , intent(in) :: iso ! isotope value (1=16O,2=D,3=18O) real(r8), intent(in) :: rbot ! density of lowest layer (kg/m3) real(r8), intent(in) :: zbot ! height of lowest level (m) real(r8), intent(in) :: wtbot ! constituents at lowest real(r8), intent(in) :: qbot ! bulk water (q) at lowest real(r8), intent(in) :: qe ! bulk evaporative flux (evp) real(r8), intent(in) :: ts ! (sea) surface temperature K real(r8), intent(in) :: rocn ! (sea) surface temperature iso ratio/Rstd real(r8), intent(in) :: ustar ! friction velocity (m/s) real(r8), intent(in) :: re ! Reynolds number ? real(r8), intent(in) :: ssq ! s.hum. saturation at Ts ! real(r8), intent(out) :: qflx ! constituentflux (kg/kg/s) ! !------------------------- Local Variables ----------------------------- real(r8) alpkn ! kinetic fractionation efficiency (m) real(r8) tau ! stress real(r8) delq ! spec. hum. difference real(r8) qstar ! spec. hum,. mixing scale real(r8) Roce ! water tracer ratio of ocean surface real(r8) alpha ! fractionation factor real(r8) R_std ! tracer ratio in evaporation !----------------------------------------------------------------------- ! !-------------------------- !calculate isotopic factors !-------------------------- ! alpha = wiso_alpl(iso,ts) !get equilibrium frac. factor ! call wiso_kmolv10(iso,ustar,alpkn) !get kinetic frac. factor call wiso_kmol(iso,rbot,zbot,ustar,alpkn) !Advanced kinetic frac. routine if(rocn .eq. 0._r8) then !no ocean model data: Roce = wiso_get_roce(iso) !set to default value else !isotopic ocean model present: R_std = wiso_get_rstd(iso) !pull ratio from ocean data Roce = R_std*rocn end if !rocn value ! !----------------------------------------------- !David Noone (Merlivat and Jouzel, 1979) version !----------------------------------------------- ! ! Compute the vapour deficit then, get the fluxes ! delq = wtbot - ssq*Roce/alpha qstar = re*delq tau = rbot * ustar * ustar qflx = tau*alpkn*qstar/ustar ! !--------------------- !Schmidt, 1999 version !--------------------- ! ! rh = qbot/ssq !calculate relative humidity ! !If RH is 100%, then assume no evaporation occurs (although isotopic equilibration does, which needs to be coded in) ! ! if(rh /= 1) then ! Rate = alpkn*(Roce/alpha - (rh*wtbot/qbot))/(1-rh) !calculate ratio in flux ! else ! Rate = 0 !Assume no evaporation occurs if RH is 100% ! end if ! ! qflx = Rate*qe !convert to specific humidity (qi) return end subroutine wiso_flxoce !======================================================================= function wiso_heff(h0) !----------------------------------------------------------------------- ! Purpose: Compute effective humidity (Jouzel type thing) ! Author: David Noone <dcn@caltech.edu> - Fri Oct 24 12:06:55 PDT 2003 !----------------------------------------------------------------------- real(r8), intent(in) :: h0 ! initial humidity real(r8) :: wiso_heff ! return humidity (subsaturation) !----------------------------------------------------------------------- wiso_heff = min(1.0_r8, fkhum*h0 + 1.0_r8-fkhum) return end function wiso_heff !======================================================================= function wiso_ssatf(tk) !----------------------------------------------------------------------- ! Purpose: Compute supersaturation based on temperature parameterization. ! Author: David Noone <dcn@caltech.edu> - Sun Jun 29 20:29:14 MDT 2003 !----------------------------------------------------------------------- real(r8), intent(in) :: tk ! temperature real(r8) :: wiso_ssatf ! return supersaturation !----------------------------------------------------------------------- #ifdef OLDWAY wiso_ssatf = max(1.0_r8, fsata + fsatb*(tk-tzero)) #else wiso_ssatf = fsata + fsatb*(tk-tzero) !! wiso_ssatf = max(wiso_ssatf, fsata) wiso_ssatf = max(wiso_ssatf, 1.0_r8) wiso_ssatf = min(wiso_ssatf, ssatmx) #endif return end function wiso_ssatf !---------------------- !Data checking routines: !---------------------- !======================================================================= function wiso_get_rstd(isp) !----------------------------------------------------------------------- ! Purpose: Retrieve internal Rstd variable, based on species index ! Author: David Noone <dcn@caltech.edu> - Sun Jun 29 20:29:14 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species index real(r8) :: wiso_get_rstd ! return isotope ratio !----------------------------------------------------------------------- wiso_get_rstd = rstd(isp) return end function wiso_get_rstd !======================================================================= function wiso_get_fisub(isp) !----------------------------------------------------------------------- ! Purpose: Retrieve internal fisub variable, based on species index ! Author: David Noone <dcn@caltech.edu> - Sun Jun 29 20:28:52 MDT 2003 !----------------------------------------------------------------------- integer , intent(in) :: isp ! species index real(r8) :: wiso_get_fisub ! return number of substitutions !----------------------------------------------------------------------- wiso_get_fisub = fisub(isp) return end function wiso_get_fisub !======================================================================= function wiso_get_ispec(name) !----------------------------------------------------------------------- ! Purpose: Retrieve speciies index, based on species name ! Author: Chuck Bardeen !----------------------------------------------------------------------- character(len=*), intent(in) :: name ! species name integer :: wiso_get_ispec ! return species index !----------------------------------------------------------------------- do wiso_get_ispec = 1, pwtspec if (name == spnam(wiso_get_ispec)) then return end if end do wiso_get_ispec = ispundef return end function wiso_get_ispec !======================================================================= function wiso_ratio(isp,qiso,qtot) !----------------------------------------------------------------------- ! Purpose: Compute isotopic ratio from masses, with numerical checks ! Author David Noone <dcn@caltech.edu> - Tue Jul 1 08:32:45 MDT 2003 !----------------------------------------------------------------------- integer, intent(in) :: isp ! species index real(r8),intent(in) :: qiso ! isotopic mass real(r8),intent(in) :: qtot ! isotopic mass real(r8) :: wiso_ratio ! return value !----------------------------------------------------------------------- ! TBD: This qtiny is different than found in the equivalent routine in ! water _tracers. Also, this value is larger than the smallest support ! mixing ratios, and probably should be made smaller so as not to ! produce incorrect ratios for small values. real(r8) :: qtiny = 1.e-16_r8 !----------------------------------------------------------------------- if (qtot > 0._r8) then wiso_ratio = qiso/(qtot+qtiny) else wiso_ratio = qiso/(qtot-qtiny) end if !! wiso_ratio = espmw(isp)*wiso_ratio/fisum(isp) ! correct! end function wiso_ratio !======================================================================= function wiso_delta(isp,qiso,qtot) !----------------------------------------------------------------------- ! Purpose: Compute isotopic delta value from masses ! Author David Noone <dcn@caltech.edu> - Tue Jul 1 08:32:45 MDT 2003 !----------------------------------------------------------------------- integer, intent(in) :: isp ! species index real(r8),intent(in) :: qiso ! isotopic mass real(r8),intent(in) :: qtot ! isotopic mass real(r8) :: wiso_delta ! return value !----------------------------------------------------------------------- wiso_delta = 1000._r8 * (wiso_ratio(isp,qiso,qtot) / Rstd(isp) - 1._r8) return end function wiso_delta !========================================================================= end module water_isotopes