suews_util_meteo.f95

Go to the documentation of this file.

!===================================================================================
MODULE meteo
 
   USE mathconstants
   IMPLICIT NONE
 
   ! REAL (KIND(1d0)),PARAMETER ::  PI=3.141592654
   REAL(kind(1d0)), PARAMETER :: rad2deg = 57.29577951
   REAL(kind(1d0)), PARAMETER :: deg2rad = 0.017453292
 
   REAL(kind(1d0)), PARAMETER :: molmass_air = 0.028965 ! kg for 1 mol dry air
   REAL(kind(1d0)), PARAMETER :: molmass_co2 = 0.04401 ! kg for 1 mol CO2
   REAL(kind(1d0)), PARAMETER :: molmass_h2o = 0.0180153 ! kg for 1 mol water vapor
   REAL(kind(1d0)), PARAMETER :: mu_h2o = molmass_air/molmass_h2o ! mol air/mol H2O
   REAL(kind(1d0)), PARAMETER :: mu_co2 = molmass_air/molmass_co2 ! mol air/mol CO2
   REAL(kind(1d0)), PARAMETER :: r_dry_mol = 8.31451 ! J/K/mol gas constant
   REAL(kind(1d0)), PARAMETER :: r_dry_mass = r_dry_mol/molmass_air ! J/K/kg GAS CONSTANT
   !REAL (KIND(1d0)),PARAMETER ::  SIGMA_SB=5.67051e-8              ! Stefan-Boltzmann constant
   REAL(kind(1d0)), PARAMETER :: epsil = 0.62197
   REAL(kind(1d0)), PARAMETER :: kb = 1.3807e-25 ! BOLTZMANN'S CONSTANT (m^3 MB K^-1)=R/A
   REAL(kind(1d0)), PARAMETER :: avogadro = 6.02252e23 ! AVOGADRO'S NUMBER (molecules/mol)
 
CONTAINS
 
   !============================================================================
   FUNCTION sat_vap_press(TK, P) RESULT(es)
      !c sg sept 99 f90
      !c     This uses eqns from Buck (1981) JAM 20, 1527-1532
      !c     units T (K) e (mb) P (mb)
      !c     f corrects for the fact that we are not dealing with pure water
      REAL(kind(1d0)) :: tk, p, tc, es, e, f
      tc = tk - 273.15
      IF (tc == 0) THEN
         tc = 0.001
      END IF
      !Valid for 50>T>-40
      e = 6.1121*exp(((18.729 - tc/227.3)*tc)/(tc + 257.87))
      f = 1.00072 + p*(3.2e-6 + 5.9e-10*tc**2)
      es = e*f
   END FUNCTION sat_vap_press
 
   REAL(kind(1d0)) function sos_dryair(tk)
      !SPEED OF SOUND IN DRY AIR, BEER (1991)
      REAL(kind(1d0)) :: tk
      sos_dryair = sqrt(1.4*r_dry_mol*tk/(molmass_air*1000.))
   END FUNCTION sos_dryair
   !============================================================================
   REAL(kind(1d0)) function potential_temp(tk, p)
      !TK = ABSOLUTE TEMPERATURE
      !P  = PRESS (hPa)
      REAL(kind(1d0)) :: tk, p
      potential_temp = tk*(1000./p)**0.286
   END FUNCTION potential_temp
 
   REAL(kind(1d0)) function latentheat_v(tk)
      !LATENT HEAT OF VAPORIZATION (J/kg) BOLTON(1980)
      !TK = ABSOLUTE TEMPERATURE
      REAL(kind(1d0)) :: tk
      latentheat_v = 2.501e6 - 2370.*(tk - 273.15)
   END FUNCTION latentheat_v
 
   REAL(kind(1d0)) function latentheat_m(tk)
      !LATENT HEAT OF MELTING (J/kg) VALID BELOW 0C BOLTON(1980)
      !TK = ABSOLUTE TEMPERATURE
      REAL(kind(1d0)) :: tk, tc
      tc = tk - 273.15
      latentheat_m = 3.3358e5 + tc*(2030.-10.46*tc)
   END FUNCTION latentheat_m
 
   REAL(kind(1d0)) function spec_heat_dryair(tk)
      ! BEER (1991) APPLIED ENVIRONMETRICS METEOROLOGICAL TABLES
      REAL(kind(1d0)) :: tk, tc
      tc = tk - 273.15
      spec_heat_dryair = 1005.+((tc + 23.15)**2)/3364.
   END FUNCTION spec_heat_dryair
 
   REAL(kind(1d0)) function spec_heat_vapor(tk, rh)
      ! BEER (1991) APPLIED ENVIRONMETRICS METEOROLOGICAL TABLES
      REAL(kind(1d0)) :: tk, tc_100, rh
      tc_100 = (tk - 273.15)/100.
      spec_heat_vapor = 1859.+0.13*rh + (19.3 + 0.569*rh)*tc_100 + (10.+0.5*rh)*tc_100**2
   END FUNCTION spec_heat_vapor
 
   REAL(kind(1d0)) function heatcapacity_air(tk, rh, p)
      REAL(kind(1d0)) :: tk, rh, p
      REAL(kind(1d0)) :: rho_d, rho_v
      REAL(kind(1d0)) :: cpd, cpv
      rho_d = density_dryair(tk, p)
      rho_v = density_vapor(tk, rh, p)
      cpd = spec_heat_dryair(tk)
      cpv = spec_heat_vapor(tk, rh)
      heatcapacity_air = rho_d*cpd + rho_v*cpv
   END FUNCTION heatcapacity_air
 
   REAL(kind(1d0)) function density_moist(tvk, p)
      ! density of moist air FROM VIRTUAL TEMPERATURE
      !TVK = VIRTUAL TEMPERATURE (K)
      != = PRESSURE (hPa)
      REAL(kind(1d0)) :: tvk, p
      density_moist = p*100./(r_dry_mass*tvk)
   END FUNCTION density_moist
 
   REAL(kind(1d0)) function density_vapor(tk, rh, p)
      !WATER VAPOR DENSITY
      REAL(kind(1d0)) :: tk, p, rh, ea
      ea = sat_vap_press(tk, p)*rh/100.
      density_vapor = (ea*100.*epsil)/(r_dry_mass*tk)
   END FUNCTION density_vapor
 
   REAL(kind(1d0)) function density_dryair(tk, p)
      REAL(kind(1d0)) :: tk, p
      density_dryair = p*100./(r_dry_mass*tk)
   END FUNCTION density_dryair
 
   REAL(kind(1d0)) function density_gas(tk, pp, molmass)
      !DENSITY FOR IDEAL GAS SPECIES GIVEN ITS PARTIAL PRESSURE (hPa) AND MOLAR MASS (kg)
      REAL(kind(1d0)) :: tk, pp, molmass
      density_gas = pp*molmass/(r_dry_mol*tk)
   END FUNCTION density_gas
 
   REAL(kind(1d0)) function partial_pressure(tk, n)
      !PARTIAL PRESSURE OF IDEAL GAS (hPa)
      REAL(kind(1d0)) :: tk, n !N IS THE NUMBER DENSITY IN mol/m3
      partial_pressure = n*kb*tk
   END FUNCTION partial_pressure
 
   REAL(kind(1d0)) function scale_height(tk)
      REAL(kind(1d0)) :: tk
      !SCALE HEIGHT FOR DRY ATMOSPHERE IN km BEER (1991)
      scale_height = r_dry_mol*tk/(molmass_air*9.81)
   END FUNCTION scale_height
 
   REAL(kind(1d0)) function vaisala_brunt_f(tk)
      !BEER (1991)
      REAL(kind(1d0)) :: tk
      vaisala_brunt_f = sqrt(0.4/1.4*9.81/scale_height(tk))
   END FUNCTION vaisala_brunt_f
 
   !=====================================================================
   ! sg sept 99 f90
   ! Uses eqns from Buck (1981) JAM 20, 1527-1532
   ! units T (deg C) e (mb) P (mb)
   ! f corrects for the fact that we are not dealing with pure water
   ! LJ Feb 2010
   !Changed to use the updated version (Buck research manual, 1996) from Buck (1981)
   !For water different equations in cold and warm temperatures
 
   FUNCTION sat_vap_press_x(Temp_c, PRESS_hPa, from, dectime) RESULT(es_hPa)
      ! USE time
      ! USE defaultnotUsed
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: temp_c, press_hpa, dectime !,pw
      REAL(kind(1d0)) :: e_mb, f, press_kpa, es_hpa
      INTEGER :: from, iv
      INTEGER, PARAMETER :: notusedi = -55
 
      es_hpa = 1000 ! initialise as 1000
 
      !If air temperature between -0.001 -
      IF (abs(temp_c) < 0.001000) THEN
         IF (from == 1) THEN ! not from determining Tw
            iv = int(press_hpa)
            CALL errorhint(29, 'Function sat_vap_press: temp_C, dectime,press_Hpa = ', temp_c, dectime, iv)
 
         END IF
         temp_c = 0.001000
      END IF
 
      press_kpa = press_hpa/10
 
      IF (temp_c < 50 .AND. temp_c > -40) THEN
         !e_mb=6.1121*EXP(((18.729-Temp_C/227.3)*Temp_C)/(Temp_C+257.87)) !Old one
         !f=1.00072+Press_hPa*(3.2E-6+5.9D-10*Temp_C**2)
 
         IF (temp_c >= 0.001000) THEN
            e_mb = 6.1121*exp(((18.678 - temp_c/234.5)*temp_c)/(temp_c + 257.14))
            f = 1.00072 + press_kpa*(3.2e-6 + 5.9e-10*temp_c**2)
            es_hpa = e_mb*f
 
         ELSEIF (temp_c <= -0.001000) THEN
            e_mb = 6.1115*exp(((23.036 - temp_c/333.7)*temp_c)/(temp_c + 279.82))
            f = 1.00022 + press_kpa*(3.83e-6 + 6.4e-10*temp_c**2)
            es_hpa = e_mb*f
         END IF
 
      ELSE
         CALL errorhint(28, 'FUNCTION sat_vap_press: [Temperature is out of range], Temp_C,dectime', temp_c, dectime, notusedi)
 
      END IF
 
      RETURN
   END FUNCTION sat_vap_press_x
 
   FUNCTION sat_vap_pressice(Temp_c, PRESS_hPa, from, dectime) RESULT(es_hPa)
      ! USE time
      ! USE defaultnotUsed
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: e_mb, f, temp_c, press_hpa, press_kpa, es_hpa, dectime !,pw
      INTEGER :: from, iv
      INTEGER, PARAMETER :: notusedi = -55
 
      ! initialisation
      es_hpa = 10
 
      !If air temperature between -0.001 -
      IF (abs(temp_c) < 0.001000) THEN
         IF (from == 1) THEN ! not from determining Tw
            iv = int(press_hpa)
            CALL errorhint(29, 'Function sat_vap_press: temp_C, dectime,press_Hpa = ', temp_c, dectime, iv)
 
         END IF
         temp_c = 0.001000
      END IF
 
      press_kpa = press_hpa/10
 
      IF (temp_c < 50 .AND. temp_c > -40) THEN
         e_mb = 6.1115*exp(((23.036 - temp_c/333.7)*temp_c)/(temp_c + 279.82))
         f = 1.00022 + press_kpa*(3.83e-6 + 6.4e-10*temp_c**2) !In hPa
         es_hpa = e_mb*f
 
      ELSE
         CALL errorhint(28, 'FUNCTION sat_vap_press: [Temperature is out of range], Temp_C,dectime', temp_c, dectime, notusedi)
 
      END IF
 
      RETURN
   END FUNCTION sat_vap_pressice
 
   !==========================================================
   !Output: specific humidity deficit in g/kg
   !Input: Dry air density and air pressure in hPa
   FUNCTION spec_hum_def(vpd_hPa, press_hPa) RESULT(dq)
      ! USE gas
      IMPLICIT NONE
      REAL(kind(1d0)) :: press_hpa, vpd_hpa, dq
      REAL(kind(1d0)), PARAMETER :: epsil_gkg = 621.97 !ratio molecular weight of water vapor/dry air in g/kg
      dq = epsil_gkg*vpd_hpa/press_hpa ! Phd Thesis II.13 p 196
   END FUNCTION spec_hum_def
 
   ! ==============================================================================
   FUNCTION spec_heat_beer(Temp_C, rh, rho_v, rho_d) RESULT(cp)
      ! Input: Air temperature, relative humidity, water vapour and dry air densities
      ! Output: heat capacity in units J kg-1 K-1
      ! Reference: Tom Beer, CSIRO, 1990. Applied Environmetrics Meteorological Tables.
      ! Can be found from SG:s office from Atmmos Moist map
      !-------------------------------------------------------------------------------
 
      ! USE defaultnotUsed
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: cp, cpd, cpm, rho_v, rho_d, rh, temp_c
 
      !Garratt equation a20 (1992)
      cpd = 1005.0 + ((temp_c + 23.16)**2)/3364.0 !Changed from 23.15 to 23.16
 
      !Beer (1990) for water vapor
      cpm = 1859 + 0.13*rh + (19.3 + 0.569*rh)*(temp_c/100.) + &
            (10.+0.5*rh)*(temp_c/100.)**2
 
      IF (abs(rho_d) < 0.000100 .OR. abs(rho_v) < 0.000100 .OR. abs(rho_d + rho_v) < 0.000100) THEN
         CALL errorhint(42, 'spec-heat_beer', rho_v, rho_d, int(temp_c))
      END IF
 
      cp = cpd*(rho_d/(rho_d + rho_v)) + cpm*(rho_v/(rho_d + rho_v))
 
      !   print*,"cp: ",cp,cpd,rho_d,rho_v,cpm,rh
   END FUNCTION spec_heat_beer
 
   !==========================================================
   !Latent_heat.f sg nov 96
   !sg sep 99 converted f90 FUNCTION
   !Added calcualation of latent heat of sublimation, LJ June 2012
 
   FUNCTION lat_vap(Temp_C, Ea_hPa, Press_hPa, cp, dectime) RESULT(lv_J_kg)
      !Input: Air temperature, Water vapour pressure, Air pressure, heat capacity
      !Output: latent heat of vaporization
 
      ! USE time
      ! USE SnowMod
      ! USE defaultnotUsed
 
      IMPLICIT NONE
      REAL(kind(1d0)) :: cp, lv_j_kg, ea_fix, tw, &
                         incr, es_tw, psyc, ea_est, press_hpa, ea_hpa, temp_c, dectime !,Temp_K
      ! REAL(KIND(1d0))::sat_vap_press,psyc_const ! functions
 
      LOGICAL :: switch1 = .false., switch2 = .false. !,debug=.true.
      INTEGER :: ii, from = 2
      REAL(kind(1d0)), PARAMETER :: notused = -55.55
 
      ea_fix = ea_hpa
      !if(debug) write(*,*)Temp_C, 'LV'
      !Temp_K=temp_C+273.16
 
      !lv=1.91846E6*(Temp_K/(Temp_K-33.91))**2
 
      lv_j_kg = (2500.25 - 2.365*temp_c)*1000 !First guess for lv in units J/kg
 
      tw = temp_c/2. !First estimate for wet bulb temperature
      incr = 3.
      DO ii = 1, 100
         IF (press_hpa < 900) THEN
            CALL errorhint(45, 'function Lat_vap', press_hpa, notused, ii)
         END IF
 
         ! if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'% 1',Tw
 
         es_tw = sat_vap_press_x(tw, press_hpa, from, dectime) !Calculate saturation vapour pressure in hPa
 
         !if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'% 2',Tw
 
         IF (press_hpa < 900) THEN
            CALL errorhint(45, 'function Lat_vap - 2', press_hpa, notused, ii)
         END IF
 
         psyc = psyc_const(cp, press_hpa, lv_j_kg) !in units hPa/K
 
         IF (press_hpa < 900) THEN
            CALL errorhint(45, 'function Lat _vap -31', press_hpa, notused, ii)
         END IF
 
         ea_est = es_tw - psyc*(temp_c - tw)
 
         lv_j_kg = (2500.25 - 2.365*tw)*1e3
 
         IF (switch1 .AND. switch2) THEN
            incr = incr/10.
            switch1 = .false.
            switch2 = .false.
         END IF
         IF (abs(ea_est - ea_fix) < 0.001000) THEN
            RETURN
         ELSEIF (ea_est > ea_fix) THEN
            tw = tw - incr
            switch1 = .true.
         ELSEIF (ea_est < ea_fix) THEN
            tw = tw + incr
            switch2 = .true.
         END IF
      END DO
 
      RETURN
   END FUNCTION lat_vap
 
   FUNCTION lat_vapsublim(Temp_C, Ea_hPa, Press_hPa, cp) RESULT(lvS_J_kg)
      !Input: Air temperature, Water vapour pressure, Air pressure, heat capacity
      !Output: latent heat of sublimation in units J/kg
 
      ! USE time
 
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: lvs_j_kg, temp_c, tw, incr, ea_hpa, press_hpa, cp
      ! REAL(KIND(1d0))::ea_fix,es_tw,psyc,ea_est,Temp_K
      ! REAL(KIND(1d0))::sat_vap_pressIce,psyc_const ! functions
      ! LOGICAL:: switch1=.FALSE.,switch2=.FALSE.!,debug=.true.
      ! INTEGER:: ii,from=2
 
      !Latent heat for sublimation
      !From Rogers&Yau (A short course in cloud physics), Wikipedia
 
      ! ea_fix=ea_hPa
 
      lvs_j_kg = (2834.1 - 0.29*temp_c)*1e3 !First guess for Ls in J/kg
 
      tw = temp_c/2. !First estimate for wet bulb temperature
      incr = 3.
      press_hpa = press_hpa
      ea_hpa = ea_hpa
      cp = cp
 
      !DO ii=1,100
 
      !    es_tw=sat_vap_pressIce(Tw,Press_hPa,from)  !Calculate saturation vapour pressure in hPa
 
      ! psyc=psyc_const(cp,Press_hPa,lv_J_kg)
 
      !   ea_est=es_tw-psyc*(temp_C-tw)
      !  lvS_J_kg=(2834.1-0.29*tw)*1e3
 
      !   IF(switch1.AND.switch2)THEN
      !      incr=incr/10.
      !     switch1=.FALSE.
      !     switch2=.FALSE.
      !    ENDIF
 
      !   IF(ABS(ea_est-ea_fix)<0.001)THEN
      !     RETURN
      !   ELSEIF(ea_est > ea_fix)THEN
      !      tw=tw-incr
      !      switch1=.TRUE.
      !   ELSEIF(ea_est< ea_fix)THEN
      !      tw=tw+incr
      !      switch2=.TRUE.
      !    ENDIF
      !   ENDDO
 
      ! RETURN
   END FUNCTION lat_vapsublim
 
   !=====================================================================
   !psyc_const.f   sg   nov 96
   !sg june 99 f90
   !calculate psyc - psychrometic constant Fritschen and Gay (1979)
 
   FUNCTION psyc_const(cp, Press_hPa, lv_J_kg) RESULT(psyc_hPa) !In units hPa/K
      USE gas
 
      IMPLICIT NONE
      REAL(kind(1d0)) :: cp, lv_j_kg, press_hpa, psyc_hpa
 
      ! cp for moist air (shuttleworth p 4.13)
      IF (cp*press_hpa < 900 .OR. lv_j_kg < 10000) THEN
         CALL errorhint(19, &
                        'in psychrometric constant calculation:  cp [J kg-1 K-1], p [hPa], Lv [J kg-1]', &
                        cp, press_hpa, int(lv_j_kg))
      END IF
 
      psyc_hpa = (cp*press_hpa)/(epsil*lv_j_kg)
      !    if(debug)write(*,*)psyc_hpa, 'g',cp,press_HPa,lv
      ! LV MJKg-1
      !www.cimis.water.ca.gov/infoEotPmEquation.jsp
      !psyc_hPa=(0.00163*(Press_hPa/10)/LV)*10
      ! write(*,*)psyc_hpa
      !psyc=psyc*100.! convert into Pa
   END FUNCTION psyc_const
 
   !==========================================================
 
   FUNCTION dewpoint(ea_hPa) RESULT(Temp_C_dew)
      ! ea = vapor pressure (hPa)
      ! td = dewpoint (oC)
      !calculates dewpoint in degC from
      ! http://www.atd.ucar.edu/weather_fl/dewpoint.html
      !     dewpoint = (237.3 * ln(e_vp/6.1078)) / (17.27 - (ln(e_vp/6.1078)))
 
      REAL(kind(1d0)) :: ea_hpa, temp_c_dew
      temp_c_dew = (237.3*log(ea_hpa/6.1078))/(17.27 - (log(ea_hpa/6.1078)))
   END FUNCTION dewpoint
   !===============================================================================
   FUNCTION slope_svp(temp_C) RESULT(s_hPa)
      !COEFFICENTS FOR CALCULATING desat/dT
      !Slope of the saturation vapor pressure vst air temperature curve
 
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: b1, b2, b3, b4, b5, b6, b7, s_hpa, temp_c
      b1 = 4.438099984d-1
      b2 = 2.857002636d-2
      b3 = 7.938054040d-4
      b4 = 1.215215065d-5
      b5 = 1.036561403d-7
      b6 = 3.532421810d-10
      b7 = -7.090244804d-13
 
      !     s - slope of saturation vapour pressure curve - Lowe (1977) -T (K)
      !     mb /K
      s_hpa = b1 + temp_c*(b2 + temp_c*(b3 + temp_c*(b4 + temp_c*(b5 + temp_c*(b6 + b7*temp_c)))))
      ! write(*,*)'s',s_hpa,temp_C
      !s_Pa=s_Pa*100  ! Pa/K
      !www.cimis.water.ca.gov/infoEotPmEquation.jsp
      ! s_hPa=(((4099 *(es_hPa/10))/(Temp_C+273.3)**2))*10
      ! if(debug)write(*,*)s_hpa
      RETURN
   END FUNCTION slope_svp
 
   !===============================================================================
   FUNCTION slopeice_svp(temp_C) RESULT(s_hPa)
      !COEFFICENTS FOR CALCULATING desat/dT
      !Slope of the saturation vapor pressure vst air temperature curve
 
      IMPLICIT NONE
 
      REAL(kind(1d0)) :: b1, b2, b3, b4, b5, b6, b7, s_hpa, temp_c
 
      b1 = 5.030305237d-1
      b2 = 3.773255020d-2
      b3 = 1.267995369d-3
      b4 = 2.477563108d-5
      b5 = 3.005693132d-7
      b6 = 2.158542548d-9
      b7 = 7.131097725d-12
 
      ! s - slope of saturation vapour pressure curve - Lowe (1977) -T (K)
      ! mb /K
      s_hpa = b1 + temp_c*(b2 + temp_c*(b3 + temp_c*(b4 + temp_c*(b5 + temp_c*(b6 + b7*temp_c)))))
 
      RETURN
   END FUNCTION slopeice_svp
 
   !------------------------------------------------------------------------
   FUNCTION qsatf(T, PMB) RESULT(qsat)
      !       MRR, 1987
      ! AT TEMPERATURE T (DEG C) AND PRESSURE PMB (MB), GET SATURATION SPECIFIC
      !       HUMIDITY (KG/KG) FROM TETEN FORMULA
 
      REAL(kind(1d0)) :: t, es, qsat, pmb
 
      REAL(kind(1d0)), PARAMETER :: &
         !Teten coefficients
         a = 6.106, &
         b = 17.27, &
         c = 237.3, &
         molar = 0.028965, & !Dry air molar fraction in kg/mol
         molar_wat_vap = 0.0180153 !Molar fraction of water vapor in kg/mol
 
      IF (t > 55) THEN
         CALL errorhint(34, 'Function qsatf', t, 0.00d0, -55)
      END IF
 
      es = a*dexp(b*t/(c + t))
      qsat = (molar_wat_vap/molar)*es/pmb !(rmh2o/rmair)*ES/PMB
   END FUNCTION qsatf
 
   FUNCTION rh2qa(RH_dec, pres_hPa, Ta_degC) RESULT(qa_gkg)
      ! convert relative humidity to specific humidity
      ! TS 31 Jul 2018: initial version
      ! Brutsaert (2005) section 2.1.2, eqn 2.2, 2.4 and 2.5.
      REAL(kind(1d0)), INTENT(in) :: rh_dec ! relative humidity in decimal
      REAL(kind(1d0)), INTENT(in) :: pres_hpa ! atmospheric pressure in hPa
      REAL(kind(1d0)), INTENT(in) :: ta_degc ! air temperature in degC
 
      REAL(kind(1d0)) :: es ! saturation vapour pressure in hPa
      REAL(kind(1d0)) :: ea ! vapour pressure in hPa
      REAL(kind(1d0)) :: qa_gkg ! specific humidity in (g kg-1)
 
      es = sat_vap_press(ta_degc + 273.15, pres_hpa)
      ea = es*rh_dec ! Brutsaert (2005) section 2.1.2, eqn 2.3
      qa_gkg = 0.622*ea/(pres_hpa - 0.378*ea)*1000 ! eqn 2.2, 2.4 and 2.5.
 
   END FUNCTION rh2qa
 
   FUNCTION qa2rh(qa_gkg, pres_hPa, Ta_degC) RESULT(RH)
      ! convert specific humidity to relative humidity
      ! TS 31 Jul 2018: initial version
      ! Brutsaert (2005) section 2.1.2, eqn 2.2, 2.4 and 2.5.
      REAL(kind(1d0)), INTENT(in) :: qa_gkg ! specific humidity in (g kg-1)
      REAL(kind(1d0)), INTENT(in) :: pres_hpa ! atmospheric pressure in hPa
      REAL(kind(1d0)), INTENT(in) :: ta_degc ! air temperature in degC
      REAL(kind(1d0)) :: rh ! relative humidity in decimal
 
      REAL(kind(1d0)) :: es ! saturation vapour pressure in hPa
      REAL(kind(1d0)) :: ea ! vapour pressure in hPa
      REAL(kind(1d0)) :: qa_kgkg !specific humidity in (kg kg-1)
 
      qa_kgkg = qa_gkg/1000
      es = sat_vap_press(ta_degc + 273.15, pres_hpa)
      ea = 500*pres_hpa*qa_kgkg/(311 + 189*qa_kgkg)
      ! qa=0.622*ea/(pres_hPa-0.378*ea)*1000 ! eqn 2.2, 2.4 and 2.5.
      rh = ea/es ! Brutsaert (2005) section 2.1.2, eqn 2.3
 
      ! set an upper limit
      IF (rh > 1) rh = 1
 
   END FUNCTION qa2rh
 
END MODULE meteo