source: trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90 @ 3523

Last change on this file since 3523 was 3489, checked in by mmaurice, 4 weeks ago

Generic PCM:

Bug fix: initialize dlnpsat to zero in function Lcpdqsat_water even if
p<psat to avoid moistadj to crash in 1D in debug mode.

MM
File size: 9.3 KB
Line 
1module watercommon_h
2
3      implicit none
4
5      real, parameter :: T_coup = 234.0
6      real, parameter :: T_h2O_ice_liq = 273.16
7      real, parameter :: T_h2O_ice_clouds = T_h2O_ice_liq-15.
8      real, parameter :: mH2O = 18.01528   
9
10      ! benjamin additions
11      real, parameter :: RLVTT = 2.257E+6 ! Latent heat of vaporization (J kg-1)
12      real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1)
13
14      real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) ! entails an energy sink but better description of albedo
15      real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3)
16      real, parameter :: rhowaterice = 9.2E+2 ! mass of water (kg/m^3)
17      real, parameter :: capcal_h2o_liq = 4181.3 ! specific heat capacity of liquid water J/kg/K
18      real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2)
19
20      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
21      real, save :: RETV
22      real, save :: RLvCp
23!$OMP THREADPRIVATE(epsi,RCPD,RCPV,RV,RVTMP2)
24     
25      contains
26
27!==================================================================
28      subroutine su_watercycle
29
30      use comcstfi_mod, only: r, cpp, mugaz
31      implicit none
32
33
34!==================================================================
35!
36!     Purpose
37!     -------
38!     Set up relevant constants and parameters for the water cycle, and water cloud properties
39!
40!     Authors
41!     -------
42!     Robin Wordsworth (2010)
43!     Jeremy Leconte (2012)
44!
45!==================================================================
46
47      epsi   = mH2O / mugaz
48      RCPD   = cpp
49
50      !RV = 1000.*R/mH2O
51      RV = 1000.*8.314/mH2O ! caution! R is R/mugaz already!
52
53      RCPV   = 1.88e3 ! specific heat capacity of water vapor at 350K
54
55      RVTMP2 = RCPV/RCPD-1. ! not currently used...
56     
57! AB : initializations added for the thermal plume model
58      RETV = RV / r - 1.
59      RLvCp = RLVTT / RCPD
60     
61      end subroutine su_watercycle
62     
63
64
65!==================================================================
66      subroutine Psat_water(T,p,psat,qsat)
67
68         implicit none
69
70!==================================================================
71!     Purpose
72!     -------
73!     Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg)
74!     for a given pressure (Pa) and temperature (K)
75!     Based on the Tetens formula from L.Li physical parametrization manual
76!
77!     Authors
78!     -------
79!     Jeremy Leconte (2012)
80!
81!==================================================================
82
83!        input
84         real, intent(in) :: T, p
85 
86!        output
87         real psat,qsat
88
89! JL12 variables for tetens formula
90         real,parameter :: Pref_solid_liquid=611.14
91         real,parameter :: Trefvaporization=35.86
92         real,parameter :: Trefsublimation=7.66
93         real,parameter :: Tmin=8.
94         real,parameter :: r3vaporization=17.269
95         real,parameter :: r3sublimation=21.875
96
97! checked vs. old watersat data 14/05/2012 by JL.
98
99         if (T.gt.T_h2O_ice_liq) then
100            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization)) ! liquid / vapour
101         else if (T.lt.Tmin) then
102            print*, "careful, T<Tmin in psat water"
103          !  psat = Pref_solid_liquid*Exp(r3sublimation*(Tmin-T_h2O_ice_liq)/(Tmin-Trefsublimation)) ! min psat 
104         ! Ehouarn: gfortran says: Error: Result of EXP underflows its kind,
105         !          so set psat to the smallest possible value instead
106            psat=tiny(psat)
107         else                 
108            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) ! solid / vapour
109         endif
110         if(psat.gt.p) then
111            qsat=1.
112         else
113            qsat=epsi*psat/(p-(1.-epsi)*psat)
114         endif
115         return
116      end subroutine Psat_water
117
118
119
120
121!==================================================================
122      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat,dlnpsat)
123
124         implicit none
125
126!==================================================================
127!     Purpose
128!     -------
129!     Compute dqsat=L/cp*d (q_sat)/d T and dlnpsat=L/cp d(ln Psat)/d T
130!     for a given temperature (K)!
131!     Based on the Tetens formula from L.Li physical parametrization manual
132!
133!     Authors
134!     -------
135!     Jeremy Leconte (2012)
136!
137!==================================================================
138
139!        input
140         real T, p, psat, qsat
141 
142!        output
143         real dqsat,dlnpsat
144
145! JL12 variables for tetens formula
146         real,parameter :: Pref_solid_liquid=611.14
147         real,parameter :: Trefvaporization=35.86
148         real,parameter :: Tmin=8.
149         real,parameter :: Trefsublimation=7.66
150         real,parameter :: r3vaporization=17.269
151         real,parameter :: r3sublimation=21.875
152
153         real :: dummy
154
155         if (psat.gt.p) then
156            dqsat=0.
157       dlnpsat=0.
158            return
159         endif
160
161         if (T.gt.T_h2O_ice_liq) then
162            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2  ! liquid / vapour
163         else if (T.lt.Tmin) then
164            print*, "careful, T<Tmin in Lcp psat water"
165            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(Tmin-Trefsublimation)**2  ! solid / vapour
166         else               
167            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2  ! solid / vapour
168         endif
169
170         dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
171         dlnpsat=RLVTT/RCPD*dummy
172         return
173      end subroutine Lcpdqsat_water
174
175
176
177
178!==================================================================
179      subroutine Tsat_water(p,Tsat)
180
181         implicit none
182
183!==================================================================
184!     Purpose
185!     -------
186!     Compute the saturation temperature
187!     for a given pressure (Pa)
188!     Based on the Tetens formula from L.Li physical parametrization manual
189!
190!     Authors
191!     -------
192!     Jeremy Leconte (2012)
193!
194!==================================================================
195
196!        input
197         real p
198 
199!        output
200         real Tsat
201
202! JL12 variables for tetens formula
203         real,parameter :: Pref_solid_liquid=611.14
204         real,parameter :: Trefvaporization=35.86
205         real,parameter :: Trefsublimation=7.66
206         real,parameter :: r3vaporization=17.269
207         real,parameter :: r3sublimation=21.875
208
209         if (p.lt.Pref_solid_liquid) then ! solid / vapour
210            Tsat =(T_h2O_ice_liq*r3sublimation- Trefsublimation*Log(p/Pref_solid_liquid))/(r3sublimation-Log(p/Pref_solid_liquid))
211         else                 ! liquid / vapour
212            Tsat =(T_h2O_ice_liq*r3vaporization- Trefvaporization*Log(p/Pref_solid_liquid))/(r3vaporization-Log(p/Pref_solid_liquid))
213         endif
214
215         return
216      end subroutine Tsat_water
217!==================================================================
218
219
220!==================================================================
221subroutine watersat(T,p,qsat)
222  implicit none
223
224!==================================================================
225!     Purpose
226!     -------
227!     Compute the water mass mixing ratio at saturation (kg/kg)
228!     for a given pressure (Pa) and temperature (K)
229!     A replacement for the old watersat.F in the Martian GCM.
230!     Based on FCTTRE.h in the LMDTERRE model.
231!
232!     JL18 watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...)
233!        which used Psat_water. This is now harmonized
234!        we put it here for archival purpose, but it is not used anymore.
235!
236!     Authors
237!     -------
238!     Robin Wordsworth (2010)
239!
240!==================================================================
241
242!   input
243  real T, p
244 
245!   output
246  real qsat
247
248! checked vs. NIST data 22/06/2010 by RW.
249! / by p gives partial pressure
250! x by epsi converts to mass mixing ratio
251
252  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
253     qsat = 100.0 * 10**(2.07023 - 0.00320991             &
254          * T - 2484.896 / T + 3.56654 * alog10(T))
255  else                 ! liquid / vapour
256     qsat = 100.0 * 10**(23.8319 - 2948.964 / T - 5.028  &
257          * alog10(T) - 29810.16 * exp( -0.0699382 * T)  &
258          + 25.21935 * exp(-2999.924/T))
259  endif
260!  qsat=epsi*qsat/p
261  if(qsat.gt.p) then
262     qsat=1.
263  else
264     qsat=epsi*qsat/(p-(1.-epsi)*qsat)
265  endif
266  return
267end subroutine watersat
268
269subroutine watersat_grad(T,qsat,dqsat)
270
271  implicit none
272
273!==================================================================
274!     Purpose
275!     -------
276!     Compute the L/cp*d (q_sat)/d T
277!     for a given temperature (K)
278!
279!     JL18: see watersat
280!
281!     Authors
282!     -------
283!     Robin Wordsworth (2010)
284!
285!==================================================================
286
287!   input
288  real T,qsat
289 
290!   output
291  real dqsat
292
293!  if (T.lt.T_coup) then ! solid / vapour !why use T_coup?????????? JL12
294  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
295     dqsat = RLVTT/RCPD*qsat*(3.56654/T             &
296          +2484.896*LOG(10.)/T**2                   &
297          -0.00320991*LOG(10.))
298  else                 ! liquid / vapour
299     dqsat = RLVTT/RCPD*qsat*LOG(10.)*              &
300          (2948.964/T**2-5.028/LOG(10.)/T           &
301          +25.21935*2999.924/T**2*EXP(-2999.924/T)  &
302          +29810.16*0.0699382*EXP(-0.0699382*T))
303  end if
304
305  return
306end subroutine watersat_grad
307
308
309
310end module watercommon_h
Note: See TracBrowser for help on using the repository browser.