source: flexpart.git/src/qvsat.f90 @ 02095e3

FPv9.3.1FPv9.3.1b_testingFPv9.3.2NetCDFdepositiondevflexpart-noresmfp9.3.1-20161214-nc4grib2nc4_repairinputlistlaptoprelease-10svn-petrasvn-trunkunivie
Last change on this file since 02095e3 was e200b7a, checked in by Matthias Langer <matthias.langer@…>, 7 years ago

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/trunk@4 ef8cc7e1-21b7-489e-abab-c1baa636049d

  • Property mode set to 100644
File size: 5.3 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22!##################################################################
23!##################################################################
24!######                                                      ######
25!######                     Developed by                     ######
26!######     Center for Analysis and Prediction of Storms     ######
27!######                University of Oklahoma                ######
28!######                                                      ######
29!##################################################################
30!##################################################################
31
32function f_qvsat( p, t )
33
34  !PURPOSE:
35  !
36  !Calculate the saturation specific humidity using enhanced Teten's
37  !formula.
38  !
39  !AUTHOR: Yuhe Liu
40  !01/08/1998
41  !
42  !MODIFICATION HISTORY:
43  !
44  !INPUT :
45  !  p        Pressure (Pascal)
46  !  t        Temperature (K)
47  !OUTPUT:
48  !  f_qvsat  Saturation water vapor specific humidity (kg/kg).
49  !
50  !Variable Declarations.
51  !
52
53  implicit none
54
55  real :: p         ! Pressure (Pascal)
56  real :: t         ! Temperature (K)
57  real :: f_qvsat   ! Saturation water vapor specific humidity (kg/kg)
58  real :: f_esl,f_esi,fespt
59
60  real,parameter ::  rd     = 287.0 ! Gas constant for dry air  (m**2/(s**2*K))
61  real,parameter ::  rv     = 461.0 ! Gas constant for water vapor  (m**2/(s**2*K)).
62  real,parameter ::  rddrv  = rd/rv
63
64
65  ! Change by A. Stohl to save computation time:
66  ! IF ( t.ge.273.15 ) THEN     ! for water
67  if ( t.ge.253.15 ) then      ! modification Petra Seibert
68                               ! (supercooled water may be present)
69    fespt=f_esl(p,t)
70  else
71    fespt=f_esi(p,t)
72  endif
73
74!!$  f_qvsat = rddrv * fespt / (p-(1.0-rddrv)*fespt)      !old
75  if (p-(1.0-rddrv)*fespt == 0.) then                     !bugfix
76     f_qvsat = 1.
77  else
78     f_qvsat = rddrv * fespt / (p-(1.0-rddrv)*fespt)
79  end if
80
81  return
82end function f_qvsat
83
84
85function f_esl( p, t )
86
87  implicit none
88
89  real :: p         ! Pressure (Pascal)
90  real :: t         ! Temperature (K)
91  real :: f_esl     ! Saturation water vapor pressure over liquid water
92
93  real :: f
94
95  !#######################################################################
96  !
97  !Saturation specific humidity parameters used in enhanced Teten's
98  !formula. (See A. Buck, JAM 1981)
99  !
100  !#######################################################################
101
102  real,parameter ::  satfwa = 1.0007
103  real,parameter ::  satfwb = 3.46e-8  ! for p in Pa
104
105  real,parameter ::  satewa = 611.21   ! es in Pa
106  real,parameter ::  satewb = 17.502
107  real,parameter ::  satewc = 32.18
108
109  real,parameter ::  satfia = 1.0003
110  real,parameter ::  satfib = 4.18e-8  ! for p in Pa
111
112  real,parameter ::  sateia = 611.15   ! es in Pa
113  real,parameter ::  sateib = 22.452
114  real,parameter ::  sateic = 0.6
115
116  f = satfwa + satfwb * p
117  f_esl = f * satewa * exp( satewb*(t-273.15)/(t-satewc) )
118
119  return
120end function f_esl
121
122function f_esi( p, t )
123
124  implicit none
125
126  real :: p         ! Pressure (Pascal)
127  real :: t         ! Temperature (K)
128  real :: f_esi     ! Saturation water vapor pressure over ice (Pa)
129
130  real :: f
131
132  !#######################################################################
133  !
134  !Saturation specific humidity parameters used in enhanced Teten's
135  !formula. (See A. Buck, JAM 1981)
136  !
137  !#######################################################################
138  !
139  real,parameter ::  satfwa = 1.0007
140  real,parameter ::  satfwb = 3.46e-8  ! for p in Pa
141
142  real,parameter ::  satewa = 611.21   ! es in Pa
143  real,parameter ::  satewb = 17.502
144  real,parameter ::  satewc = 32.18
145
146  real,parameter ::  satfia = 1.0003
147  real,parameter ::  satfib = 4.18e-8  ! for p in Pa
148
149  real,parameter ::  sateia = 611.15   ! es in Pa
150  real,parameter ::  sateib = 22.452
151  real,parameter ::  sateic = 0.6
152
153  f = satfia + satfib * p
154  f_esi = f * sateia * exp( sateib*(t-273.15)/(t-sateic) )
155
156  return
157end function f_esi
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG