SUEWS API Site
Documentation of SUEWS source code
suews_util_meteo.f95
Go to the documentation of this file.
1 !===================================================================================
2 MODULE meteo
3 
4  USE mathconstants
5  IMPLICIT NONE
6 
7  ! REAL (KIND(1d0)),PARAMETER :: PI=3.141592654
8  REAL(KIND(1d0)), PARAMETER :: rad2deg = 57.29577951
9  REAL(KIND(1d0)), PARAMETER :: deg2rad = 0.017453292
10 
11  REAL(KIND(1d0)), PARAMETER :: molmass_air = 0.028965 ! kg for 1 mol dry air
12  REAL(KIND(1d0)), PARAMETER :: molmass_co2 = 0.04401 ! kg for 1 mol CO2
13  REAL(KIND(1d0)), PARAMETER :: molmass_h2o = 0.0180153 ! kg for 1 mol water vapor
14  REAL(KIND(1d0)), PARAMETER :: mu_h2o = molmass_air/molmass_h2o ! mol air/mol H2O
15  REAL(KIND(1d0)), PARAMETER :: mu_co2 = molmass_air/molmass_co2 ! mol air/mol CO2
16  REAL(KIND(1d0)), PARAMETER :: r_dry_mol = 8.31451 ! J/K/mol gas constant
17  REAL(KIND(1D0)), PARAMETER :: r_dry_mass = r_dry_mol/molmass_air ! J/K/kg GAS CONSTANT
18  !REAL (KIND(1d0)),PARAMETER :: SIGMA_SB=5.67051e-8 ! Stefan-Boltzmann constant
19  REAL(KIND(1d0)), PARAMETER :: epsil = 0.62197
20  REAL(KIND(1d0)), PARAMETER :: kb = 1.3807e-25 ! BOLTZMANN'S CONSTANT (m^3 MB K^-1)=R/A
21  REAL(KIND(1d0)), PARAMETER :: avogadro = 6.02252e23 ! AVOGADRO'S NUMBER (molecules/mol)
22 
23 CONTAINS
24 
25  !============================================================================
26  FUNCTION sat_vap_press(TK, P) RESULT(es)
27  !c sg sept 99 f90
28  !c This uses eqns from Buck (1981) JAM 20, 1527-1532
29  !c units T (K) e (mb) P (mb)
30  !c f corrects for the fact that we are not dealing with pure water
31  REAL(KIND(1d0)) :: TK, P, TC, es, e, f
32  tc = tk - 273.15
33  IF (tc == 0) THEN
34  tc = 0.001
35  ENDIF
36  !Valid for 50>T>-40
37  e = 6.1121*exp(((18.729 - tc/227.3)*tc)/(tc + 257.87))
38  f = 1.00072 + p*(3.2e-6 + 5.9e-10*tc**2)
39  es = e*f
40  END FUNCTION sat_vap_press
41 
42  REAL(KIND(1d0)) FUNCTION sos_dryair(TK)
43  !SPEED OF SOUND IN DRY AIR, BEER (1991)
44  REAL(KIND(1d0)) ::TK
45  sos_dryair = sqrt(1.4*r_dry_mol*tk/(molmass_air*1000.))
46  END FUNCTION sos_dryair
47  !============================================================================
48  REAL(KIND(1d0)) FUNCTION potential_temp(TK, P)
49  !TK = ABSOLUTE TEMPERATURE
50  !P = PRESS (hPa)
51  REAL(KIND(1d0)) ::TK, P
52  potential_temp = tk*(1000./p)**0.286
53  END FUNCTION potential_temp
54 
55  REAL(KIND(1d0)) FUNCTION latentheat_v(TK)
56  !LATENT HEAT OF VAPORIZATION (J/kg) BOLTON(1980)
57  !TK = ABSOLUTE TEMPERATURE
58  REAL(KIND(1d0)) ::TK
59  latentheat_v = 2.501e6 - 2370.*(tk - 273.15)
60  END FUNCTION latentheat_v
61 
62  REAL(KIND(1d0)) FUNCTION latentheat_m(TK)
63  !LATENT HEAT OF MELTING (J/kg) VALID BELOW 0C BOLTON(1980)
64  !TK = ABSOLUTE TEMPERATURE
65  REAL(KIND(1d0)) ::TK, TC
66  tc = tk - 273.15
67  latentheat_m = 3.3358e5 + tc*(2030.-10.46*tc)
68  END FUNCTION latentheat_m
69 
70  REAL(KIND(1d0)) FUNCTION spec_heat_dryair(TK)
71  ! BEER (1991) APPLIED ENVIRONMETRICS METEOROLOGICAL TABLES
72  REAL(KIND(1d0)) ::TK, TC
73  tc = tk - 273.15
74  spec_heat_dryair = 1005.+((tc + 23.15)**2)/3364.
75  END FUNCTION spec_heat_dryair
76 
77  REAL(KIND(1d0)) FUNCTION spec_heat_vapor(TK, RH)
78  ! BEER (1991) APPLIED ENVIRONMETRICS METEOROLOGICAL TABLES
79  REAL(KIND(1d0)) ::TK, TC_100, RH
80  tc_100 = (tk - 273.15)/100.
81  spec_heat_vapor = 1859.+0.13*rh + (19.3 + 0.569*rh)*tc_100 + (10.+0.5*rh)*tc_100**2
82  END FUNCTION spec_heat_vapor
83 
84  REAL(KIND(1d0)) FUNCTION heatcapacity_air(TK, RH, P)
85  REAL(KIND(1d0)) ::TK, RH, P
86  REAL(KIND(1d0)) ::RHO_D, RHO_V
87  REAL(KIND(1d0)) ::CPD, CPV
88  rho_d = density_dryair(tk, p)
89  rho_v = density_vapor(tk, rh, p)
90  cpd = spec_heat_dryair(tk)
91  cpv = spec_heat_vapor(tk, rh)
92  heatcapacity_air = rho_d*cpd + rho_v*cpv
93  END FUNCTION heatcapacity_air
94 
95  REAL(KIND(1d0)) FUNCTION density_moist(TVK, P)
96  ! density of moist air FROM VIRTUAL TEMPERATURE
97  !TVK = VIRTUAL TEMPERATURE (K)
98  != = PRESSURE (hPa)
99  REAL(KIND(1d0)) ::TVK, P
100  density_moist = p*100./(r_dry_mass*tvk)
101  END FUNCTION density_moist
102 
103  REAL(KIND(1d0)) FUNCTION density_vapor(TK, RH, P)
104  !WATER VAPOR DENSITY
105  REAL(KIND(1d0)) ::TK, P, RH, EA
106  ea = sat_vap_press(tk, p)*rh/100.
107  density_vapor = (ea*100.*epsil)/(r_dry_mass*tk)
108  END FUNCTION density_vapor
109 
110  REAL(KIND(1d0)) FUNCTION density_dryair(TK, P)
111  REAL(KIND(1d0)) ::TK, P
112  density_dryair = p*100./(r_dry_mass*tk)
113  END FUNCTION density_dryair
114 
115  REAL(KIND(1d0)) FUNCTION density_gas(TK, PP, MOLMASS)
116  !DENSITY FOR IDEAL GAS SPECIES GIVEN ITS PARTIAL PRESSURE (hPa) AND MOLAR MASS (kg)
117  REAL(KIND(1d0)) ::TK, PP, MOLMASS
118  density_gas = pp*molmass/(r_dry_mol*tk)
119  END FUNCTION density_gas
120 
121  REAL(KIND(1d0)) FUNCTION partial_pressure(TK, N)
122  !PARTIAL PRESSURE OF IDEAL GAS (hPa)
123  REAL(KIND(1d0)) ::TK, N !N IS THE NUMBER DENSITY IN mol/m3
124  partial_pressure = n*kb*tk
125  END FUNCTION partial_pressure
126 
127  REAL(KIND(1d0)) FUNCTION scale_height(TK)
128  REAL(KIND(1d0)) ::TK
129  !SCALE HEIGHT FOR DRY ATMOSPHERE IN km BEER (1991)
131  END FUNCTION scale_height
132 
133  REAL(KIND(1d0)) FUNCTION vaisala_brunt_f(TK)
134  !BEER (1991)
135  REAL(KIND(1d0)) ::TK
136  vaisala_brunt_f = sqrt(0.4/1.4*9.81/scale_height(tk))
137  END FUNCTION vaisala_brunt_f
138 
139  !=====================================================================
140  ! sg sept 99 f90
141  ! Uses eqns from Buck (1981) JAM 20, 1527-1532
142  ! units T (deg C) e (mb) P (mb)
143  ! f corrects for the fact that we are not dealing with pure water
144  ! LJ Feb 2010
145  !Changed to use the updated version (Buck research manual, 1996) from Buck (1981)
146  !For water different equations in cold and warm temperatures
147 
148  FUNCTION sat_vap_press_x(Temp_c, PRESS_hPa, from, dectime) RESULT(es_hPa)
149  ! USE time
150  ! USE defaultnotUsed
151  IMPLICIT NONE
152 
153  REAL(KIND(1d0))::temp_C, press_hpa, dectime!,pw
154  REAL(KIND(1d0))::e_mb, f, press_kpa, es_hPA
155  INTEGER:: from, iv
156  INTEGER, PARAMETER::notUsedI = -55
157 
158  es_hpa = 1000 ! initialise as 1000
159 
160  !If air temperature between -0.001 -
161  IF (abs(temp_c) < 0.001000) THEN
162  IF (from == 1) THEN ! not from determining Tw
163  iv = int(press_hpa)
164  CALL errorhint(29, 'Function sat_vap_press: temp_C, dectime,press_Hpa = ', temp_c, dectime, iv)
165 
166  ENDIF
167  temp_c = 0.001000
168  ENDIF
169 
170  press_kpa = press_hpa/10
171 
172  IF (temp_c < 50 .AND. temp_c > -40) THEN
173  !e_mb=6.1121*EXP(((18.729-Temp_C/227.3)*Temp_C)/(Temp_C+257.87)) !Old one
174  !f=1.00072+Press_hPa*(3.2E-6+5.9D-10*Temp_C**2)
175 
176  IF (temp_c >= 0.001000) THEN
177  e_mb = 6.1121*exp(((18.678 - temp_c/234.5)*temp_c)/(temp_c + 257.14))
178  f = 1.00072 + press_kpa*(3.2e-6 + 5.9e-10*temp_c**2)
179  es_hpa = e_mb*f
180 
181  ELSEIF (temp_c <= -0.001000) THEN
182  e_mb = 6.1115*exp(((23.036 - temp_c/333.7)*temp_c)/(temp_c + 279.82))
183  f = 1.00022 + press_kpa*(3.83e-6 + 6.4e-10*temp_c**2)
184  es_hpa = e_mb*f
185  ENDIF
186 
187  ELSE
188  CALL errorhint(28, 'FUNCTION sat_vap_press: [Temperature is out of range], Temp_C,dectime', temp_c, dectime, notusedi)
189 
190  ENDIF
191 
192  RETURN
193  END FUNCTION sat_vap_press_x
194 
195  FUNCTION sat_vap_pressice(Temp_c, PRESS_hPa, from, dectime) RESULT(es_hPa)
196  ! USE time
197  ! USE defaultnotUsed
198  IMPLICIT NONE
199 
200  REAL(KIND(1d0))::e_mb, f, temp_C, press_hpa, press_kpa, es_hPA, dectime!,pw
201  INTEGER:: from, iv
202  INTEGER, PARAMETER::notUsedI = -55
203 
204  ! initialisation
205  es_hpa = 10
206 
207  !If air temperature between -0.001 -
208  IF (abs(temp_c) < 0.001000) THEN
209  IF (from == 1) THEN ! not from determining Tw
210  iv = int(press_hpa)
211  CALL errorhint(29, 'Function sat_vap_press: temp_C, dectime,press_Hpa = ', temp_c, dectime, iv)
212 
213  ENDIF
214  temp_c = 0.001000
215  ENDIF
216 
217  press_kpa = press_hpa/10
218 
219  IF (temp_c < 50 .AND. temp_c > -40) THEN
220  e_mb = 6.1115*exp(((23.036 - temp_c/333.7)*temp_c)/(temp_c + 279.82))
221  f = 1.00022 + press_kpa*(3.83e-6 + 6.4e-10*temp_c**2) !In hPa
222  es_hpa = e_mb*f
223 
224  ELSE
225  CALL errorhint(28, 'FUNCTION sat_vap_press: [Temperature is out of range], Temp_C,dectime', temp_c, dectime, notusedi)
226 
227  ENDIF
228 
229  RETURN
230  END FUNCTION sat_vap_pressice
231 
232  !==========================================================
233  !Output: specific humidity deficit in g/kg
234  !Input: Dry air density and air pressure in hPa
235  FUNCTION spec_hum_def(vpd_hPa, press_hPa) RESULT(dq)
236  ! USE gas
237  IMPLICIT NONE
238  REAL(KIND(1d0)) :: press_hPa, vpd_hPa, dq
239  REAL(KIND(1d0)), PARAMETER :: epsil_gkg = 621.97 !ratio molecular weight of water vapor/dry air in g/kg
240  dq = epsil_gkg*vpd_hpa/press_hpa ! Phd Thesis II.13 p 196
241  END FUNCTION spec_hum_def
242 
243  ! ==============================================================================
244  FUNCTION spec_heat_beer(Temp_C, rh, rho_v, rho_d) RESULT(cp)
245  ! Input: Air temperature, relative humidity, water vapour and dry air densities
246  ! Output: heat capacity in units J kg-1 K-1
247  ! Reference: Tom Beer, CSIRO, 1990. Applied Environmetrics Meteorological Tables.
248  ! Can be found from SG:s office from Atmmos Moist map
249  !-------------------------------------------------------------------------------
250 
251  ! USE defaultnotUsed
252  IMPLICIT NONE
253 
254  REAL(KIND(1d0))::cp, cpd, cpm, rho_v, rho_d, rh, temp_C
255 
256  !Garratt equation a20 (1992)
257  cpd = 1005.0 + ((temp_c + 23.16)**2)/3364.0 !Changed from 23.15 to 23.16
258 
259  !Beer (1990) for water vapor
260  cpm = 1859 + 0.13*rh + (19.3 + 0.569*rh)*(temp_c/100.) + &
261  (10.+0.5*rh)*(temp_c/100.)**2
262 
263  IF (abs(rho_d) < 0.000100 .OR. abs(rho_v) < 0.000100 .OR. abs(rho_d + rho_v) < 0.000100) THEN
264  CALL errorhint(42, 'spec-heat_beer', rho_v, rho_d, int(temp_c))
265  ENDIF
266 
267  cp = cpd*(rho_d/(rho_d + rho_v)) + cpm*(rho_v/(rho_d + rho_v))
268 
269  ! print*,"cp: ",cp,cpd,rho_d,rho_v,cpm,rh
270  END FUNCTION spec_heat_beer
271 
272  !==========================================================
273  !Latent_heat.f sg nov 96
274  !sg sep 99 converted f90 FUNCTION
275  !Added calcualation of latent heat of sublimation, LJ June 2012
276 
277  FUNCTION lat_vap(Temp_C, Ea_hPa, Press_hPa, cp, dectime) RESULT(lv_J_kg)
278  !Input: Air temperature, Water vapour pressure, Air pressure, heat capacity
279  !Output: latent heat of vaporization
280 
281  ! USE time
282  ! USE SnowMod
283  ! USE defaultnotUsed
284 
285  IMPLICIT NONE
286  REAL(KIND(1d0))::cp, lv_J_kg, ea_fix, tw, &
287  incr, es_tw, psyc, ea_est, press_hPa, ea_HPa, temp_C, dectime!,Temp_K
288  ! REAL(KIND(1d0))::sat_vap_press,psyc_const ! functions
289 
290  LOGICAL:: switch1 = .false., switch2 = .false.!,debug=.true.
291  INTEGER:: ii, from = 2
292  REAL(KIND(1d0)), PARAMETER::notUsed = -55.55
293 
294  ea_fix = ea_hpa
295  !if(debug) write(*,*)Temp_C, 'LV'
296  !Temp_K=temp_C+273.16
297 
298  !lv=1.91846E6*(Temp_K/(Temp_K-33.91))**2
299 
300  lv_j_kg = (2500.25 - 2.365*temp_c)*1000 !First guess for lv in units J/kg
301 
302  tw = temp_c/2. !First estimate for wet bulb temperature
303  incr = 3.
304  DO ii = 1, 100
305  IF (press_hpa < 900) THEN
306  CALL errorhint(45, 'function Lat_vap', press_hpa, notused, ii)
307  ENDIF
308 
309  ! if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'% 1',Tw
310 
311  es_tw = sat_vap_press_x(tw, press_hpa, from, dectime) !Calculate saturation vapour pressure in hPa
312 
313  !if(debug.and.dectime>55.13.and.dectime<55.2)write(35,*)'% 2',Tw
314 
315  IF (press_hpa < 900) THEN
316  CALL errorhint(45, 'function Lat_vap - 2', press_hpa, notused, ii)
317  ENDIF
318 
319  psyc = psyc_const(cp, press_hpa, lv_j_kg) !in units hPa/K
320 
321  IF (press_hpa < 900) THEN
322  CALL errorhint(45, 'function Lat _vap -31', press_hpa, notused, ii)
323  ENDIF
324 
325  ea_est = es_tw - psyc*(temp_c - tw)
326 
327  lv_j_kg = (2500.25 - 2.365*tw)*1e3
328 
329  IF (switch1 .AND. switch2) THEN
330  incr = incr/10.
331  switch1 = .false.
332  switch2 = .false.
333  ENDIF
334  IF (abs(ea_est - ea_fix) < 0.001000) THEN
335  RETURN
336  ELSEIF (ea_est > ea_fix) THEN
337  tw = tw - incr
338  switch1 = .true.
339  ELSEIF (ea_est < ea_fix) THEN
340  tw = tw + incr
341  switch2 = .true.
342  ENDIF
343  ENDDO
344 
345  RETURN
346  END FUNCTION lat_vap
347 
348  FUNCTION lat_vapsublim(Temp_C, Ea_hPa, Press_hPa, cp) RESULT(lvS_J_kg)
349  !Input: Air temperature, Water vapour pressure, Air pressure, heat capacity
350  !Output: latent heat of sublimation in units J/kg
351 
352  ! USE time
353 
354  IMPLICIT NONE
355 
356  REAL(KIND(1d0))::lvS_J_kg, temp_C, tw, incr, Ea_hPa, Press_hPa, cp
357  ! REAL(KIND(1d0))::ea_fix,es_tw,psyc,ea_est,Temp_K
358  ! REAL(KIND(1d0))::sat_vap_pressIce,psyc_const ! functions
359  ! LOGICAL:: switch1=.FALSE.,switch2=.FALSE.!,debug=.true.
360  ! INTEGER:: ii,from=2
361 
362  !Latent heat for sublimation
363  !From Rogers&Yau (A short course in cloud physics), Wikipedia
364 
365  ! ea_fix=ea_hPa
366 
367  lvs_j_kg = (2834.1 - 0.29*temp_c)*1e3 !First guess for Ls in J/kg
368 
369  tw = temp_c/2. !First estimate for wet bulb temperature
370  incr = 3.
371  press_hpa = press_hpa
372  ea_hpa = ea_hpa
373  cp = cp
374 
375  !DO ii=1,100
376 
377  ! es_tw=sat_vap_pressIce(Tw,Press_hPa,from) !Calculate saturation vapour pressure in hPa
378 
379  ! psyc=psyc_const(cp,Press_hPa,lv_J_kg)
380 
381  ! ea_est=es_tw-psyc*(temp_C-tw)
382  ! lvS_J_kg=(2834.1-0.29*tw)*1e3
383 
384  ! IF(switch1.AND.switch2)THEN
385  ! incr=incr/10.
386  ! switch1=.FALSE.
387  ! switch2=.FALSE.
388  ! ENDIF
389 
390  ! IF(ABS(ea_est-ea_fix)<0.001)THEN
391  ! RETURN
392  ! ELSEIF(ea_est > ea_fix)THEN
393  ! tw=tw-incr
394  ! switch1=.TRUE.
395  ! ELSEIF(ea_est< ea_fix)THEN
396  ! tw=tw+incr
397  ! switch2=.TRUE.
398  ! ENDIF
399  ! ENDDO
400 
401  ! RETURN
402  END FUNCTION lat_vapsublim
403 
404  !=====================================================================
405  !psyc_const.f sg nov 96
406  !sg june 99 f90
407  !calculate psyc - psychrometic constant Fritschen and Gay (1979)
408 
409  FUNCTION psyc_const(cp, Press_hPa, lv_J_kg) RESULT(psyc_hPa) !In units hPa/K
410  USE gas
411 
412  IMPLICIT NONE
413  REAL(KIND(1d0))::cp, lv_J_kg, press_hPa, psyc_hpa
414 
415  ! cp for moist air (shuttleworth p 4.13)
416  IF (cp*press_hpa < 900 .OR. lv_j_kg < 10000) THEN
417  CALL errorhint(19, &
418  'in psychrometric constant calculation: cp [J kg-1 K-1], p [hPa], Lv [J kg-1]', &
419  cp, press_hpa, int(lv_j_kg))
420  ENDIF
421 
422  psyc_hpa = (cp*press_hpa)/(epsil*lv_j_kg)
423  ! if(debug)write(*,*)psyc_hpa, 'g',cp,press_HPa,lv
424  ! LV MJKg-1
425  !www.cimis.water.ca.gov/infoEotPmEquation.jsp
426  !psyc_hPa=(0.00163*(Press_hPa/10)/LV)*10
427  ! write(*,*)psyc_hpa
428  !psyc=psyc*100.! convert into Pa
429  END FUNCTION psyc_const
430 
431  !==========================================================
432 
433  FUNCTION dewpoint(ea_hPa) RESULT(Temp_C_dew)
434  ! ea = vapor pressure (hPa)
435  ! td = dewpoint (oC)
436  !calculates dewpoint in degC from
437  ! http://www.atd.ucar.edu/weather_fl/dewpoint.html
438  ! dewpoint = (237.3 * ln(e_vp/6.1078)) / (17.27 - (ln(e_vp/6.1078)))
439 
440  REAL(KIND(1d0))::ea_hPa, temp_C_dew
441  temp_c_dew = (237.3*log(ea_hpa/6.1078))/(17.27 - (log(ea_hpa/6.1078)))
442  END FUNCTION dewpoint
443  !===============================================================================
444  FUNCTION slope_svp(temp_C) RESULT(s_hPa)
445  !COEFFICENTS FOR CALCULATING desat/dT
446  !Slope of the saturation vapor pressure vst air temperature curve
447 
448  IMPLICIT NONE
449 
450  REAL(KIND(1d0)):: b1, b2, b3, b4, b5, b6, b7, S_hPa, temp_C
451  b1 = 4.438099984d-1
452  b2 = 2.857002636d-2
453  b3 = 7.938054040d-4
454  b4 = 1.215215065d-5
455  b5 = 1.036561403d-7
456  b6 = 3.532421810d-10
457  b7 = -7.090244804d-13
458 
459  ! s - slope of saturation vapour pressure curve - Lowe (1977) -T (K)
460  ! mb /K
461  s_hpa = b1 + temp_c*(b2 + temp_c*(b3 + temp_c*(b4 + temp_c*(b5 + temp_c*(b6 + b7*temp_c)))))
462  ! write(*,*)'s',s_hpa,temp_C
463  !s_Pa=s_Pa*100 ! Pa/K
464  !www.cimis.water.ca.gov/infoEotPmEquation.jsp
465  ! s_hPa=(((4099 *(es_hPa/10))/(Temp_C+273.3)**2))*10
466  ! if(debug)write(*,*)s_hpa
467  RETURN
468  END FUNCTION slope_svp
469 
470  !===============================================================================
471  FUNCTION slopeice_svp(temp_C) RESULT(s_hPa)
472  !COEFFICENTS FOR CALCULATING desat/dT
473  !Slope of the saturation vapor pressure vst air temperature curve
474 
475  IMPLICIT NONE
476 
477  REAL(KIND(1d0)):: b1, b2, b3, b4, b5, b6, b7, S_hPa, temp_C
478 
479  b1 = 5.030305237d-1
480  b2 = 3.773255020d-2
481  b3 = 1.267995369d-3
482  b4 = 2.477563108d-5
483  b5 = 3.005693132d-7
484  b6 = 2.158542548d-9
485  b7 = 7.131097725d-12
486 
487  ! s - slope of saturation vapour pressure curve - Lowe (1977) -T (K)
488  ! mb /K
489  s_hpa = b1 + temp_c*(b2 + temp_c*(b3 + temp_c*(b4 + temp_c*(b5 + temp_c*(b6 + b7*temp_c)))))
490 
491  RETURN
492  END FUNCTION slopeice_svp
493 
494  !------------------------------------------------------------------------
495  FUNCTION qsatf(T, PMB) RESULT(qsat)
496  ! MRR, 1987
497  ! AT TEMPERATURE T (DEG C) AND PRESSURE PMB (MB), GET SATURATION SPECIFIC
498  ! HUMIDITY (KG/KG) FROM TETEN FORMULA
499 
500  REAL(KIND(1D0))::T, es, qsat, PMB
501 
502  REAL(KIND(1D0)), PARAMETER:: &
503  !Teten coefficients
504  a = 6.106, &
505  b = 17.27, &
506  c = 237.3, &
507  molar = 0.028965, & !Dry air molar fraction in kg/mol
508  molar_wat_vap = 0.0180153 !Molar fraction of water vapor in kg/mol
509 
510  IF (t > 55) THEN
511  CALL errorhint(34, 'Function qsatf', t, 0.00d0, -55)
512  ENDIF
513 
514  es = a*dexp(b*t/(c + t))
515  qsat = (molar_wat_vap/molar)*es/pmb!(rmh2o/rmair)*ES/PMB
516  END FUNCTION qsatf
517 
518  FUNCTION rh2qa(RH_dec, pres_hPa, Ta_degC) RESULT(qa_gkg)
519  ! convert relative humidity to specific humidity
520  ! TS 31 Jul 2018: initial version
521  ! Brutasert (2005) section 2.1.2, eqn 2.2, 2.4 and 2.5.
522  REAL(KIND(1D0)), INTENT(in) :: RH_dec ! relative humidity in decimal
523  REAL(KIND(1D0)), INTENT(in) :: pres_hPa ! atmospheric pressure in hPa
524  REAL(KIND(1D0)), INTENT(in) :: Ta_degC ! air temperature in degC
525 
526  REAL(KIND(1d0)) ::es ! saturation vapour pressure in hPa
527  REAL(KIND(1d0)) ::ea ! vapour pressure in hPa
528  REAL(KIND(1d0)) ::qa_gkg ! specific humidity in (g kg-1)
529 
530  es = sat_vap_press(ta_degc + 273.15, pres_hpa)
531  ea = es*rh_dec ! Brutasert (2005) section 2.1.2, eqn 2.3
532  qa_gkg = 0.622*ea/(pres_hpa - 0.378*ea)*1000 ! eqn 2.2, 2.4 and 2.5.
533 
534  END FUNCTION rh2qa
535 
536  FUNCTION qa2rh(qa_gkg, pres_hPa, Ta_degC) RESULT(RH)
537  ! convert specific humidity to relative humidity
538  ! TS 31 Jul 2018: initial version
539  ! Brutasert (2005) section 2.1.2, eqn 2.2, 2.4 and 2.5.
540  REAL(KIND(1d0)), INTENT(in) :: qa_gkg ! specific humidity in (g kg-1)
541  REAL(KIND(1D0)), INTENT(in) :: pres_hPa ! atmospheric pressure in hPa
542  REAL(KIND(1D0)), INTENT(in) :: Ta_degC ! air temperature in degC
543  REAL(KIND(1D0)) :: RH ! relative humidity in decimal
544 
545  REAL(KIND(1d0)) ::es ! saturation vapour pressure in hPa
546  REAL(KIND(1d0)) ::ea ! vapour pressure in hPa
547  REAL(KIND(1d0)) ::qa_kgkg !specific humidity in (kg kg-1)
548 
549  qa_kgkg = qa_gkg/1000
550  es = sat_vap_press(ta_degc + 273.15, pres_hpa)
551  ea = 500*pres_hpa*qa_kgkg/(311 + 189*qa_kgkg)
552  ! qa=0.622*ea/(pres_hPa-0.378*ea)*1000 ! eqn 2.2, 2.4 and 2.5.
553  rh = ea/es ! Brutasert (2005) section 2.1.2, eqn 2.3
554 
555  END FUNCTION qa2rh
556 
557 END MODULE meteo
real(kind(1d0)), parameter r_dry_mass
real(kind(1d0)) function vaisala_brunt_f(TK)
real(kind(1d0)) function potential_temp(TK, P)
real(kind(1d0)) function heatcapacity_air(TK, RH, P)
real(kind(1d0)) epsil
real(kind(1d0)), parameter mu_co2
real(kind(1d0)), parameter molmass_h2o
real(kind(1d0)) function density_dryair(TK, P)
real(kind(1d0)) function qa2rh(qa_gkg, pres_hPa, Ta_degC)
real(kind(1d0)) function spec_hum_def(vpd_hPa, press_hPa)
real(kind(1d0)) function slope_svp(temp_C)
real(kind(1d0)) function spec_heat_beer(Temp_C, rh, rho_v, rho_d)
real(kind(1d0)) function latentheat_m(TK)
real(kind(1d0)), parameter epsil
real(kind(1d0)), parameter molmass_co2
real(kind(1d0)) function partial_pressure(TK, N)
real(kind(1d0)), parameter molmass_air
real(kind(1d0)), parameter mu_h2o
real(kind(1d0)) function scale_height(TK)
real(kind(1d0)), parameter avogadro
real(kind(1d0)) function density_vapor(TK, RH, P)
real(kind(1d0)) function latentheat_v(TK)
real(kind(1d0)) function slopeice_svp(temp_C)
real(kind(1d0)) function density_moist(TVK, P)
real(kind(1d0)) function dewpoint(ea_hPa)
real(kind(1d0)) function lat_vapsublim(Temp_C, Ea_hPa, Press_hPa, cp)
real(kind(1d0)), parameter r_dry_mol
real(kind(1d0)), parameter deg2rad
real(kind(1d0)) function density_gas(TK, PP, MOLMASS)
real(kind(1d0)) function rh2qa(RH_dec, pres_hPa, Ta_degC)
real(kind(1d0)) function psyc_const(cp, Press_hPa, lv_J_kg)
real(kind(1d0)) function spec_heat_vapor(TK, RH)
real(kind(1d0)) function sat_vap_pressice(Temp_c, PRESS_hPa, from, dectime)
real(kind(1d0)) function spec_heat_dryair(TK)
real(kind(1d0)) function sat_vap_press_x(Temp_c, PRESS_hPa, from, dectime)
real(kind(1d0)) function sos_dryair(TK)
real(kind(1d0)) function lat_vap(Temp_C, Ea_hPa, Press_hPa, cp, dectime)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
real(kind(1d0)) function qsatf(T, PMB)
real(kind(1d0)), parameter rad2deg
real(kind(1d0)) function sat_vap_press(TK, P)
real(kind(1d0)), parameter kb