Line data Source code
1 : !-------------------------------------------------------------------------
2 : !$Id$
3 : !===============================================================================
4 : module saturation
5 :
6 : ! Description:
7 : ! Contains functions that compute saturation with respect
8 : ! to liquid or ice.
9 : !-----------------------------------------------------------------------
10 :
11 : use model_flags, only: & ! h1g, 2010-06-18
12 : I_sat_sphum
13 :
14 : use clubb_precision, only: &
15 : core_rknd ! Variable(s)
16 :
17 : use model_flags, only: &
18 : saturation_bolton, &
19 : saturation_gfdl, &
20 : saturation_flatau, &
21 : saturation_lookup
22 :
23 : implicit none
24 :
25 : private ! Change default so all items private
26 :
27 : public :: sat_mixrat_liq, sat_mixrat_ice, rcm_sat_adj, sat_vapor_press_liq
28 :
29 : private :: sat_vapor_press_ice_flatau, sat_vapor_press_ice_bolton
30 :
31 : interface sat_mixrat_liq
32 : module procedure sat_mixrat_liq_k ! Works over a single vertical level
33 : module procedure sat_mixrat_liq_2D ! Works over all vertical levels and columns
34 : end interface sat_mixrat_liq
35 :
36 : interface sat_mixrat_ice
37 : module procedure sat_mixrat_ice_k ! Works over a single vertical level
38 : module procedure sat_mixrat_ice_2D ! Works over all vertical levels and columns
39 : end interface sat_mixrat_ice
40 :
41 : ! Lookup table of values for saturation
42 : real( kind = core_rknd ), parameter, private, dimension(188:343) :: &
43 : svp_liq_lookup_table = &
44 : (/ 0.049560547_core_rknd, 0.059753418_core_rknd, 0.070129395_core_rknd, 0.083618164_core_rknd, &
45 : 0.09814453_core_rknd, 0.11444092_core_rknd, 0.13446045_core_rknd, 0.15686035_core_rknd, &
46 : 0.18218994_core_rknd, 0.21240234_core_rknd, 0.24725342_core_rknd, 0.28668213_core_rknd, &
47 : 0.33184814_core_rknd, 0.3826294_core_rknd, 0.4416504_core_rknd, 0.50775146_core_rknd, &
48 : 0.58343506_core_rknd, 0.6694946_core_rknd, 0.7668457_core_rknd, 0.87750244_core_rknd, &
49 : 1.0023804_core_rknd, 1.1434937_core_rknd, 1.3028564_core_rknd, 1.482544_core_rknd, &
50 : 1.6847534_core_rknd, 1.9118042_core_rknd, 2.1671143_core_rknd, 2.4535522_core_rknd, &
51 : 2.774231_core_rknd, 3.1330566_core_rknd, 3.5343628_core_rknd, 3.9819336_core_rknd, &
52 : 4.480713_core_rknd, 5.036072_core_rknd, 5.6540527_core_rknd, 6.340088_core_rknd, &
53 : 7.1015015_core_rknd, 7.9450684_core_rknd, 8.8793335_core_rknd, 9.91217_core_rknd, &
54 : 11.053528_core_rknd, 12.313049_core_rknd, 13.70166_core_rknd, 15.231018_core_rknd, &
55 : 16.91394_core_rknd, 18.764038_core_rknd, 20.795898_core_rknd, 23.025574_core_rknd, &
56 : 25.470093_core_rknd, 28.147766_core_rknd, 31.078003_core_rknd, 34.282043_core_rknd, &
57 : 37.782593_core_rknd, 41.60382_core_rknd, 45.771606_core_rknd, 50.31366_core_rknd, &
58 : 55.259644_core_rknd, 60.641174_core_rknd, 66.492004_core_rknd, 72.84802_core_rknd, &
59 : 79.74756_core_rknd, 87.23126_core_rknd, 95.34259_core_rknd, 104.12747_core_rknd, &
60 : 113.634796_core_rknd, 123.91641_core_rknd, 135.02725_core_rknd, 147.02563_core_rknd, &
61 : 159.97308_core_rknd, 173.93488_core_rknd, 188.97995_core_rknd, 205.18109_core_rknd, &
62 : 222.61517_core_rknd, 241.36334_core_rknd, 261.51108_core_rknd, 283.14853_core_rknd, &
63 : 306.37054_core_rknd, 331.27698_core_rknd, 357.97278_core_rknd, 386.56842_core_rknd, &
64 : 417.17978_core_rknd, 449.9286_core_rknd, 484.94254_core_rknd, 522.3556_core_rknd, &
65 : 562.30804_core_rknd, 604.947_core_rknd, 650.42645_core_rknd, 698.9074_core_rknd, &
66 : 750.55835_core_rknd, 805.55554_core_rknd, 864.0828_core_rknd, 926.3325_core_rknd, &
67 : 992.5052_core_rknd, 1062.8102_core_rknd, 1137.4657_core_rknd, 1216.6995_core_rknd, &
68 : 1300.7483_core_rknd, 1389.8594_core_rknd, 1484.2896_core_rknd, 1584.3064_core_rknd, &
69 : 1690.1881_core_rknd, 1802.224_core_rknd, 1920.7146_core_rknd, 2045.9724_core_rknd, &
70 : 2178.3218_core_rknd, 2318.099_core_rknd, 2465.654_core_rknd, 2621.3489_core_rknd, &
71 : 2785.5596_core_rknd, 2958.6758_core_rknd, 3141.101_core_rknd, 3333.2534_core_rknd, &
72 : 3535.5657_core_rknd, 3748.4863_core_rknd, 3972.4792_core_rknd, 4208.024_core_rknd, &
73 : 4455.616_core_rknd, 4715.7686_core_rknd, 4989.0127_core_rknd, 5275.8945_core_rknd, &
74 : 5576.9795_core_rknd, 5892.8535_core_rknd, 6224.116_core_rknd, 6571.3926_core_rknd, &
75 : 6935.3213_core_rknd, 7316.5674_core_rknd, 7715.8105_core_rknd, 8133.755_core_rknd, &
76 : 8571.125_core_rknd, 9028.667_core_rknd, 9507.15_core_rknd, 10007.367_core_rknd, &
77 : 10530.132_core_rknd, 11076.282_core_rknd, 11646.683_core_rknd, 12242.221_core_rknd, &
78 : 12863.808_core_rknd, 13512.384_core_rknd, 14188.913_core_rknd, 14894.385_core_rknd, &
79 : 15629.823_core_rknd, 16396.268_core_rknd, 17194.799_core_rknd, 18026.516_core_rknd, &
80 : 18892.55_core_rknd, 19794.07_core_rknd, 20732.262_core_rknd, 21708.352_core_rknd, &
81 : 22723.592_core_rknd, 23779.273_core_rknd, 24876.709_core_rknd, 26017.258_core_rknd, &
82 : 27202.3_core_rknd, 28433.256_core_rknd, 29711.578_core_rknd, 31038.766_core_rknd /)
83 : !$acc declare create( svp_liq_lookup_table )
84 :
85 :
86 : contains
87 :
88 : !-------------------------------------------------------------------------
89 : ! Wrapped in interface sat_mixrat_liq
90 6368645631 : function sat_mixrat_liq_k( p_in_Pa, T_in_K, &
91 : saturation_formula )
92 : !$acc routine seq
93 : ! Description:
94 : ! Used to compute the saturation mixing ratio of liquid water.
95 :
96 : ! References:
97 : ! Formula from Emanuel 1994, 4.4.14
98 : !-------------------------------------------------------------------------
99 :
100 : use constants_clubb, only: &
101 : ep
102 :
103 : use constants_clubb, only: &
104 : T_freeze_K
105 :
106 : implicit none
107 :
108 : ! -------------------- Input Variables --------------------
109 : real( kind = core_rknd ), intent(in) :: &
110 : p_in_Pa, & ! Pressure [Pa]
111 : T_in_K ! Temperature [K]
112 :
113 : integer, intent(in) :: &
114 : saturation_formula
115 :
116 : ! -------------------- Output Variables --------------------
117 : real( kind = core_rknd ) :: &
118 : sat_mixrat_liq_k
119 :
120 : ! -------------------- Local Variables --------------------
121 : real( kind = core_rknd ) :: &
122 : T_in_C, &
123 : T_in_C_sqd, &
124 : T_in_K_clipped
125 :
126 : integer :: &
127 : T_in_K_int
128 :
129 : ! Constant parameters
130 :
131 : ! Relative error norm expansion (-50 to 50 deg_C) from
132 : ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
133 : ! (The 100 coefficient converts from mb to Pa)
134 : ! real, dimension(7), parameter :: a = &
135 : ! 100.* (/ 6.11176750, 0.443986062, 0.143053301E-01, &
136 : ! 0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
137 : ! 0.638780966E-10 /)
138 :
139 : ! Relative error norm expansion (-85 to 70 deg_C) from
140 : ! Table 4 of pp. 1511 of Flatau et al.
141 : !real( kind = core_rknd ), dimension(9), parameter :: a = &
142 : !100._core_rknd * &
143 : ! Commented out because the form has been redone, causing these number to no longer be needed,
144 : ! leaving them in for now for reference.
145 : ! (/ 6.11583699_core_rknd, 0.444606896_core_rknd, 0.143177157E-01_core_rknd, &
146 : ! 0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
147 : ! 0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
148 :
149 : real( kind = core_rknd ), parameter :: &
150 : min_T_in_C = -85._core_rknd, & ! [deg_C]
151 : min_T_in_K = 173.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
152 :
153 : ! ---------------------- Output Variables ----------------------
154 : real( kind = core_rknd ) :: &
155 : esat ! Saturation vapor pressure over water [Pa]
156 :
157 : ! -------------------- Begin Code --------------------
158 :
159 : ! Calculate the SVP for water vapor.
160 6368645631 : select case ( saturation_formula )
161 : case ( saturation_flatau )
162 :
163 : ! Using the Flatau, et al. polynomial approximation for SVP over vapor
164 :
165 : ! Determine deg K - 273.15
166 0 : T_in_C = T_in_K - T_freeze_K
167 :
168 : ! Since this approximation is only good out to -85 degrees Celsius we
169 : ! truncate the result here (Flatau, et al. 1992)
170 0 : T_in_C = max( T_in_C, min_T_in_C )
171 :
172 : ! Polynomial approx. (Flatau, et al. 1992)
173 :
174 : ! This is the generalized formula but is not computationally efficient.
175 : ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
176 : ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
177 :
178 : ! esat = a(1)
179 :
180 : ! do i = 2, size( a ) , 1
181 : ! esat = esat + a(i) * ( T_in_C )**(i-1)
182 : ! end do
183 :
184 : ! The 8th order polynomial fit. When running deep
185 : ! convective cases I noticed that absolute temperature often dips below
186 : ! -50 deg_C at higher altitudes, where the 6th order approximation is
187 : ! not accurate. -dschanen 20 Nov 2008
188 : !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
189 : !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
190 :
191 :
192 : ! Factoring the polynomial above and changing it into this form allows the cpu
193 : ! to complete the calculations out of order. This is because modern cpus can complete
194 : ! multiple instructions at once if they do not depend on eachother, in the above case
195 : ! each instruction relies on the result of the last. In this version however, the terms
196 : ! in the parentheses could potentially be calculated in parallel by different execution
197 : ! units in the cpu, then only when those terms are being multiplied together do the
198 : ! instructions need to be done one at a time. See clubb issue 834 for more info.
199 : ! - Gunther Huebler, Aug 2018
200 0 : T_in_C_sqd = T_in_C**2
201 :
202 : esat = &
203 : - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
204 : * ( T_in_C + 90.72381630364440_core_rknd ) &
205 : * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
206 : * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
207 0 : * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
208 :
209 : case ( saturation_bolton )
210 :
211 : ! Using the Bolton 1980 approximations for SVP over vapor
212 : ! Generally this more computationally expensive than the Flatau polnomial expansion
213 : esat = 611.2_core_rknd &
214 : * exp( (17.67_core_rknd *(T_in_K-T_freeze_K)) &
215 0 : / (T_in_K-29.65_core_rknd) ) ! Known magic number
216 :
217 : ! ---> h1g
218 : case ( saturation_gfdl )
219 :
220 : ! Using GFDL polynomial approximation for SVP with respect to liquid
221 : ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
222 : ! we threshold the temperature. This will yield a minimal saturation at
223 : ! cold temperatures.
224 6368645631 : T_in_K_clipped = max( min_T_in_K, T_in_K )
225 :
226 : ! Goff Gratch equation, uncertain below -70 C
227 :
228 : esat = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
229 : 5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
230 : 1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
231 : *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
232 : 8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
233 : *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
234 6368645631 : log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
235 :
236 : ! <--- h1g
237 :
238 : case ( saturation_lookup )
239 :
240 0 : T_in_K_int = int( anint( T_in_K ) )
241 :
242 : ! Since this approximation is only good out to -85 degrees Celsius we
243 : ! truncate the result here
244 0 : T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
245 :
246 : ! Use the lookup table to determine the saturation vapor pressure.
247 0 : esat = svp_liq_lookup_table( T_in_K_int )
248 :
249 : case default
250 :
251 : ! Undefined approximation
252 6368645631 : esat = -99999.999_core_rknd
253 :
254 : end select
255 :
256 : ! If esat exceeds the air pressure, then assume esat~=0.5*pressure
257 : ! and set rsat = ep = 0.622
258 6368645631 : if ( p_in_Pa-esat < 1.0_core_rknd ) then
259 : sat_mixrat_liq_k = ep
260 : else
261 :
262 6368645631 : if ( I_sat_sphum ) then ! h1g, 2010-06-18 begin mod
263 :
264 : ! GFDL uses specific humidity
265 : ! Formula for Saturation Specific Humidity
266 : sat_mixrat_liq_k = ep * ( esat / ( p_in_Pa &
267 : - (1.0_core_rknd-ep) * esat ) )
268 : else
269 :
270 : ! Formula for Saturation Mixing Ratio:
271 : !
272 : ! rs = (epsilon) * [ esat / ( p - esat ) ];
273 : ! where epsilon = R_d / R_v
274 6368645631 : sat_mixrat_liq_k = ep * esat / ( p_in_Pa - esat )
275 :
276 : endif ! h1g, 2010-06-18 end mod
277 :
278 : end if
279 :
280 : return
281 :
282 : end function sat_mixrat_liq_k
283 :
284 : !-------------------------------------------------------------------------
285 : !
286 80415504 : function sat_mixrat_liq_2D( nz, ngrdcol, p_in_Pa, T_in_K, &
287 : saturation_formula, &
288 : start_index_in )
289 : !
290 : ! Description:
291 : ! Used to compute the saturation mixing ratio of liquid water.
292 : !
293 : ! References:
294 : ! Formula from Emanuel 1994, 4.4.14
295 : !-------------------------------------------------------------------------
296 :
297 : use constants_clubb, only: &
298 : ep
299 :
300 : use constants_clubb, only: T_freeze_K
301 :
302 : implicit none
303 :
304 : ! -------------------- Input Variables --------------------
305 : integer, intent(in) :: &
306 : nz, &
307 : ngrdcol
308 :
309 : integer, intent(in) :: &
310 : saturation_formula
311 :
312 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
313 : p_in_Pa, & ! Pressure [Pa]
314 : T_in_K ! Temperature [K]
315 :
316 : ! -------------------- Output Variables --------------------
317 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
318 : sat_mixrat_liq_2D
319 :
320 : ! -------------------- Optional Input --------------------
321 : integer, intent(in), optional :: &
322 : start_index_in
323 :
324 : ! -------------------- Local Variables --------------------
325 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
326 80415504 : esat
327 :
328 : real( kind = core_rknd ) :: &
329 : T_in_C, &
330 : T_in_C_sqd, &
331 : T_in_K_clipped
332 :
333 : integer :: &
334 : T_in_K_int
335 :
336 : integer :: i,k
337 :
338 : ! Constant parameters
339 :
340 : ! Relative error norm expansion (-50 to 50 deg_C) from
341 : ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
342 : ! (The 100 coefficient converts from mb to Pa)
343 : ! real, dimension(7), parameter :: a = &
344 : ! 100.* (/ 6.11176750, 0.443986062, 0.143053301E-01, &
345 : ! 0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
346 : ! 0.638780966E-10 /)
347 :
348 : ! Relative error norm expansion (-85 to 70 deg_C) from
349 : ! Table 4 of pp. 1511 of Flatau et al.
350 : !real( kind = core_rknd ), dimension(9), parameter :: a = &
351 : !100._core_rknd * &
352 : ! Commented out because the form has been redone, causing these number to no longer be needed,
353 : ! leaving them in for now for reference.
354 : ! (/ 6.11583699_core_rknd, 0.444606896_core_rknd, 0.143177157E-01_core_rknd, &
355 : ! 0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
356 : ! 0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
357 :
358 : real( kind = core_rknd ), parameter :: &
359 : min_T_in_C = -85._core_rknd, & ! [deg_C]
360 : min_T_in_K = 173.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
361 :
362 : integer :: &
363 : start_index
364 :
365 : ! -------------------- Begin Code --------------------
366 :
367 : !$acc data create(esat) &
368 : !$acc copyin(p_in_Pa,T_in_K), &
369 : !$acc copyout(sat_mixrat_liq_2D)
370 :
371 : ! start_index is an optional argument and
372 : ! used for choosing the sub-arrays
373 80415504 : if ( present(start_index_in) ) then
374 17870112 : start_index = start_index_in
375 : else
376 : start_index = 1
377 : end if
378 :
379 80415504 : select case ( saturation_formula )
380 : case ( saturation_flatau )
381 :
382 : !$acc parallel loop gang vector collapse(2) default(present)
383 0 : do k = start_index, nz
384 0 : do i = 1, ngrdcol
385 :
386 : ! Determine deg K - 273.15
387 0 : T_in_C = T_in_K(i,k) - T_freeze_K
388 :
389 : ! Since this approximation is only good out to -85 degrees Celsius we
390 : ! truncate the result here (Flatau, et al. 1992)
391 0 : T_in_C = max( T_in_C, min_T_in_C )
392 :
393 : ! Polynomial approx. (Flatau, et al. 1992)
394 :
395 : ! This is the generalized formula but is not computationally efficient.
396 : ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
397 : ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
398 :
399 : ! esat = a(1)
400 :
401 : ! do i = 2, size( a ) , 1
402 : ! esat = esat + a(i) * ( T_in_C )**(i-1)
403 : ! end do
404 :
405 : ! The 8th order polynomial fit. When running deep
406 : ! convective cases I noticed that absolute temperature often dips below
407 : ! -50 deg_C at higher altitudes, where the 6th order approximation is
408 : ! not accurate. -dschanen 20 Nov 2008
409 : !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
410 : !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
411 :
412 :
413 : ! Factoring the polynomial above and changing it into this form allows the cpu
414 : ! to complete the calculations out of order. This is because modern cpus can complete
415 : ! multiple instructions at once if they do not depend on eachother, in the above case
416 : ! each instruction relies on the result of the last. In this version however, the terms
417 : ! in the parentheses could potentially be calculated in parallel by different execution
418 : ! units in the cpu, then only when those terms are being multiplied together do the
419 : ! instructions need to be done one at a time. See clubb issue 834 for more info.
420 : ! - Gunther Huebler, Aug 2018
421 0 : T_in_C_sqd = T_in_C**2
422 :
423 : esat(i,k) = &
424 : - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
425 : * ( T_in_C + 90.72381630364440_core_rknd ) &
426 : * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
427 : * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
428 0 : * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
429 : end do
430 : end do
431 : !$acc end parallel loop
432 :
433 : case ( saturation_bolton )
434 :
435 : ! Using the Bolton 1980 approximations for SVP over vapor
436 : ! Generally this more computationally expensive than the Flatau polnomial expansion
437 : !$acc parallel loop gang vector collapse(2) default(present)
438 0 : do k = start_index, nz
439 0 : do i = 1, ngrdcol
440 0 : esat(i,k) = 611.2_core_rknd &
441 : * exp( (17.67_core_rknd *(T_in_K(i,k)-T_freeze_K)) &
442 0 : / (T_in_K(i,k)-29.65_core_rknd) ) ! Known magic number
443 : end do
444 : end do
445 : !$acc end parallel loop
446 :
447 : ! ---> h1g
448 : case ( saturation_gfdl )
449 :
450 : ! Using GFDL polynomial approximation for SVP with respect to liquid
451 : !$acc parallel loop gang vector collapse(2) default(present)
452 6888928176 : do k = start_index, nz
453 >11376*10^7 : do i = 1, ngrdcol
454 :
455 : ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
456 : ! we threshold the temperature. This will yield a minimal saturation at
457 : ! cold temperatures.
458 >10687*10^7 : T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
459 :
460 : ! Goff Gratch equation, uncertain below -70 C
461 :
462 : esat(i,k) = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
463 : 5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
464 : 1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
465 : *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
466 : 8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
467 : *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
468 >11368*10^7 : log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
469 : end do
470 : end do
471 : !$acc end parallel loop
472 :
473 : ! <--- h1g
474 :
475 : case ( saturation_lookup )
476 :
477 : !$acc parallel loop gang vector collapse(2) default(present)
478 0 : do k = start_index, nz
479 0 : do i = 1, ngrdcol
480 0 : T_in_K_int = int( anint( T_in_K(i,k) ) )
481 :
482 : ! Since this approximation is only good out to -85 degrees Celsius we
483 : ! truncate the result here
484 0 : T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
485 :
486 : ! Use the lookup table to determine the saturation vapor pressure.
487 0 : esat(i,k) = svp_liq_lookup_table( T_in_K_int )
488 : end do
489 : end do
490 : !$acc end parallel loop
491 :
492 : case default
493 :
494 : ! Undefined approximation
495 80415504 : esat = -99999.999_core_rknd
496 :
497 : end select
498 :
499 : !$acc parallel loop gang vector collapse(2) default(present)
500 6888928176 : do k = start_index, nz
501 >11376*10^7 : do i = 1, ngrdcol
502 :
503 : ! If esat exceeds the air pressure, then assume esat~=0.5*pressure
504 : ! and set rsat = ep = 0.622
505 >11368*10^7 : if ( p_in_Pa(i,k)-esat(i,k) < 1.0_core_rknd ) then
506 2040744207 : sat_mixrat_liq_2D(i,k) = ep
507 : else
508 :
509 : if ( I_sat_sphum ) then ! h1g, 2010-06-18 begin mod
510 :
511 : ! GFDL uses specific humidity
512 : ! Formula for Saturation Specific Humidity
513 : sat_mixrat_liq_2D(i,k) = ep * ( esat(i,k) &
514 : / ( p_in_Pa(i,k) - (1.0_core_rknd-ep * esat(i,k) ) ) )
515 : else
516 :
517 : ! Formula for Saturation Mixing Ratio:
518 : !
519 : ! rs = (epsilon) * [ esat / ( p - esat ) ];
520 : ! where epsilon = R_d / R_v
521 >10483*10^7 : sat_mixrat_liq_2D(i,k) = ep * esat(i,k) / ( p_in_Pa(i,k) - esat(i,k) )
522 :
523 : endif ! h1g, 2010-06-18 end mod
524 :
525 : end if
526 :
527 : end do
528 : end do
529 : !$acc end parallel loop
530 :
531 : !$acc end data
532 :
533 : end function sat_mixrat_liq_2D
534 :
535 : !-----------------------------------------------------------------
536 : ! Wrapped in interface sat_vapor_press_liq
537 0 : subroutine sat_vapor_press_liq( T_in_K, &
538 : saturation_formula, &
539 : esat )
540 :
541 : ! Description:
542 : ! Computes SVP for water vapor. Calls one of the other functions
543 : ! that calculate an approximation to SVP.
544 :
545 : ! References:
546 : ! None
547 : !-----------------------------------------------------------------
548 :
549 : use clubb_precision, only: &
550 : core_rknd ! Variable(s)
551 :
552 : implicit none
553 :
554 : ! ------------------------ Input Variables ------------------------
555 : real( kind = core_rknd ), intent(in) :: &
556 : T_in_K ! Temperature [K]
557 :
558 : integer, intent(in) :: &
559 : saturation_formula
560 :
561 : ! ------------------------ Output Variables ------------------------
562 : real( kind = core_rknd ), intent(out) :: &
563 : esat ! Saturation Vapor Pressure over Water [Pa]
564 :
565 : ! ------------------------ Being Code ------------------------
566 :
567 : ! Saturation Vapor Pressure, esat, can be found to be approximated
568 : ! in many different ways.
569 0 : select case ( saturation_formula )
570 : case ( saturation_bolton )
571 :
572 : ! Using the Bolton 1980 approximations for SVP over vapor
573 : call sat_vapor_press_liq_bolton( T_in_K, &
574 0 : esat )
575 :
576 : !Earthworks case
577 : case ( saturation_flatau )
578 :
579 : ! Using the Flatau, et al. polynomial approximation for SVP over vapor
580 : call sat_vapor_press_liq_flatau( T_in_K, &
581 0 : esat )
582 :
583 : ! ---> h1g
584 : case ( saturation_gfdl )
585 :
586 : ! Using GFDL polynomial approximation for SVP with respect to liquid
587 : call sat_vapor_press_liq_gfdl( T_in_K, &
588 0 : esat )
589 :
590 : ! <--- h1g
591 : case ( saturation_lookup )
592 :
593 : ! Use the lookup table to determine the saturation vapor pressure.
594 0 : esat = sat_vapor_press_liq_lookup( T_in_K )
595 :
596 : case default
597 :
598 : ! Undefined approximation
599 0 : esat = -99999.999_core_rknd
600 :
601 : end select
602 :
603 0 : return
604 :
605 : end subroutine sat_vapor_press_liq
606 :
607 : !------------------------------------------------------------------------
608 0 : subroutine sat_vapor_press_liq_flatau( T_in_K, &
609 : esat )
610 :
611 : ! Description:
612 : ! Computes SVP for water vapor.
613 :
614 : ! References:
615 : ! ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
616 : ! and Cotton. (1992) Journal of Applied Meteorology, Vol. 31,
617 : ! pp. 1507--1513
618 : !------------------------------------------------------------------------
619 :
620 : use constants_clubb, only: T_freeze_K
621 :
622 : use clubb_precision, only: &
623 : core_rknd ! Variable(s)
624 :
625 : implicit none
626 :
627 : ! Constant parameters
628 :
629 : ! Relative error norm expansion (-50 to 50 deg_C) from
630 : ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
631 : ! (The 100 coefficient converts from mb to Pa)
632 : ! real, dimension(7), parameter :: a = &
633 : ! 100.* (/ 6.11176750, 0.443986062, 0.143053301E-01, &
634 : ! 0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
635 : ! 0.638780966E-10 /)
636 :
637 : ! Relative error norm expansion (-85 to 70 deg_C) from
638 : ! Table 4 of pp. 1511 of Flatau et al.
639 : !real( kind = core_rknd ), dimension(9), parameter :: a = &
640 : !100._core_rknd * &
641 : ! Commented out because the form has been redone, causing these number to no longer be needed,
642 : ! leaving them in for now for reference.
643 : ! (/ 6.11583699_core_rknd, 0.444606896_core_rknd, 0.143177157E-01_core_rknd, &
644 : ! 0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
645 : ! 0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
646 :
647 : real( kind = core_rknd ), parameter :: min_T_in_C = -85._core_rknd ! [deg_C]
648 :
649 : ! ---------------------- Input Variables ----------------------
650 : real( kind = core_rknd ), intent(in) :: &
651 : T_in_K ! Temperature [K]
652 :
653 : ! ---------------------- Output Variables ----------------------
654 : real( kind = core_rknd ), intent(out) :: &
655 : esat ! Saturation vapor pressure over water [Pa]
656 :
657 : ! ---------------------- Local Variables ----------------------
658 : real( kind = core_rknd ) :: T_in_C, T_in_C_sqd
659 :
660 : ! ---------------------- Begin Code ----------------------
661 :
662 : ! Determine deg K - 273.15
663 0 : T_in_C = T_in_K - T_freeze_K
664 :
665 : ! Since this approximation is only good out to -85 degrees Celsius we
666 : ! truncate the result here (Flatau, et al. 1992)
667 0 : T_in_C = max( T_in_C, min_T_in_C )
668 :
669 : ! Polynomial approx. (Flatau, et al. 1992)
670 :
671 : ! This is the generalized formula but is not computationally efficient.
672 : ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
673 : ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
674 :
675 : ! esat = a(1)
676 :
677 : ! do i = 2, size( a ) , 1
678 : ! esat = esat + a(i) * ( T_in_C )**(i-1)
679 : ! end do
680 :
681 : ! The 8th order polynomial fit. When running deep
682 : ! convective cases I noticed that absolute temperature often dips below
683 : ! -50 deg_C at higher altitudes, where the 6th order approximation is
684 : ! not accurate. -dschanen 20 Nov 2008
685 : !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
686 : !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
687 :
688 :
689 : ! Factoring the polynomial above and changing it into this form allows the cpu
690 : ! to complete the calculations out of order. This is because modern cpus can complete
691 : ! multiple instructions at once if they do not depend on eachother, in the above case
692 : ! each instruction relies on the result of the last. In this version however, the terms
693 : ! in the parentheses could potentially be calculated in parallel by different execution
694 : ! units in the cpu, then only when those terms are being multiplied together do the
695 : ! instructions need to be done one at a time. See clubb issue 834 for more info.
696 : ! - Gunther Huebler, Aug 2018
697 0 : T_in_C_sqd = T_in_C**2
698 :
699 : esat = &
700 : - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
701 : * ( T_in_C + 90.72381630364440_core_rknd ) &
702 : * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
703 : * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
704 0 : * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
705 :
706 0 : return
707 : end subroutine sat_vapor_press_liq_flatau
708 :
709 :
710 : !------------------------------------------------------------------------
711 0 : subroutine sat_vapor_press_liq_bolton( T_in_K, &
712 : esat )
713 : ! Description:
714 : ! Computes SVP for water vapor.
715 : ! References:
716 : ! Bolton 1980
717 : !------------------------------------------------------------------------
718 :
719 : use constants_clubb, only: T_freeze_K
720 :
721 : use clubb_precision, only: &
722 : core_rknd ! Variable(s)
723 :
724 : implicit none
725 :
726 : ! --------------------- Input Variables ---------------------
727 : real( kind = core_rknd ), intent(in) :: &
728 : T_in_K ! Temperature [K]
729 :
730 : ! --------------------- Output Variables ---------------------
731 : real( kind = core_rknd ), intent(out) :: &
732 : esat ! Saturation vapor pressure over water [Pa]
733 :
734 :
735 : ! --------------------------- Begin Code ---------------------------
736 :
737 : ! (Bolton 1980) approx.
738 : ! Generally this more computationally expensive than the Flatau polnomial expansion
739 : esat = 611.2_core_rknd &
740 : * exp( (17.67_core_rknd *(T_in_K-T_freeze_K)) &
741 0 : / (T_in_K-29.65_core_rknd) ) ! Known magic number
742 :
743 0 : return
744 : end subroutine sat_vapor_press_liq_bolton
745 :
746 :
747 : ! ---> h1g, 2010-06-16
748 : !------------------------------------------------------------------------
749 0 : subroutine sat_vapor_press_liq_gfdl( T_in_K, &
750 : esat )
751 : ! Description:
752 : ! copy from "GFDL polysvp.F90"
753 : ! Compute saturation vapor pressure with respect to liquid by using
754 : ! function from Goff and Gratch (1946)
755 :
756 : ! Polysvp returned in units of pa.
757 : ! T_in_K is input in units of K.
758 : !------------------------------------------------------------------------
759 :
760 : use clubb_precision, only: &
761 : core_rknd ! Variable(s)
762 :
763 : implicit none
764 :
765 : ! --------------------------- Input Variables ---------------------------
766 : real( kind = core_rknd ), intent(in) :: &
767 : T_in_K ! Absolute temperature [K]
768 :
769 : ! --------------------------- Output Variables ---------------------------
770 : real( kind = core_rknd ), intent(out) :: &
771 : esat ! Saturation vapor pressure over water [Pa]
772 :
773 : ! --------------------------- Local Variables ---------------------------
774 : real( kind = core_rknd ), parameter :: &
775 : min_T_in_K = 203.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
776 :
777 : real( kind = core_rknd ) :: &
778 : T_in_K_clipped ! Absolute temperature with minimum threshold applied [K]
779 :
780 : ! --------------------------- Begin Code ---------------------------
781 :
782 : ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
783 : ! we threshold the temperature. This will yield a minimal saturation at
784 : ! cold temperatures.
785 0 : T_in_K_clipped = max( min_T_in_K, T_in_K )
786 :
787 : ! Goff Gratch equation, uncertain below -70 C
788 :
789 : esat = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
790 : 5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
791 : 1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
792 : *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
793 : 8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
794 : *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
795 0 : log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
796 :
797 0 : return
798 : end subroutine sat_vapor_press_liq_gfdl
799 : ! <--- h1g, 2010-06-16
800 :
801 : !------------------------------------------------------------------------
802 0 : elemental function sat_vapor_press_liq_lookup( T_in_K ) result ( esat )
803 :
804 : ! Description:
805 : ! Computes SVP for water vapor, using a lookup table.
806 : !
807 : ! The lookup table was constructed using the Flatau approximation.
808 :
809 : ! References:
810 : ! ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
811 : ! and Cotton. (1992) Journal of Applied Meteorology, Vol. 31,
812 : ! pp. 1507--1513
813 : !------------------------------------------------------------------------
814 :
815 : implicit none
816 :
817 : ! External
818 : intrinsic :: max, min, int, anint
819 :
820 : ! Input Variables
821 : real( kind = core_rknd ), intent(in) :: T_in_K ! Temperature [K]
822 :
823 : ! Output Variables
824 : real( kind = core_rknd ) :: esat ! Saturation vapor pressure over water [Pa]
825 :
826 : ! Local Variables
827 : integer :: T_in_K_int
828 :
829 : ! ---- Begin Code ----
830 :
831 0 : T_in_K_int = int( anint( T_in_K ) )
832 :
833 : ! Since this approximation is only good out to -85 degrees Celsius we
834 : ! truncate the result here
835 0 : T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
836 :
837 : ! Use the lookup table to determine the saturation vapor pressure.
838 0 : esat = svp_liq_lookup_table( T_in_K_int )
839 :
840 : return
841 : end function sat_vapor_press_liq_lookup
842 :
843 : !------------------------------------------------------------------------
844 : ! Wrapped in interface sat_mixrat_ice
845 0 : function sat_mixrat_ice_k( p_in_Pa, T_in_K, saturation_formula )
846 :
847 : ! Description:
848 : ! Used to compute the saturation mixing ratio of ice.
849 :
850 : ! References:
851 : ! Formula from Emanuel 1994, 4.4.15
852 : !-------------------------------------------------------------------------
853 :
854 : use constants_clubb, only: &
855 : ep ! Variable(s)
856 :
857 : use clubb_precision, only: &
858 : core_rknd ! Variable(s)
859 :
860 : implicit none
861 :
862 : ! ------------------------ Input Variables ------------------------
863 : real( kind = core_rknd ), intent(in) :: &
864 : p_in_Pa, & ! Pressure [Pa]
865 : T_in_K ! Temperature [K]
866 :
867 : integer, intent(in) :: &
868 : saturation_formula
869 :
870 : ! ------------------------ Output Variables ------------------------
871 : real( kind = core_rknd ) :: &
872 : sat_mixrat_ice_k
873 :
874 : ! ------------------------ Local Variables ------------------------
875 : real( kind = core_rknd ), dimension(1,1) :: &
876 : p_in_Pa_col, &
877 : T_in_K_col
878 :
879 : real( kind = core_rknd ), dimension(1,1) :: &
880 : sat_mixrat_ice_col
881 :
882 : integer :: i, k ! Loop indices
883 :
884 : ! ------------------------ Begin Code ------------------------
885 :
886 : ! Copy inputs to 2D arrays
887 0 : p_in_Pa_col(1,1) = p_in_Pa
888 0 : T_in_K_col(1,1) = T_in_K
889 :
890 : ! Call 2D version
891 0 : sat_mixrat_ice_col = sat_mixrat_ice_2D( 1, 1, p_in_Pa_col, T_in_K_col, saturation_formula )
892 :
893 : ! Copy 2D result into output
894 0 : sat_mixrat_ice_k = sat_mixrat_ice_col( 1, 1 )
895 :
896 : return
897 : end function sat_mixrat_ice_k
898 :
899 : !------------------------------------------------------------------------
900 : ! Wrapped in interface sat_mixrat_ice
901 35740224 : function sat_mixrat_ice_2D( nz, ngrdcol, p_in_Pa, T_in_K, saturation_formula )
902 :
903 : ! Description:
904 : ! Used to compute the saturation mixing ratio of ice.
905 :
906 : ! References:
907 : ! Formula from Emanuel 1994, 4.4.15
908 : !-------------------------------------------------------------------------
909 :
910 : use constants_clubb, only: &
911 : T_freeze_K, & ! Variable(s)
912 : ep
913 :
914 : use clubb_precision, only: &
915 : core_rknd ! Variable(s)
916 :
917 : implicit none
918 :
919 : ! ------------------------ Input Variables ------------------------
920 : integer, intent(in) :: &
921 : nz, &
922 : ngrdcol
923 :
924 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
925 : p_in_Pa, & ! Pressure [Pa]
926 : T_in_K ! Temperature [K]
927 :
928 : integer, intent(in) :: &
929 : saturation_formula
930 :
931 : ! ------------------------ Output Variables ------------------------
932 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
933 : sat_mixrat_ice_2D
934 :
935 : ! ------------------------ Local Variables ------------------------
936 : ! Relative error norm expansion (-90 to 0 deg_C) from
937 : ! Table 4 of pp. 1511 of Flatau et al. 1992 (Ice)
938 : real( kind = core_rknd ), dimension(9), parameter :: a = &
939 : 100._core_rknd * (/ 6.09868993_core_rknd, 0.499320233_core_rknd, 0.184672631E-01_core_rknd, &
940 : 0.402737184E-03_core_rknd, 0.565392987E-05_core_rknd, 0.521693933E-07_core_rknd, &
941 : 0.307839583E-09_core_rknd, 0.105785160E-11_core_rknd, 0.161444444E-14_core_rknd /)
942 :
943 : real( kind = core_rknd ), parameter :: &
944 : min_T_in_C = -90._core_rknd, & ! [deg_C]
945 : min_T_in_K = 173.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
946 :
947 : real( kind = core_rknd ) :: &
948 : T_in_C, & ! Temperature [deg_C]
949 : T_in_K_clipped ! Absolute temperature with minimum threshold applied [K]
950 :
951 : real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
952 35740224 : esat_ice
953 :
954 : integer :: i, k ! Loop indices
955 :
956 : ! ------------------------ Begin Code ------------------------
957 :
958 : !$acc data create( esat_ice ) &
959 : !$acc copyin( p_in_Pa, T_in_K ) &
960 : !$acc copyout( sat_mixrat_ice_2D )
961 :
962 : ! Determine the SVP for the given temperature
963 35740224 : select case ( saturation_formula )
964 : case ( saturation_bolton )
965 :
966 : ! Using the Bolton 1980 approximations for SVP over ice
967 : !$acc parallel loop gang vector collapse(2) default(present)
968 0 : do i = 1, ngrdcol
969 0 : do k = 1, nz
970 :
971 : ! Exponential approx.
972 0 : esat_ice(i,k) = 100.0_core_rknd * exp( 23.33086_core_rknd - &
973 0 : (6111.72784_core_rknd/T_in_K(i,k)) + (0.15215_core_rknd*log( T_in_K(i,k) )) )
974 :
975 : end do
976 : end do
977 : !$acc end parallel loop
978 :
979 : case ( saturation_flatau )
980 :
981 : ! Using the Flatau, et al. polynomial approximation for SVP over ice
982 : !$acc parallel loop gang vector collapse(2) default(present)
983 0 : do i = 1, ngrdcol
984 0 : do k = 1, nz
985 :
986 : ! Determine deg K - 273.15
987 0 : T_in_C = T_in_K(i,k) - T_freeze_K
988 :
989 : ! Since this approximation is only good out to -90 degrees Celsius we
990 : ! truncate the result here (Flatau, et al. 1992)
991 0 : T_in_C = max( T_in_C, min_T_in_C )
992 :
993 : ! Polynomial approx. (Flatau, et al. 1992)
994 : ! esati = a(1)
995 :
996 : ! do i = 2, size( a ), 1
997 : ! esati = esati + a(i) * ( T_in_C )**(i-1)
998 : ! end do
999 :
1000 : esat_ice(i,k) = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
1001 0 : *( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*a(9) ) ) ) ) ) ) )
1002 :
1003 : end do
1004 : end do
1005 : !$acc end parallel loop
1006 :
1007 : ! ---> h1g, 2010-06-16
1008 : case ( saturation_gfdl )
1009 :
1010 : ! Using GFDL polynomial approximation for SVP with respect to ice
1011 : !$acc parallel loop gang vector collapse(2) default(present)
1012 596778624 : do i = 1, ngrdcol
1013 48285042624 : do k = 1, nz
1014 :
1015 : ! Since the Goff-Gratch ice approximation is valid only down to -100 degrees Celsius,
1016 : ! we threshold the temperature. This will yield a minimal saturation at
1017 : ! cold temperatures.
1018 47688264000 : T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
1019 :
1020 : ! Goff Gratch equation (good down to -100 C)
1021 :
1022 : esat_ice(i,k) = 10._core_rknd**(-9.09718_core_rknd* &
1023 : (273.16_core_rknd/T_in_K_clipped-1._core_rknd)-3.56654_core_rknd* &
1024 : log10(273.16_core_rknd/T_in_K_clipped)+0.876793_core_rknd* &
1025 : (1._core_rknd-T_in_K_clipped/273.16_core_rknd)+ &
1026 48249302400 : log10(6.1071_core_rknd))*100._core_rknd ! Known magic number
1027 : end do
1028 : end do
1029 : !$acc end parallel loop
1030 :
1031 : ! <--- h1g, 2010-06-16
1032 :
1033 : case default
1034 :
1035 : ! Undefined approximation
1036 35740224 : esat_ice = -99999.999_core_rknd
1037 :
1038 : end select
1039 :
1040 :
1041 : !$acc parallel loop gang vector collapse(2) default(present)
1042 3073659264 : do k = 1, nz
1043 50761923264 : do i = 1, ngrdcol
1044 :
1045 : ! If esat_ice exceeds the air pressure, then assume esat_ice~=0.5*pressure
1046 : ! and set rsat = ep = 0.622
1047 50726183040 : if ( p_in_Pa(i,k)-esat_ice(i,k) < 1.0_core_rknd ) then
1048 890022639 : sat_mixrat_ice_2D(i,k) = ep
1049 : else
1050 :
1051 : if ( I_sat_sphum ) then ! h1g, 2010-06-18 begin mod
1052 :
1053 : ! GFDL uses specific humidity
1054 : ! Formula for Saturation Specific Humidity
1055 : sat_mixrat_ice_2D(i,k) = ep * ( esat_ice(i,k) &
1056 : / ( p_in_Pa(i,k) - (1.0_core_rknd-ep) * esat_ice(i,k) ) )
1057 : else
1058 :
1059 : ! Formula for Saturation Mixing Ratio:
1060 : !
1061 : ! rs = (epsilon) * [ esat / ( p - esat ) ];
1062 : ! where epsilon = R_d / R_v
1063 46798241361 : sat_mixrat_ice_2D(i,k) = ep * esat_ice(i,k) / ( p_in_Pa(i,k) - esat_ice(i,k) )
1064 :
1065 : endif ! h1g, 2010-06-18 end mod
1066 : end if
1067 :
1068 : end do
1069 : end do
1070 : !$acc end parallel loop
1071 :
1072 : !$acc end data
1073 :
1074 35740224 : return
1075 : end function sat_mixrat_ice_2D
1076 :
1077 : !------------------------------------------------------------------------
1078 : subroutine sat_vapor_press_ice( nz, ngrdcol, T_in_K, saturation_formula, &
1079 : esat_ice )
1080 : !
1081 : ! Description:
1082 : ! Computes SVP for ice, using one of the various approximations.
1083 : !
1084 : ! References:
1085 : ! None
1086 : !------------------------------------------------------------------------
1087 :
1088 :
1089 : use clubb_precision, only: &
1090 : core_rknd ! Variable(s)
1091 :
1092 : implicit none
1093 :
1094 : ! ---------------------- Input Variable ----------------------
1095 : integer, intent(in) :: &
1096 : nz, &
1097 : ngrdcol
1098 :
1099 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1100 : T_in_K ! Temperature [K]
1101 :
1102 : integer, intent(in) :: &
1103 : saturation_formula
1104 :
1105 : ! ---------------------- Output Variable ----------------------
1106 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1107 : esat_ice ! Saturation Vapor Pressure over Ice [Pa]
1108 :
1109 : ! ---------------------- Begin Code ----------------------
1110 :
1111 : select case ( saturation_formula )
1112 : case ( saturation_bolton )
1113 :
1114 : ! Using the Bolton 1980 approximations for SVP over ice
1115 : call sat_vapor_press_ice_bolton( nz, ngrdcol, T_in_K, &
1116 : esat_ice )
1117 :
1118 : case ( saturation_flatau )
1119 :
1120 : ! Using the Flatau, et al. polynomial approximation for SVP over ice
1121 : call sat_vapor_press_ice_flatau( nz, ngrdcol, T_in_K, &
1122 : esat_ice )
1123 :
1124 : ! ---> h1g, 2010-06-16
1125 : case ( saturation_gfdl )
1126 :
1127 : ! Using GFDL polynomial approximation for SVP with respect to ice
1128 : call sat_vapor_press_ice_gfdl( nz, ngrdcol, T_in_K, &
1129 : esat_ice )
1130 : ! <--- h1g, 2010-06-16
1131 :
1132 : case default
1133 :
1134 : ! Undefined approximation
1135 : esat_ice = -99999.999_core_rknd
1136 :
1137 : end select
1138 :
1139 : return
1140 :
1141 : end subroutine sat_vapor_press_ice
1142 :
1143 : !------------------------------------------------------------------------
1144 : subroutine sat_vapor_press_ice_flatau( nz, ngrdcol, T_in_K, &
1145 : esati )
1146 : !
1147 : ! Description:
1148 : ! Computes SVP for ice.
1149 : !
1150 : ! References:
1151 : ! ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
1152 : ! and Cotton. (1992) Journal of Applied Meteorology, Vol. 31,
1153 : ! pp. 1507--1513
1154 : !------------------------------------------------------------------------
1155 : use constants_clubb, only: T_freeze_K
1156 :
1157 : use clubb_precision, only: &
1158 : core_rknd ! Variable(s)
1159 :
1160 : implicit none
1161 :
1162 : ! ------------------------ Input Variables ------------------------
1163 : integer, intent(in) :: &
1164 : nz, &
1165 : ngrdcol
1166 :
1167 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1168 : T_in_K ! Temperature [deg_K]
1169 :
1170 : ! ------------------------ Output Variables ------------------------
1171 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1172 : esati ! Saturation vapor pressure over ice [Pa]
1173 :
1174 : ! ------------------------ Local Variables ------------------------
1175 : ! Relative error norm expansion (-90 to 0 deg_C) from
1176 : ! Table 4 of pp. 1511 of Flatau et al. 1992 (Ice)
1177 : real( kind = core_rknd ), dimension(9), parameter :: a = &
1178 : 100._core_rknd * (/ 6.09868993_core_rknd, 0.499320233_core_rknd, 0.184672631E-01_core_rknd, &
1179 : 0.402737184E-03_core_rknd, 0.565392987E-05_core_rknd, 0.521693933E-07_core_rknd, &
1180 : 0.307839583E-09_core_rknd, 0.105785160E-11_core_rknd, 0.161444444E-14_core_rknd /)
1181 :
1182 : real( kind = core_rknd ), parameter :: min_T_in_C = -90._core_rknd ! [deg_C]
1183 :
1184 : real( kind = core_rknd ) :: &
1185 : T_in_C ! Temperature [deg_C]
1186 :
1187 : integer :: i, k
1188 :
1189 : ! ------------------------ Begin Code ------------------------
1190 :
1191 : do i = 1, ngrdcol
1192 : do k = 1, nz
1193 :
1194 : ! Determine deg K - 273.15
1195 : T_in_C = T_in_K(i,k) - T_freeze_K
1196 :
1197 : ! Since this approximation is only good out to -90 degrees Celsius we
1198 : ! truncate the result here (Flatau, et al. 1992)
1199 : T_in_C = max( T_in_C, min_T_in_C )
1200 :
1201 : ! Polynomial approx. (Flatau, et al. 1992)
1202 : ! esati = a(1)
1203 :
1204 : ! do i = 2, size( a ), 1
1205 : ! esati = esati + a(i) * ( T_in_C )**(i-1)
1206 : ! end do
1207 :
1208 : esati(i,k) = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
1209 : *( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*a(9) ) ) ) ) ) ) )
1210 :
1211 : end do
1212 : end do
1213 :
1214 : return
1215 :
1216 : end subroutine sat_vapor_press_ice_flatau
1217 :
1218 : !------------------------------------------------------------------------
1219 : subroutine sat_vapor_press_ice_bolton( nz, ngrdcol, T_in_K, &
1220 : esati )
1221 : !
1222 : ! Description:
1223 : ! Computes SVP for ice.
1224 : !
1225 : ! References:
1226 : ! Bolton 1980
1227 : !------------------------------------------------------------------------
1228 :
1229 : use clubb_precision, only: &
1230 : core_rknd ! Variable(s)
1231 :
1232 : implicit none
1233 :
1234 : ! ------------------------ Input Variables ------------------------
1235 : integer, intent(in) :: &
1236 : nz, &
1237 : ngrdcol
1238 :
1239 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1240 : T_in_K ! Temperature [deg_K]
1241 :
1242 : ! ------------------------ Output Variables ------------------------
1243 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1244 : esati ! Saturation vapor pressure over ice [Pa]
1245 :
1246 : ! ------------------------ Local Variables ------------------------
1247 : integer :: i, k
1248 :
1249 : ! ------------------------ Begin Code ------------------------
1250 :
1251 : do i = 1, ngrdcol
1252 : do k = 1, nz
1253 :
1254 : ! Exponential approx.
1255 : esati(i,k) = 100.0_core_rknd * exp( 23.33086_core_rknd - &
1256 : (6111.72784_core_rknd/T_in_K(i,k)) + (0.15215_core_rknd*log( T_in_K(i,k) )) )
1257 :
1258 : end do
1259 : end do
1260 :
1261 : return
1262 :
1263 : end subroutine sat_vapor_press_ice_bolton
1264 :
1265 :
1266 : ! ---> h1g, 2010-06-16
1267 : !------------------------------------------------------------------------
1268 : subroutine sat_vapor_press_ice_gfdl( nz, ngrdcol, T_in_K, &
1269 : esati )
1270 : ! Description:
1271 : ! copy from "GFDL polysvp.F90"
1272 : ! Compute saturation vapor pressure with respect to liquid by using
1273 : ! function from Goff and Gratch (1946)
1274 : !
1275 : ! Polysvp returned in units of pa.
1276 : ! T_in_K is input in units of K.
1277 : !------------------------------------------------------------------------
1278 :
1279 : use clubb_precision, only: &
1280 : core_rknd ! Variable(s)
1281 :
1282 : implicit none
1283 :
1284 : ! ------------------------ Input Variables ------------------------
1285 : integer, intent(in) :: &
1286 : nz, &
1287 : ngrdcol
1288 :
1289 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
1290 : T_in_K ! Temperature [deg_K]
1291 :
1292 : ! ------------------------ Output Variables ------------------------
1293 : real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
1294 : esati ! Saturation vapor pressure over ice [Pa]
1295 :
1296 : ! ------------------------ Local Variables ------------------------
1297 : real( kind = core_rknd ), parameter :: &
1298 : min_T_in_K = 173.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
1299 :
1300 : real( kind = core_rknd ) :: &
1301 : T_in_K_clipped ! Absolute temperature with minimum threshold applied [K]
1302 :
1303 : integer :: i, k
1304 :
1305 : ! ------------------------ Begin Code ------------------------
1306 :
1307 : do i = 1, ngrdcol
1308 : do k = 1, nz
1309 :
1310 : ! Since the Goff-Gratch ice approximation is valid only down to -100 degrees Celsius,
1311 : ! we threshold the temperature. This will yield a minimal saturation at
1312 : ! cold temperatures.
1313 : T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
1314 :
1315 : ! Goff Gratch equation (good down to -100 C)
1316 :
1317 : esati(i,k) = 10._core_rknd**(-9.09718_core_rknd* &
1318 : (273.16_core_rknd/T_in_K_clipped-1._core_rknd)-3.56654_core_rknd* &
1319 : log10(273.16_core_rknd/T_in_K_clipped)+0.876793_core_rknd* &
1320 : (1._core_rknd-T_in_K_clipped/273.16_core_rknd)+ &
1321 : log10(6.1071_core_rknd))*100._core_rknd ! Known magic number
1322 : end do
1323 : end do
1324 :
1325 : return
1326 :
1327 : end subroutine sat_vapor_press_ice_gfdl
1328 : ! <--- h1g, 2010-06-16
1329 :
1330 : !-------------------------------------------------------------------------
1331 0 : function rcm_sat_adj( thlm, rtm, p_in_Pa, exner, saturation_formula ) result ( rcm )
1332 :
1333 : ! Description:
1334 : !
1335 : ! This function uses an iterative method to find the value of rcm
1336 : ! from an initial profile that has saturation at some point.
1337 : !
1338 : ! References:
1339 : ! None
1340 : !-------------------------------------------------------------------------
1341 :
1342 : use clubb_precision, only: &
1343 : core_rknd ! Variable(s)
1344 :
1345 : use constants_clubb, only: &
1346 : Cp, & ! Variable(s)
1347 : Lv, &
1348 : zero_threshold
1349 :
1350 : implicit none
1351 :
1352 : ! Local Constant(s)
1353 : real( kind = core_rknd ), parameter :: &
1354 : tolerance = 0.001_core_rknd ! Tolerance on theta calculation [K]
1355 :
1356 : integer, parameter :: &
1357 : itermax = 1000000 ! Maximum interations
1358 :
1359 : ! External
1360 : intrinsic :: max, abs
1361 :
1362 : ! Input Variable(s)
1363 : real( kind = core_rknd ), intent(in) :: &
1364 : thlm, & ! Liquid Water Potential Temperature [K]
1365 : rtm, & ! Total Water Mixing Ratio [kg/kg]
1366 : p_in_Pa, & ! Pressure [Pa]
1367 : exner ! Exner function [-]
1368 :
1369 : integer, intent(in) :: &
1370 : saturation_formula
1371 :
1372 : ! Output Variable(s)
1373 : real( kind = core_rknd ) :: rcm ! Cloud water mixing ratio [kg/kg]
1374 :
1375 : ! Local Variable(s)
1376 : real( kind = core_rknd ) :: &
1377 : theta, answer, too_low, too_high, & ! [K]
1378 : rsat
1379 :
1380 : integer :: iteration
1381 :
1382 : ! ----- Begin Code -----
1383 :
1384 : ! Default initialization
1385 0 : theta = thlm
1386 0 : too_high = 0.0_core_rknd
1387 0 : too_low = 0.0_core_rknd
1388 :
1389 0 : do iteration = 1, itermax, 1
1390 :
1391 : answer = theta - (Lv/(Cp*exner)) &
1392 0 : *(MAX( rtm - sat_mixrat_liq(p_in_Pa,theta*exner, saturation_formula), zero_threshold ))
1393 :
1394 0 : if ( ABS(answer - thlm) <= tolerance ) then
1395 : exit
1396 0 : else if ( answer - thlm > tolerance ) then
1397 : too_high = theta
1398 0 : else if ( thlm - answer > tolerance ) THEN
1399 0 : too_low = theta
1400 : end if
1401 :
1402 : ! For the first timestep, be sure to set a "too_high"
1403 : ! that is "way too high."
1404 0 : if ( iteration == 1 ) then
1405 0 : too_high = theta + 20.0_core_rknd
1406 : end if
1407 :
1408 0 : theta = (too_low + too_high)/2.0_core_rknd
1409 :
1410 : end do ! 1..itermax
1411 :
1412 0 : if ( iteration == itermax ) then
1413 : ! Magic Eric Raut added to remove compiler warning (clearly this value is not used)
1414 0 : rcm = 0.0_core_rknd
1415 :
1416 0 : error stop "Error in rcm_sat_adj: could not determine rcm"
1417 : else
1418 0 : rcm = MAX( rtm - sat_mixrat_liq( p_in_Pa, theta*exner, saturation_formula), zero_threshold )
1419 : return
1420 : end if
1421 :
1422 : end function rcm_sat_adj
1423 :
1424 : end module saturation
|
