atmmoiststab_module Module Reference

Functions/Subroutines
subroutine	cal_atmmoist (temp_c, press_hpa, avrh, dectime, lv_j_kg, lvs_j_kg, es_hpa, ea_hpa, vpd_hpa, vpd_pa, dq, dens_dry, avcp, air_dens)
subroutine	cal_stab (stabilitymethod, zzd, z0m, zdm, avu1, temp_c, qh_init, avdens, avcp, l_mod, tstar, ustar, zl)
real(kind(1d0)) function	stab_psi_mom (stabilitymethod, zl)
real(kind(1d0)) function	stab_psi_heat (stabilitymethod, zl)
real(kind(1d0)) function	stab_phi_mom (stabilitymethod, zl)
real(kind(1d0)) function	stab_phi_heat (stabilitymethod, zl)
real(kind(1d0)) function	psi_mom_j12 (zl)
real(kind(1d0)) function	phi_mom_j12 (zl)
real(kind(1d0)) function	psi_heat_j12 (zl)
real(kind(1d0)) function	phi_heat_j12 (zl)
real(kind(1d0)) function	psi_mom_g00 (zl)
real(kind(1d0)) function	psi_heat_g00 (zl)
real(kind(1d0)) function	phi_mom_g00 (zl)
real(kind(1d0)) function	phi_heat_g00 (zl)
real(kind(1d0)) function	psi_conv (zl, ax)
real(kind(1d0)) function	phi_conv (zl, ax)
real(kind(1d0)) function	dpsi_dzl_g00 (zl, psik, phik, psic, phic)
real(kind(1d0)) function	psi_cb05 (zl, k1, k2)
real(kind(1d0)) function	psi_mom_cb05 (zl)
real(kind(1d0)) function	psi_heat_cb05 (zl)
real(kind(1d0)) function	phi_cb05 (zl, k1, k2)
real(kind(1d0)) function	phi_mom_cb05 (zl)
real(kind(1d0)) function	phi_heat_cb05 (zl)
real(kind(1d0)) function	phi_mom_k75 (zl)
real(kind(1d0)) function	phi_heat_k75 (zl)
real(kind(1d0)) function	psi_mom_k75 (zl)
real(kind(1d0)) function	psi_heat_k75 (zl)
real(kind(1d0)) function	phi_mom_b71 (zl)
real(kind(1d0)) function	phi_heat_b71 (zl)
real(kind(1d0)) function	psi_mom_b71 (zl)
real(kind(1d0)) function	psi_heat_b71 (zl)

Variables
real(kind(1d0)), parameter	neut_limit = 1.E-2
real(kind(1d0)), parameter	k = 0.4
real(kind(1d0)), parameter	grav = 9.80665
integer, parameter	j12 = 2
integer, parameter	k75 = 3
integer, parameter	b71 = 4

Function/Subroutine Documentation

cal_atmmoist()

subroutine atmmoiststab_module::cal_atmmoist	(	real(kind(1d0)), intent(in)	temp_c,
		real(kind(1d0)), intent(in)	press_hpa,
		real(kind(1d0)), intent(in)	avrh,
		real(kind(1d0)), intent(in)	dectime,
		real(kind(1d0)), intent(out)	lv_j_kg,
		real(kind(1d0)), intent(out)	lvs_j_kg,
		real(kind(1d0)), intent(out)	es_hpa,
		real(kind(1d0)), intent(out)	ea_hpa,
		real(kind(1d0)), intent(out)	vpd_hpa,
		real(kind(1d0)), intent(out)	vpd_pa,
		real(kind(1d0)), intent(out)	dq,
		real(kind(1d0)), intent(out)	dens_dry,
		real(kind(1d0)), intent(out)	avcp,
		real(kind(1d0)), intent(out)	air_dens )

Definition at line 16 of file suews_phys_atmmoiststab.f95.

 
      USE meteo, ONLY: &
         sat_vap_press_x, spec_heat_beer, &
         lat_vap, lat_vapsublim, spec_hum_def
 
      IMPLICIT NONE
      REAL(KIND(1D0)) :: vap_dens
 
      REAL(KIND(1D0)), INTENT(in) :: &
         temp_c, &
         press_hpa, &
         avrh, dectime
      REAL(KIND(1D0)), INTENT(out) :: &
         lv_j_kg, & !Latent heat of vaporization in [J kg-1]
         lvs_j_kg, & !Latent heat of sublimation in J/kg
         es_hpa, & !Saturation vapour pressure over water in hPa
         ea_hpa, & !Vapour pressure of water in hPa
         vpd_hpa, & !vapour pressure deficit in hPa
         vpd_pa, & !vapour pressure deficit in Pa
         dq, & !Specific humidity deficit in g/kg
         dens_dry, & !Vap density or absolute humidity [kg m-3]
         avcp, & !specific heat capacity in J kg-1 K-1
         air_dens !Air density in kg/m3
 
      REAL(KIND(1D0)), PARAMETER :: &
         !  comp          = 0.9995, &
         !  epsil         = 0.62197,&           !ratio molecular weight of water vapor/dry air (kg/mol/kg/mol)
         !  epsil_gkg     = 621.97, &           !ratio molecular weight of water vapor/dry air in g/kg
         !  dry_gas       = 8.31451,&           !Dry gas constant (J/k/mol)
         !  gas_ct_wat    = 461.05,&            !Gas constant for water (J/kg/K)
         !  molar         = 0.028965,&          !Dry air molar fraction in kg/mol
         !  molar_wat_vap = 0.0180153,&         !Molar fraction of water vapor in kg/mol
         gas_ct_dry = 8.31451/0.028965, & !j/kg/k=dry_gas/molar
         gas_ct_wv = 8.31451/0.0180153 !j/kg/k=dry_gas/molar_wat_vap
      !  waterDens     = 999.8395            !Density of water in 0 cel deg
      INTEGER :: from = 1
 
      !Saturation vapour pressure over water in hPa
      es_hpa = sat_vap_press_x(temp_c, press_hpa, from, dectime) ! dectime is more or less unnecessary here
 
      !Vapour pressure of water in hPa
      ea_hpa = avrh/100*es_hpa
 
      ! if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'%',Temp_C
 
      vpd_hpa = es_hpa - ea_hpa !vapour pressure deficit in hPa
      vpd_pa = (es_hpa*100) - (ea_hpa*100) !vapour pressure deficit in Pa
 
      dq = spec_hum_def(vpd_hpa, press_hpa) !Specific humidity deficit in g/kg
 
      !Vap density or absolute humidity         (kg/m3)
      vap_dens = (ea_hpa*100/((temp_c + 273.16)*gas_ct_wv))
 
      !density Dry Air Beer(1990)        kg/m3
      dens_dry = ((press_hpa - ea_hpa)*100)/(gas_ct_dry*(273.16 + temp_c))
 
      !Air density in kg/m3
      air_dens = (press_hpa*100)/(gas_ct_dry*(temp_c + 273.16))
 
      !Calculate specific heat capacity in J kg-1 K-1
      avcp = spec_heat_beer(temp_c, avrh, vap_dens, dens_dry)
 
      !Latent heat of vaporization in [J kg-1]
      lv_j_kg = lat_vap(temp_c, ea_hpa, press_hpa, avcp, dectime)
 
      !Latent heat of sublimation in J/kg
      IF (temp_c < 0.000) THEN
         lvs_j_kg = lat_vapsublim(temp_c, ea_hpa, press_hpa, avcp)
      ELSE
         lvs_j_kg = 2834000
      END IF
      !if(debug)write(*,*)lv_J_kg,Temp_C,'lv2'
      IF (press_hpa < 900) THEN
         CALL errorhint(46, 'Function LUMPS_cal_AtmMoist', press_hpa, -55.55d0, -55)
      END IF
      RETURN

References errorhint(), meteo::lat_vap(), meteo::lat_vapsublim(), meteo::sat_vap_press_x(), meteo::spec_heat_beer(), and meteo::spec_hum_def().

Referenced by suews_driver::suews_cal_biogenco2(), suews_driver::suews_cal_biogenco2_dts(), suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the call graph for this function:

Here is the caller graph for this function:

cal_stab()

subroutine atmmoiststab_module::cal_stab	(	integer, intent(in)	stabilitymethod,
		real(kind(1d0)), intent(in)	zzd,
		real(kind(1d0)), intent(in)	z0m,
		real(kind(1d0)), intent(in)	zdm,
		real(kind(1d0)), intent(in)	avu1,
		real(kind(1d0)), intent(in)	temp_c,
		real(kind(1d0)), intent(in)	qh_init,
		real(kind(1d0)), intent(in)	avdens,
		real(kind(1d0)), intent(in)	avcp,
		real(kind(1d0)), intent(out)	l_mod,
		real(kind(1d0)), intent(out)	tstar,
		real(kind(1d0)), intent(out)	ustar,
		real(kind(1d0)), intent(out)	zl )

Definition at line 111 of file suews_phys_atmmoiststab.f95.

 
      IMPLICIT NONE
      INTEGER, INTENT(in) :: StabilityMethod
 
      REAL(KIND(1D0)), INTENT(in) :: zzd !Active measurement height (meas. height-displac. height)
      REAL(KIND(1D0)), INTENT(in) :: z0m !Aerodynamic roughness length
      REAL(KIND(1D0)), INTENT(in) :: zdm !Displacement height
      REAL(KIND(1D0)), INTENT(in) :: avU1 !Average wind speed
      REAL(KIND(1D0)), INTENT(in) :: Temp_C !Air temperature
      REAL(KIND(1D0)), INTENT(in) :: QH_init !sensible heat flux [W m-2]
      REAL(KIND(1D0)), INTENT(in) :: avdens ! air density [kg m-3]
      REAL(KIND(1D0)), INTENT(in) :: avcp ! volumetric heat capacity [J m-3 K-1]
 
      REAL(KIND(1D0)), INTENT(out) :: L_MOD !Obukhov length
      REAL(KIND(1D0)), INTENT(out) :: TStar !T*
      REAL(KIND(1D0)), INTENT(out) :: UStar !Friction velocity
      REAL(KIND(1D0)), INTENT(out) :: zL ! Stability scale
      ! REAL(KIND(1d0)),INTENT(out)::psim   !Stability function of momentum
 
      REAL(KIND(1D0)) :: G_T_K, &
                         kuz, &
                         lold, &
                         psim, &
                         z0l, &
                         psimz0, &
                         h_init, &
                         h
      INTEGER, PARAMETER :: notUsedI = -55
 
      INTEGER :: i
 
      LOGICAL :: debug = .false.
 
      !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
      h_init = qh_init/(avdens*avcp)
 
      IF (debug) WRITE (*, *) stabilitymethod, z0m, avu1, h_init, ustar, l_mod
      g_t_k = (grav/(temp_c + 273.16))*k !gravity constant/(Temperature*Von Karman Constant)
      kuz = k*avu1 !Von Karman constant*mean wind speed
      IF (zzd < 0) CALL errorhint(32, &
                                  'Windspeed Ht too low relative to zdm [Stability calc]- values [z-zdm, zdm]', &
                                  zzd, zdm, notusedi)
 
      ustar = kuz/log(zzd/z0m) !Initial setting of u* and calc. of L_MOD (neutral situation)
      IF (abs(h_init) < 0.001) THEN ! prevent zero TStar
         h = 0.001
      ELSE
         h = h_init
      END IF
      tstar = (-h/ustar)
      l_mod = (ustar**2)/(g_t_k*tstar)
 
      IF (log(zzd/z0m) < 0.001000) THEN
         ! PRINT*, 1/(z0m-z0m)
         CALL errorhint(17, 'In stability subroutine, (z-zd) < z0.', zzd, z0m, notusedi)
      END IF
      i = 1
      lold = -999.
      z0l = z0m/l_mod !z0m roughness length
      DO WHILE ((abs(lold - l_mod) > 0.01) .AND. (i < 330)) !NT: add error threshold !Iteration starts
         lold = l_mod
         zl = zzd/l_mod
         z0l = z0m/l_mod !z0m roughness length
 
         ! IF (zL>2)THEN
         !    CALL ErrorHint(73,'LUMPS_atmos_functions_stab.f95: stability scale (z/L), UStar',zL,UStar,notUsedI)
         !    RETURN !MO-theory not necessarily valid above zL>2. Still causes problematic UStar values and correct limit might be 0.3.
         !    !Needs more investigations.
         ! END IF
 
         psim = stab_psi_mom(stabilitymethod, zl)
         psimz0 = stab_psi_mom(stabilitymethod, z0l)
 
         ustar = kuz/(log(zzd/z0m) - psim + psimz0) !Friction velocity in non-neutral situation
 
         tstar = (-h/ustar)
         l_mod = (ustar**2)/(g_t_k*tstar)
 
         ! IF(UStar<0.001000)THEN       !If u* too small
         !    UStar=KUZ/(LOG(Zzd/z0m))
         !    PRINT*, 'UStar info',UStar,KUZ,(LOG(Zzd/z0m)),Zzd,z0m
         !    CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] zl,dectime',zl,dectime,notUsedI)
         !    CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] z0l,UStar',z0l,UStar,notUsedI)
         !    CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] psim,psimz0',psim,psimz0,notUsedI)
         !    CALL ErrorHint(30,'SUBROUTINE STAB_lumps:[ u*< 0.001] AVU1,log(zzd/z0m)',AVU1,LOG(zzd/z0m),notUsedI)
         !
         !    ! RETURN
         ! ENDIF
 
         i = i + 1
      END DO
 
      ! limit zL to be within [-5,5]
      ! IF (zL < -5 .OR. zL > 5) THEN
      ! zL = MIN(5., MAX(-5., zL))
      ! limit other output variables as well as z/L
      ! L_MOD = zzd/zL
      ! limit L_mod to 3.e4 for consistency with output
      ! IF (ABS(L_MOD) > 3E4) L_MOD = L_MOD/ABS(L_MOD)*3E4
      ! z0L = z0m/L_MOD
      ! psim = stab_psi_mom(StabilityMethod, zL)
      ! psimz0 = stab_psi_mom(StabilityMethod, z0L)
 
      ! TS 01 Aug 2018:
      ! TS: limit UStar and TStar to reasonable values
      ! to prevent potential issues in other stability-related calcualtions
      ! 02 Aug 2018: set a low limit at 0.15 m/s (Schumann 1988, BLM, DOI: 10.1007/BF00123019)
      ! UStar = KUZ/(LOG(Zzd/z0m) - psim + psimz0)
      ! UStar = MAX(0.15, UStar)
      ! TStar = (-H/UStar)
 
      ! IF (UStar < 0.0001) THEN       !If u* still too small after iteration, then force quit simulation and write out error info
      !    ! UStar=KUZ/(LOG(Zzd/z0m))
      !    PRINT *, 'UStar', UStar, KUZ, (LOG(Zzd/z0m)), Zzd, z0m
      !    CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] zl,L_MOD', zl, L_MOD, notUsedI)
      !    CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] z0l,UStar', z0l, UStar, notUsedI)
      !    CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] psim,psimz0', psim, psimz0, notUsedI)
      !    CALL ErrorHint(30, 'SUBROUTINE cal_Stab:[ u*< 0.0001] AVU1,log(zzd/z0m)', AVU1, LOG(zzd/z0m), notUsedI)
 
      !    ! RETURN
      ! ENDIF
 
      ! TS 11 Feb 2021: limit UStar and TStar to reasonable ranges
      ! under all conditions, min(UStar)==0.001 m s-1 (Jimenez et al 2012, MWR, https://doi.org/10.1175/mwr-d-11-00056.1
      ustar = max(0.001, ustar)
      ! under convective/unstable condition, min(UStar)==0.15 m s-1: (Schumann 1988, BLM, https://doi.org/10.1007/BF00123019)
      IF (z0l < -neut_limit) ustar = max(0.15, ustar)
 
      ! update associated variables
      tstar = (-h/ustar)
      l_mod = (ustar**2)/(g_t_k*tstar)
      zl = zzd/l_mod
 
      ! if ( L_MOD<-990 ) then
      !   print*, 'L_MOD too low',L_MOD
      !   print*, 1/(L_MOD-L_MOD)
      !
      ! end if
 

References errorhint(), grav, k, neut_limit, and stab_psi_mom().

Referenced by suews_driver::suews_cal_resistance(), and suews_driver::suews_cal_resistance_dts().

Here is the call graph for this function:

Here is the caller graph for this function:

dpsi_dzl_g00()

real(kind(1d0)) function atmmoiststab_module::dpsi_dzl_g00	(	real(kind(1d0))	zl,
		real(kind(1d0))	psik,
		real(kind(1d0))	phik,
		real(kind(1d0))	psic,
		real(kind(1d0))	phic )

Definition at line 562 of file suews_phys_atmmoiststab.f95.

      ! Grachev et al. (2000), eqn 15
      ! note: Psi is a zL-weighted sum of Kansas and convective components
      ! derived from d(Psi)/d(zL)=d((psiK+zL**2*psiC)/(1+zL**2))/d(zL)
 
      REAL(KIND(1D0)) :: zl, dPsi, psiC, phiC, psiK, phiK
      REAL(KIND(1D0)) :: x_zl, x_psiC, x_phiC, x_psiK, x_phiK
 
      ! Kansas part
      x_psik = -(2*psik*zl)/(1 + zl**2)**2
      x_phik = -phik/(zl*(1 + zl**2))
 
      ! convective contribution
      x_psic = (2*psic*zl)/(1 + zl**2)**2
      x_phic = -(phic*zl)/(1 + zl**2)
 
      ! zL related
      x_zl = 1/zl
 
      ! totoal
      dpsi = x_psic + x_psik + x_phik + x_phic + x_zl
 

phi_cb05()

real(kind(1d0)) function atmmoiststab_module::phi_cb05	(	real(kind(1d0))	zl,
		real(kind(1d0))	k1,
		real(kind(1d0))	k2 )

Definition at line 626 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 22/24
      REAL(KIND(1D0)) :: zl, phi, k1, k2, zlk2
      zlk2 = zl**k2
      phi = 1 + k1*( &
            (zl + zlk2*(1 + zlk2)**((1 - k2)/k2)) &
            /(zl + (1 + zlk2)**(1/k2)) &
            )

Referenced by phi_heat_cb05(), and phi_mom_cb05().

Here is the caller graph for this function:

phi_conv()

real(kind(1d0)) function atmmoiststab_module::phi_conv	(	real(kind(1d0))	zl,
		real(kind(1d0))	ax )

Definition at line 548 of file suews_phys_atmmoiststab.f95.

      ! Grachev et al. (2000), eqn 15
      REAL(KIND(1D0)) :: zl, phiC, ax, dPsi_dzL
      REAL(KIND(1D0)), PARAMETER :: pi = acos(-1.)
 
      ! d(psi_conv)/d(zL)
      dpsi_dzl = (9 + 4*(-3 + sqrt(3.))*ax*zl - 9*(1 - ax*zl)**(1/3.))/ &
                 (6.*zl*(2 - 2*ax*zl + (1 - ax*zl)**(2/3.)))
 
      ! phi=1-zL*d(Psi)/d(zL)
      phic = 1 - zl*dpsi_dzl
 

Referenced by phi_heat_g00(), and phi_mom_g00().

Here is the caller graph for this function:

phi_heat_b71()

real(kind(1d0)) function atmmoiststab_module::phi_heat_b71

(

real(kind(1d0))

zl

)

Definition at line 744 of file suews_phys_atmmoiststab.f95.

      ! Businger et al. (1971)
      ! modified by Hogstrom (1988): see Table VII
      REAL(KIND(1D0)) :: zl, phih
 
      IF (abs(zl) < neut_limit) THEN
         phih = 1
      ELSEIF (zl < -neut_limit) THEN
         phih = 0.95*(1.-11.6*zl)**(-0.5)
      ELSEIF (zl > neut_limit) THEN
         phih = 0.95 + 7.8*zl
      END IF
 

References neut_limit.

Referenced by phi_heat_g00(), and stab_phi_heat().

Here is the caller graph for this function:

phi_heat_cb05()

real(kind(1d0)) function atmmoiststab_module::phi_heat_cb05

(

real(kind(1d0))

zl

)

Definition at line 650 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 24
      REAL(KIND(1D0)), PARAMETER :: c = 5.3
      REAL(KIND(1D0)), PARAMETER :: d = 1.1
      REAL(KIND(1D0)) :: zl, phih, zld
 
      IF (abs(zl) < neut_limit) THEN
         phih = 1
      ELSEIF (zl > neut_limit) THEN
         zld = zl**d
         phih = phi_cb05(zl, c, d)
      END IF
 

References neut_limit, and phi_cb05().

Referenced by phi_heat_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_heat_g00()

real(kind(1d0)) function atmmoiststab_module::phi_heat_g00

(

real(kind(1d0))

zl

)

Definition at line 512 of file suews_phys_atmmoiststab.f95.

      REAL(KIND(1D0)), PARAMETER :: ah = 34
      REAL(KIND(1D0)) :: zl, phih
      REAL(KIND(1D0)) :: phih_K ! Kansas part
      ! REAL(KIND(1D0)):: psih_K ! Kansas part
      REAL(KIND(1D0)) :: phih_C ! convective part
      ! REAL(KIND(1D0)):: psih_C ! convective part
 
      IF (abs(zl) < neut_limit) THEN
         phih = 1
      ELSEIF (zl < -neut_limit) THEN
         ! kansas
         ! psih_K = psi_heat_B71(ZL)
         phih_k = phi_heat_b71(zl)
 
         ! convective
         ! psih_C = psi_conv(ZL, ah)
         phih_c = phi_conv(zl, ah)
 
         phih = dot_product([phih_k, phih_c], [1d0, zl**2])
         phih = phih/sum([1d0, zl**2])
 
      END IF
 

References neut_limit, phi_conv(), and phi_heat_b71().

Referenced by phi_heat_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_heat_j12()

real(kind(1d0)) function atmmoiststab_module::phi_heat_j12

(

real(kind(1d0))

zl

)

Definition at line 410 of file suews_phys_atmmoiststab.f95.

      REAL(KIND(1D0)) :: zl, phih
 
      IF (abs(zl) < neut_limit) THEN
         phih = 1
      ELSEIF (zl < -neut_limit) THEN
         phih = phi_heat_g00(zl)
      ELSEIF (zl > neut_limit) THEN
         phih = phi_heat_cb05(zl)
      END IF
 

References neut_limit, phi_heat_cb05(), and phi_heat_g00().

Referenced by stab_phi_heat().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_heat_k75()

real(kind(1d0)) function atmmoiststab_module::phi_heat_k75

(

real(kind(1d0))

zl

)

Definition at line 683 of file suews_phys_atmmoiststab.f95.

      ! Kondo (1975) adopted by Campbell & Norman eqn 7.22 and 7.23 p 97
      REAL(KIND(1D0)) :: zl, phih
 
      IF (abs(zl) < neut_limit) THEN
         phih = 1
      ELSEIF (zl < -neut_limit) THEN
         phih = phi_mom_k75(zl)**2
      ELSEIF (zl > neut_limit) THEN
         phih = phi_mom_k75(zl)
      END IF
 

References neut_limit, and phi_mom_k75().

Referenced by stab_phi_heat().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_mom_b71()

real(kind(1d0)) function atmmoiststab_module::phi_mom_b71

(

real(kind(1d0))

zl

)

Definition at line 729 of file suews_phys_atmmoiststab.f95.

      ! Businger et al. (1971)
      ! modified by Hogstrom (1988): see Table VI
      REAL(KIND(1D0)) :: zl, phim
 
      IF (abs(zl) < neut_limit) THEN
         phim = 1
      ELSEIF (zl < -neut_limit) THEN
         phim = (1.-19.3*zl)**(-0.25)
      ELSEIF (zl > neut_limit) THEN
         phim = 1 + 6*zl
      END IF
 

References neut_limit.

Referenced by phi_mom_g00(), and stab_phi_mom().

Here is the caller graph for this function:

phi_mom_cb05()

real(kind(1d0)) function atmmoiststab_module::phi_mom_cb05

(

real(kind(1d0))

zl

)

Definition at line 636 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 22
      REAL(KIND(1D0)), PARAMETER :: a = 6.1
      REAL(KIND(1D0)), PARAMETER :: b = 2.5
      REAL(KIND(1D0)) :: zl, phim
 
      IF (abs(zl) < neut_limit) THEN
         phim = 1
      ELSEIF (zl > neut_limit) THEN
         phim = phi_cb05(zl, a, b)
      END IF
 

References neut_limit, and phi_cb05().

Referenced by phi_mom_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_mom_g00()

real(kind(1d0)) function atmmoiststab_module::phi_mom_g00

(

real(kind(1d0))

zl

)

Definition at line 482 of file suews_phys_atmmoiststab.f95.

      REAL(KIND(1D0)), PARAMETER :: am = 10
      REAL(KIND(1D0)) :: zl, phim
      REAL(KIND(1D0)) :: phim_K ! Kansas part
      ! REAL(KIND(1D0)):: psim_K ! Kansas part
      REAL(KIND(1D0)) :: phim_C ! convective part
      ! REAL(KIND(1D0)):: psim_C ! convective part
 
      IF (abs(zl) < neut_limit) THEN
         phim = 1
      ELSEIF (zl < -neut_limit) THEN
         ! kansas
         ! psim_K = psi_mom_B71(ZL)
         phim_k = phi_mom_b71(zl)
 
         ! convective
         ! psim_C = psi_conv(ZL, am)
         phim_c = phi_conv(zl, am)
 
         ! by definition correct but not used due to bizzare numeric behaviour
         ! phim = 1 - zL*dPsi_dzL_G00(zL, psim_K, phim_K, psim_C, phim_C)
 
         ! weighted sum:
         phim = dot_product([phim_k, phim_c], [1d0, zl**2])
         phim = phim/sum([1d0, zl**2])
 
      END IF
 

References neut_limit, phi_conv(), and phi_mom_b71().

Referenced by phi_mom_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_mom_j12()

real(kind(1d0)) function atmmoiststab_module::phi_mom_j12

(

real(kind(1d0))

zl

)

Definition at line 376 of file suews_phys_atmmoiststab.f95.

      ! not directly provided by J12
      ! equations below are derived from related original formulations
      REAL(KIND(1D0)) :: zl, phim
 
      IF (abs(zl) < neut_limit) THEN
         phim = 1
      ELSEIF (zl < -neut_limit) THEN
         ! derivative of eqn 17 in Jimenez et al. (2012)
         phim = phi_mom_g00(zl)
      ELSEIF (zl > neut_limit) THEN
         ! derivative of eqn 18 in Jimenez et al. (2012)
         phim = phi_mom_cb05(zl)
      END IF
 

References neut_limit, phi_mom_cb05(), and phi_mom_g00().

Referenced by stab_phi_mom().

Here is the call graph for this function:

Here is the caller graph for this function:

phi_mom_k75()

real(kind(1d0)) function atmmoiststab_module::phi_mom_k75

(

real(kind(1d0))

zl

)

Definition at line 669 of file suews_phys_atmmoiststab.f95.

      ! Kondo (1975) adopted by Campbell & Norman eqn 7.22 and 7.23 p 97
      REAL(KIND(1D0)) :: zl, phim
 
      IF (abs(zl) < neut_limit) THEN
         phim = 1
      ELSEIF (zl < -neut_limit) THEN
         phim = 1/(1 - 16*zl)**.25
      ELSEIF (zl > neut_limit) THEN
         phim = 1 + 6*zl/(1 + zl)
      END IF
 

References neut_limit.

Referenced by phi_heat_k75(), and stab_phi_mom().

Here is the caller graph for this function:

psi_cb05()

real(kind(1d0)) function atmmoiststab_module::psi_cb05	(	real(kind(1d0))	zl,
		real(kind(1d0))	k1,
		real(kind(1d0))	k2 )

Definition at line 592 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 21/23
      REAL(KIND(1D0)) :: zl, psi, k1, k2
      psi = -k1*log(zl + (1 + zl**k2)**(1/k2))

Referenced by psi_heat_cb05(), and psi_mom_cb05().

Here is the caller graph for this function:

psi_conv()

real(kind(1d0)) function atmmoiststab_module::psi_conv	(	real(kind(1d0))	zl,
		real(kind(1d0))	ax )

Definition at line 538 of file suews_phys_atmmoiststab.f95.

      ! Grachev et al. (2000), eqn 12
      REAL(KIND(1D0)) :: zl, psiC, y, ax
      REAL(KIND(1D0)), PARAMETER :: pi = acos(-1.)
 
      y = (1 - ax*zl)**(1/3.)
      psic = 3./2*log(y**2 + y + 1/3.) - sqrt(3.)*atan(2*y + 1/sqrt(3.)) + pi/sqrt(3.)
 

Referenced by psi_heat_g00(), and psi_mom_g00().

Here is the caller graph for this function:

psi_heat_b71()

real(kind(1d0)) function atmmoiststab_module::psi_heat_b71

(

real(kind(1d0))

zl

)

Definition at line 777 of file suews_phys_atmmoiststab.f95.

      ! integral form of phi_heat_B71
      REAL(KIND(1D0)) :: zl, psih, x
 
      IF (abs(zl) < neut_limit) THEN
         psih = 0
      ELSEIF (zl < -neut_limit) THEN
         x = 0.95*(1.-11.6*zl)**(0.5)
         psih = 2*log((1 + x)/2)
      ELSEIF (zl > neut_limit) THEN
         IF (zl < 1) THEN
            psih = (-7.8)*zl
         ELSE
            psih = (-7.8)*(1 + log(zl))
         END IF
      END IF
 

References neut_limit.

Referenced by psi_heat_g00(), and stab_psi_heat().

Here is the caller graph for this function:

psi_heat_cb05()

real(kind(1d0)) function atmmoiststab_module::psi_heat_cb05

(

real(kind(1d0))

zl

)

Definition at line 612 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 23
      REAL(KIND(1D0)), PARAMETER :: c = 5.3
      REAL(KIND(1D0)), PARAMETER :: d = 1.1
      REAL(KIND(1D0)) :: zl, psih
 
      IF (abs(zl) < neut_limit) THEN
         psih = 0
      ELSEIF (zl > neut_limit) THEN
         psih = psi_cb05(zl, c, d)
      END IF
 

References neut_limit, and psi_cb05().

Referenced by psi_heat_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_heat_g00()

real(kind(1d0)) function atmmoiststab_module::psi_heat_g00

(

real(kind(1d0))

zl

)

Definition at line 451 of file suews_phys_atmmoiststab.f95.

      ! Grachev et al. (2000), eqn 18
      ! determination of `ah` see sect. 4.1 of G00
      REAL(KIND(1D0)), PARAMETER :: ah = 34
 
      REAL(KIND(1D0)) :: zl, psih
      REAL(KIND(1D0)) :: psih_k
      REAL(KIND(1D0)) :: psih_c
 
      IF (abs(zl) < neut_limit) THEN
         psih = 0
      ELSEIF (zl < -neut_limit) THEN
         ! Kansas part:
         psih_k = psi_heat_b71(zl)
         ! convective contribution:
         psih_c = psi_conv(zl, ah)
         ! zL weighted sum:
         psih = dot_product([psih_k, psih_c], [1d0, zl**2])
         psih = psih/sum([1d0, zl**2])
 
      END IF
 

References neut_limit, psi_conv(), and psi_heat_b71().

Referenced by psi_heat_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_heat_j12()

real(kind(1d0)) function atmmoiststab_module::psi_heat_j12

(

real(kind(1d0))

zl

)

Definition at line 393 of file suews_phys_atmmoiststab.f95.

      ! Jimenez et al. (2012), eqn 17 and 19
      REAL(KIND(1D0)) :: zl, psih
 
      IF (abs(zl) < neut_limit) THEN
         psih = 0
      ELSEIF (zl < -neut_limit) THEN
         ! Jimenez et al. (2012), eqn 17
         psih = psi_heat_g00(zl)
      ELSEIF (zl > neut_limit) THEN
         ! Jimenez et al. (2012), eqn 19:
         ! Cheng and Brutsaert (2005)
         psih = psi_heat_cb05(zl)
      END IF
 

References neut_limit, psi_heat_cb05(), and psi_heat_g00().

Referenced by stab_psi_heat().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_heat_k75()

real(kind(1d0)) function atmmoiststab_module::psi_heat_k75

(

real(kind(1d0))

zl

)

Definition at line 711 of file suews_phys_atmmoiststab.f95.

      ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 and 7.27 p 97
      REAL(KIND(1D0)) :: zl, psih
 
      IF (abs(zl) < neut_limit) THEN
         psih = 0
      ELSEIF (zl < -neut_limit) THEN
         psih = (2)*log((1 + (1 - 16*zl)**0.5)/2)
      ELSEIF (zl > neut_limit) THEN
         psih = (-6)*log(1 + zl)
      END IF
 

References neut_limit.

Referenced by psi_mom_k75(), and stab_psi_heat().

Here is the caller graph for this function:

psi_mom_b71()

real(kind(1d0)) function atmmoiststab_module::psi_mom_b71

(

real(kind(1d0))

zl

)

Definition at line 759 of file suews_phys_atmmoiststab.f95.

      ! integral form of phi_mom_B71
      REAL(KIND(1D0)), PARAMETER :: PIOVER2 = acos(-1.)/2.
      REAL(KIND(1D0)) :: zl, psim, x, x2
 
      IF (abs(zl) < neut_limit) THEN
         psim = 0
      ELSEIF (zl < -neut_limit) THEN
         x = (1 - 19.3*zl)**(0.25)
         x2 = log((1 + (x**2.))/2.)
         psim = (2.*log((1 + x)/2.)) + x2 - (2.*atan(x)) + piover2
 
      ELSEIF (zl > neut_limit) THEN
         psim = (-6)*zl
      END IF
 

References neut_limit.

Referenced by psi_mom_g00(), and stab_psi_mom().

Here is the caller graph for this function:

psi_mom_cb05()

real(kind(1d0)) function atmmoiststab_module::psi_mom_cb05

(

real(kind(1d0))

zl

)

Definition at line 598 of file suews_phys_atmmoiststab.f95.

      ! Cheng and Brutsaert (2005), eqn 21
      REAL(KIND(1D0)), PARAMETER :: a = 6.1
      REAL(KIND(1D0)), PARAMETER :: b = 2.5
      REAL(KIND(1D0)) :: zl, psim
 
      IF (abs(zl) < neut_limit) THEN
         psim = 0
      ELSEIF (zl > neut_limit) THEN
         psim = psi_cb05(zl, a, b)
      END IF
 

References neut_limit, and psi_cb05().

Referenced by psi_mom_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_mom_g00()

real(kind(1d0)) function atmmoiststab_module::psi_mom_g00

(

real(kind(1d0))

zl

)

Definition at line 429 of file suews_phys_atmmoiststab.f95.

      ! Grachev et al. (2000), eqn 18
      ! determination of `am` see sect. 4.1 of G00
      REAL(KIND(1D0)), PARAMETER :: am = 10
      REAL(KIND(1D0)) :: zl, psim
      REAL(KIND(1D0)) :: psim_k
      REAL(KIND(1D0)) :: psim_c
 
      IF (abs(zl) < neut_limit) THEN
         psim = 0
      ELSEIF (zl < -neut_limit) THEN
         ! Kansas part:
         psim_k = psi_mom_b71(zl)
         ! convective contribution:
         psim_c = psi_conv(zl, am)
         ! zL weighted sum:
         psim = dot_product([psim_k, psim_c], [1d0, zl**2])
         psim = psim/sum([1d0, zl**2])
      END IF
 

References neut_limit, psi_conv(), and psi_mom_b71().

Referenced by psi_mom_j12().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_mom_j12()

real(kind(1d0)) function atmmoiststab_module::psi_mom_j12

(

real(kind(1d0))

zl

)

Definition at line 357 of file suews_phys_atmmoiststab.f95.

      ! Jimenez et al. (2012), eqn 17 and 18
      REAL(KIND(1D0)), PARAMETER :: a = 6.1
      REAL(KIND(1D0)), PARAMETER :: b = 2.5
      REAL(KIND(1D0)) :: zl, psim
 
      IF (abs(zl) < neut_limit) THEN
         psim = 0
      ELSEIF (zl < -neut_limit) THEN
         ! Jimenez et al. (2012), eqn 17
         psim = psi_mom_g00(zl)
      ELSEIF (zl > neut_limit) THEN
         ! Jimenez et al. (2012), eqn 18:
         ! Cheng and Brutsaert (2005)
         psim = psi_mom_cb05(zl)
      END IF
 

References neut_limit, psi_mom_cb05(), and psi_mom_g00().

Referenced by stab_psi_mom().

Here is the call graph for this function:

Here is the caller graph for this function:

psi_mom_k75()

real(kind(1d0)) function atmmoiststab_module::psi_mom_k75

(

real(kind(1d0))

zl

)

Definition at line 697 of file suews_phys_atmmoiststab.f95.

      ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 and 7.27 p 97
      REAL(KIND(1D0)) :: zl, psim
 
      IF (abs(zl) < neut_limit) THEN
         psim = 0
      ELSEIF (zl < -neut_limit) THEN
         psim = 0.6*psi_heat_k75(zl)
      ELSEIF (zl > neut_limit) THEN
         psim = psi_heat_k75(zl)
      END IF
 

References neut_limit, and psi_heat_k75().

Referenced by stab_psi_mom().

Here is the call graph for this function:

Here is the caller graph for this function:

stab_phi_heat()

real(kind(1d0)) function atmmoiststab_module::stab_phi_heat	(	integer	stabilitymethod,
		real(kind(1d0))	zl )

Definition at line 334 of file suews_phys_atmmoiststab.f95.

      IMPLICIT NONE
 
      REAL(KIND(1D0)) :: zl, phih
      INTEGER :: StabilityMethod
 
      SELECT CASE (stabilitymethod)
      CASE (j12)
         phih = phi_heat_j12(zl)
 
      CASE (k75)
         phih = phi_heat_k75(zl)
 
      CASE (b71)
         phih = phi_heat_b71(zl)
 
      END SELECT
 

References b71, j12, k75, phi_heat_b71(), phi_heat_j12(), and phi_heat_k75().

Referenced by rsl_module::cal_beta_lc(), rsl_module::cal_ch(), and rsl_module::cal_psih_hat().

Here is the call graph for this function:

Here is the caller graph for this function:

stab_phi_mom()

real(kind(1d0)) function atmmoiststab_module::stab_phi_mom	(	integer	stabilitymethod,
		real(kind(1d0))	zl )

Definition at line 313 of file suews_phys_atmmoiststab.f95.

      IMPLICIT NONE
 
      REAL(KIND(1D0)) :: zl, phim
      INTEGER :: StabilityMethod
 
      SELECT CASE (stabilitymethod)
      CASE (j12)
         phim = phi_mom_j12(zl)
 
      CASE (k75)
         phim = phi_mom_k75(zl)
 
      CASE (b71)
         phim = phi_mom_b71(zl)
 
      END SELECT
 

References b71, j12, k75, phi_mom_b71(), phi_mom_j12(), and phi_mom_k75().

Referenced by rsl_module::cal_cm(), and rsl_module::cal_psim_hat().

Here is the call graph for this function:

Here is the caller graph for this function:

stab_psi_heat()

real(kind(1d0)) function atmmoiststab_module::stab_psi_heat	(	integer	stabilitymethod,
		real(kind(1d0))	zl )

Definition at line 292 of file suews_phys_atmmoiststab.f95.

      IMPLICIT NONE
 
      REAL(KIND(1D0)) :: zl, psih
      INTEGER :: StabilityMethod
 
      SELECT CASE (stabilitymethod)
      CASE (j12)
         psih = psi_heat_j12(zl)
 
      CASE (k75)
         psih = psi_heat_k75(zl)
 
      CASE (b71)
         psih = psi_heat_b71(zl)
 
      END SELECT
 

References b71, j12, k75, psi_heat_b71(), psi_heat_j12(), and psi_heat_k75().

Referenced by resist_module::aerodynamicresistance(), rsl_module::rslprofile(), and rsl_module::rslprofile_dts().

Here is the call graph for this function:

Here is the caller graph for this function:

stab_psi_mom()

real(kind(1d0)) function atmmoiststab_module::stab_psi_mom	(	integer	stabilitymethod,
		real(kind(1d0))	zl )

Definition at line 271 of file suews_phys_atmmoiststab.f95.

      IMPLICIT NONE
 
      REAL(KIND(1D0)) :: zl, psim
      INTEGER :: StabilityMethod
 
      SELECT CASE (stabilitymethod)
      CASE (j12)
         psim = psi_mom_j12(zl)
 
      CASE (k75)
         psim = psi_mom_k75(zl)
 
      CASE (b71)
         psim = psi_mom_b71(zl)
 
      END SELECT
 

References b71, j12, k75, psi_mom_b71(), psi_mom_j12(), and psi_mom_k75().

Referenced by resist_module::aerodynamicresistance(), cal_stab(), rsl_module::cal_z0_rsl(), rsl_module::rslprofile(), and rsl_module::rslprofile_dts().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

b71

integer, parameter atmmoiststab_module::b71 = 4

Definition at line 10 of file suews_phys_atmmoiststab.f95.

   INTEGER, PARAMETER :: B71 = 4 ! Businger et al (1971) modifed  Hogstrom (1988)

Referenced by stab_phi_heat(), stab_phi_mom(), stab_psi_heat(), and stab_psi_mom().

grav

real(kind(1d0)), parameter atmmoiststab_module::grav = 9.80665

Definition at line 5 of file suews_phys_atmmoiststab.f95.

   REAL(KIND(1D0)), PARAMETER :: grav = 9.80665 !g - gravity - physics today august 1987

Referenced by cal_stab().

j12

integer, parameter atmmoiststab_module::j12 = 2

Definition at line 8 of file suews_phys_atmmoiststab.f95.

   INTEGER, PARAMETER :: J12 = 2 ! Cheng and Brutsaert (2005)

Referenced by stab_phi_heat(), stab_phi_mom(), stab_psi_heat(), and stab_psi_mom().

k

real(kind(1d0)), parameter atmmoiststab_module::k = 0.4

Definition at line 4 of file suews_phys_atmmoiststab.f95.

   REAL(KIND(1D0)), PARAMETER :: k = 0.4 !Von Karman's contant

Referenced by cal_stab(), and bluews_module::cbl_initial().

k75

integer, parameter atmmoiststab_module::k75 = 3

Definition at line 9 of file suews_phys_atmmoiststab.f95.

   INTEGER, PARAMETER :: k75 = 3 ! Kondo (1975) adopted by Campbell & Norman eqn 7.26 p 97

Referenced by stab_phi_heat(), stab_phi_mom(), stab_psi_heat(), and stab_psi_mom().

neut_limit

real(kind(1d0)), parameter atmmoiststab_module::neut_limit = 1.E-2

Definition at line 3 of file suews_phys_atmmoiststab.f95.

   REAL(KIND(1D0)), PARAMETER :: neut_limit = 1.e-2 !Limit for neutral stability

Referenced by cal_stab(), phi_heat_b71(), phi_heat_cb05(), phi_heat_g00(), phi_heat_j12(), phi_heat_k75(), phi_mom_b71(), phi_mom_cb05(), phi_mom_g00(), phi_mom_j12(), phi_mom_k75(), psi_heat_b71(), psi_heat_cb05(), psi_heat_g00(), psi_heat_j12(), psi_heat_k75(), psi_mom_b71(), psi_mom_cb05(), psi_mom_g00(), psi_mom_j12(), psi_mom_k75(), rsl_module::rslprofile(), and rsl_module::rslprofile_dts().