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

Last change on this file since 3436 was 2060, checked in by aboissinot, 6 years ago

Add the thermal plume model (cf. Rio et al. 2010) extended to gas giant.
Specific parameters are set in thermcell_mod.
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            return
158         endif
159
160         if (T.gt.T_h2O_ice_liq) then
161            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2  ! liquid / vapour
162         else if (T.lt.Tmin) then
163            print*, "careful, T<Tmin in Lcp psat water"
164            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(Tmin-Trefsublimation)**2  ! solid / vapour
165         else               
166            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2  ! solid / vapour
167         endif
168
169         dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
170         dlnpsat=RLVTT/RCPD*dummy
171         return
172      end subroutine Lcpdqsat_water
173
174
175
176
177!==================================================================
178      subroutine Tsat_water(p,Tsat)
179
180         implicit none
181
182!==================================================================
183!     Purpose
184!     -------
185!     Compute the saturation temperature
186!     for a given pressure (Pa)
187!     Based on the Tetens formula from L.Li physical parametrization manual
188!
189!     Authors
190!     -------
191!     Jeremy Leconte (2012)
192!
193!==================================================================
194
195!        input
196         real p
197 
198!        output
199         real Tsat
200
201! JL12 variables for tetens formula
202         real,parameter :: Pref_solid_liquid=611.14
203         real,parameter :: Trefvaporization=35.86
204         real,parameter :: Trefsublimation=7.66
205         real,parameter :: r3vaporization=17.269
206         real,parameter :: r3sublimation=21.875
207
208         if (p.lt.Pref_solid_liquid) then ! solid / vapour
209            Tsat =(T_h2O_ice_liq*r3sublimation- Trefsublimation*Log(p/Pref_solid_liquid))/(r3sublimation-Log(p/Pref_solid_liquid))
210         else                 ! liquid / vapour
211            Tsat =(T_h2O_ice_liq*r3vaporization- Trefvaporization*Log(p/Pref_solid_liquid))/(r3vaporization-Log(p/Pref_solid_liquid))
212         endif
213
214         return
215      end subroutine Tsat_water
216!==================================================================
217
218
219!==================================================================
220subroutine watersat(T,p,qsat)
221  implicit none
222
223!==================================================================
224!     Purpose
225!     -------
226!     Compute the water mass mixing ratio at saturation (kg/kg)
227!     for a given pressure (Pa) and temperature (K)
228!     A replacement for the old watersat.F in the Martian GCM.
229!     Based on FCTTRE.h in the LMDTERRE model.
230!
231!     JL18 watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...)
232!        which used Psat_water. This is now harmonized
233!        we put it here for archival purpose, but it is not used anymore.
234!
235!     Authors
236!     -------
237!     Robin Wordsworth (2010)
238!
239!==================================================================
240
241!   input
242  real T, p
243 
244!   output
245  real qsat
246
247! checked vs. NIST data 22/06/2010 by RW.
248! / by p gives partial pressure
249! x by epsi converts to mass mixing ratio
250
251  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
252     qsat = 100.0 * 10**(2.07023 - 0.00320991             &
253          * T - 2484.896 / T + 3.56654 * alog10(T))
254  else                 ! liquid / vapour
255     qsat = 100.0 * 10**(23.8319 - 2948.964 / T - 5.028  &
256          * alog10(T) - 29810.16 * exp( -0.0699382 * T)  &
257          + 25.21935 * exp(-2999.924/T))
258  endif
259!  qsat=epsi*qsat/p
260  if(qsat.gt.p) then
261     qsat=1.
262  else
263     qsat=epsi*qsat/(p-(1.-epsi)*qsat)
264  endif
265  return
266end subroutine watersat
267
268subroutine watersat_grad(T,qsat,dqsat)
269
270  implicit none
271
272!==================================================================
273!     Purpose
274!     -------
275!     Compute the L/cp*d (q_sat)/d T
276!     for a given temperature (K)
277!
278!     JL18: see watersat
279!
280!     Authors
281!     -------
282!     Robin Wordsworth (2010)
283!
284!==================================================================
285
286!   input
287  real T,qsat
288 
289!   output
290  real dqsat
291
292!  if (T.lt.T_coup) then ! solid / vapour !why use T_coup?????????? JL12
293  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
294     dqsat = RLVTT/RCPD*qsat*(3.56654/T             &
295          +2484.896*LOG(10.)/T**2                   &
296          -0.00320991*LOG(10.))
297  else                 ! liquid / vapour
298     dqsat = RLVTT/RCPD*qsat*LOG(10.)*              &
299          (2948.964/T**2-5.028/LOG(10.)/T           &
300          +25.21935*2999.924/T**2*EXP(-2999.924/T)  &
301          +29810.16*0.0699382*EXP(-0.0699382*T))
302  end if
303
304  return
305end subroutine watersat_grad
306
307
308
309end module watercommon_h
Note: See TracBrowser for help on using the repository browser.