************************************************************************
*
*  Tests the "qs" function
*       qs_CC.gs
*
************************************************************************

function test_qs (args)

P  = subwrd(args,1)
T  = subwrd(args,2)
*P  = 100000
*T  = 300

say ' '
say 'P,T = 'P','T
say ' '

qsat_CC = qs_CC(P,T)

say ' '
say "From Clausius-Clapeyron's formula:"
say ' qs(kg/kg) = 'qsat_CC
say ' qs(g/kg)  = 'qsat_CC*1000.
say ' '

return

*************************************************************************

function qs_CC(P,T)

*-----------------------------------------------------------------------
*
*  Calculates the saturation specific humidity (qs, in (kg/kg))
*  given the air pressure (Pa) and temperature (K)
*
*  Use Clausius-Clapeyron formula in the definition of "qs" (using
*  Equation 2.13, page 19, in "A Short Course in Cloud Physics", 2nd
*  Edition, by R. R. Rogers):
*
*         qs = eps*es(T)/(P-(1-eps)*es(T))
*
*     where
*
*         eps   = 0.622
*         Tfr   = freezing temperature (K) (here specified as 273.15 K)
*         es(T) = 611*exp((Lv/Rv)*(T-Tfr)/(Tfr*T)) = Clausius-Clapeyron
*         P     = air pressure (Pa)
*         T     = air temperature (K)
*
*-----------------------------------------------------------------------
*
say 'Function "qs_CC":'
*
* constants required for function
*
eps = 0.622
Tfr = 273.15

say '  P(Pa) = 'P
say '  P(mb) = 'P*0.01
say '  T(K)  = 'T
say '  T(C)  = 'T-Tfr

esat = es(Tfr,T)
qs = eps*esat/(P-(1-eps)*esat)

say '  qs(kg/kg) = 'qs
say '  qs(g/kg)  = 'qs*1000.

return(qs)

*************************************************************************

function es(Tfr,T)

*-----------------------------------------------------------------------
*
*  Calculates the saturation vapor pressure ("es", in Pa) given the air
*  temperature (K)
*
*  Solve Clausius-Clapeyron formula
*
*         es = 611*exp((Lv/Rv)*(T-Tfr)/(Tfr*T))
*
*  for the saturation vapor pressure, es:
*
*     where
*
*         T   = air temperature (K)
*         Tfr = freezing temperature (K)
*         Lv  = Latent heat of vaporization  (J/kg)
*         Rv  = Gas constant for water vapor (J/(deg*kg))
*
*-----------------------------------------------------------------------
*
*say 'Function "es":'
*
* constants required for function
*
Lv = 2500000.
Rv = 461.5

*say 'T(K) = 'T
*say 'T(C) = 'T-Tfr

es = 611.*math_exp((Lv/Rv)*(T-Tfr)/(Tfr*T))

*say 'es(Pa) = 'es
*say 'es(mb) = 'es*100.

return(es)

*************************************************************************
