source: flexpart.git/src/qvsat.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 4.0 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4!##################################################################
5!##################################################################
6!######                                                      ######
7!######                     Developed by                     ######
8!######     Center for Analysis and Prediction of Storms     ######
9!######                University of Oklahoma                ######
10!######                                                      ######
11!##################################################################
12!##################################################################
13
14function f_qvsat( p, t )
15
16  !PURPOSE:
17  !
18  !Calculate the saturation specific humidity using enhanced Teten's
19  !formula.
20  !
21  !AUTHOR: Yuhe Liu
22  !01/08/1998
23  !
24  !MODIFICATION HISTORY:
25  !
26  !INPUT :
27  !  p        Pressure (Pascal)
28  !  t        Temperature (K)
29  !OUTPUT:
30  !  f_qvsat  Saturation water vapor specific humidity (kg/kg).
31  !
32  !Variable Declarations.
33  !
34
35  implicit none
36
37  real :: p         ! Pressure (Pascal)
38  real :: t         ! Temperature (K)
39  real :: f_qvsat   ! Saturation water vapor specific humidity (kg/kg)
40  real :: f_esl,f_esi,fespt
41
42  real,parameter ::  rd     = 287.0 ! Gas constant for dry air  (m**2/(s**2*K))
43  real,parameter ::  rv     = 461.0 ! Gas constant for water vapor  (m**2/(s**2*K)).
44  real,parameter ::  rddrv  = rd/rv
45
46
47  ! Change by A. Stohl to save computation time:
48  ! IF ( t.ge.273.15 ) THEN     ! for water
49  if ( t.ge.253.15 ) then      ! modification Petra Seibert
50                               ! (supercooled water may be present)
51    fespt=f_esl(p,t)
52  else
53    fespt=f_esi(p,t)
54  endif
55
56!!$  f_qvsat = rddrv * fespt / (p-(1.0-rddrv)*fespt)      !old
57  if (p-(1.0-rddrv)*fespt == 0.) then                     !bugfix
58     f_qvsat = 1.
59  else
60     f_qvsat = rddrv * fespt / (p-(1.0-rddrv)*fespt)
61  end if
62
63  return
64end function f_qvsat
65
66
67function f_esl( p, t )
68
69  implicit none
70
71  real :: p         ! Pressure (Pascal)
72  real :: t         ! Temperature (K)
73  real :: f_esl     ! Saturation water vapor pressure over liquid water
74
75  real :: f
76
77  !#######################################################################
78  !
79  !Saturation specific humidity parameters used in enhanced Teten's
80  !formula. (See A. Buck, JAM 1981)
81  !
82  !#######################################################################
83
84  real,parameter ::  satfwa = 1.0007
85  real,parameter ::  satfwb = 3.46e-8  ! for p in Pa
86
87  real,parameter ::  satewa = 611.21   ! es in Pa
88  real,parameter ::  satewb = 17.502
89  real,parameter ::  satewc = 32.18
90
91  real,parameter ::  satfia = 1.0003
92  real,parameter ::  satfib = 4.18e-8  ! for p in Pa
93
94  real,parameter ::  sateia = 611.15   ! es in Pa
95  real,parameter ::  sateib = 22.452
96  real,parameter ::  sateic = 0.6
97
98  f = satfwa + satfwb * p
99  f_esl = f * satewa * exp( satewb*(t-273.15)/(t-satewc) )
100
101  return
102end function f_esl
103
104function f_esi( p, t )
105
106  implicit none
107
108  real :: p         ! Pressure (Pascal)
109  real :: t         ! Temperature (K)
110  real :: f_esi     ! Saturation water vapor pressure over ice (Pa)
111
112  real :: f
113
114  !#######################################################################
115  !
116  !Saturation specific humidity parameters used in enhanced Teten's
117  !formula. (See A. Buck, JAM 1981)
118  !
119  !#######################################################################
120  !
121  real,parameter ::  satfwa = 1.0007
122  real,parameter ::  satfwb = 3.46e-8  ! for p in Pa
123
124  real,parameter ::  satewa = 611.21   ! es in Pa
125  real,parameter ::  satewb = 17.502
126  real,parameter ::  satewc = 32.18
127
128  real,parameter ::  satfia = 1.0003
129  real,parameter ::  satfib = 4.18e-8  ! for p in Pa
130
131  real,parameter ::  sateia = 611.15   ! es in Pa
132  real,parameter ::  sateib = 22.452
133  real,parameter ::  sateic = 0.6
134
135  f = satfia + satfib * p
136  f_esi = f * sateia * exp( sateib*(t-273.15)/(t-sateic) )
137
138  return
139end function f_esi
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG