module_gfs_funcphys.F
References to this file elsewhere.
1 !-------------------------------------------------------------------------------
2 module module_gfs_funcphys
3 !$$$ Module Documentation Block
4 !
5 ! Module: funcphys API for basic thermodynamic physics
6 ! Author: Iredell Org: W/NX23 Date: 1999-03-01
7 !
8 ! Abstract: This module provides an Application Program Interface
9 ! for computing basic thermodynamic physics functions, in particular
10 ! (1) saturation vapor pressure as a function of temperature,
11 ! (2) dewpoint temperature as a function of vapor pressure,
12 ! (3) equivalent potential temperature as a function of temperature
13 ! and scaled pressure to the kappa power,
14 ! (4) temperature and specific humidity along a moist adiabat
15 ! as functions of equivalent potential temperature and
16 ! scaled pressure to the kappa power,
17 ! (5) scaled pressure to the kappa power as a function of pressure, and
18 ! (6) temperature at the lifting condensation level as a function
19 ! of temperature and dewpoint depression.
20 ! The entry points required to set up lookup tables start with a "g".
21 ! All the other entry points are functions starting with an "f" or
22 ! are subroutines starting with an "s". These other functions and
23 ! subroutines are elemental; that is, they return a scalar if they
24 ! are passed only scalars, but they return an array if they are passed
25 ! an array. These other functions and subroutines can be inlined, too.
26 !
27 ! Program History Log:
28 ! 1999-03-01 Mark Iredell
29 ! 1999-10-15 Mark Iredell SI unit for pressure (Pascals)
30 ! 2001-02-26 Mark Iredell Ice phase changes of Hong and Moorthi
31 !
32 ! Public Variables:
33 ! krealfp Integer parameter kind or length of reals (=kind_phys)
34 !
35 ! Public Subprograms:
36 ! gpvsl Compute saturation vapor pressure over liquid table
37 !
38 ! fpvsl Elementally compute saturation vapor pressure over liquid
39 ! function result Real(krealfp) saturation vapor pressure in Pascals
40 ! t Real(krealfp) temperature in Kelvin
41 !
42 ! fpvslq Elementally compute saturation vapor pressure over liquid
43 ! function result Real(krealfp) saturation vapor pressure in Pascals
44 ! t Real(krealfp) temperature in Kelvin
45 !
46 ! fpvslx Elementally compute saturation vapor pressure over liquid
47 ! function result Real(krealfp) saturation vapor pressure in Pascals
48 ! t Real(krealfp) temperature in Kelvin
49 !
50 ! gpvsi Compute saturation vapor pressure over ice table
51 !
52 ! fpvsi Elementally compute saturation vapor pressure over ice
53 ! function result Real(krealfp) saturation vapor pressure in Pascals
54 ! t Real(krealfp) temperature in Kelvin
55 !
56 ! fpvsiq Elementally compute saturation vapor pressure over ice
57 ! function result Real(krealfp) saturation vapor pressure in Pascals
58 ! t Real(krealfp) temperature in Kelvin
59 !
60 ! fpvsix Elementally compute saturation vapor pressure over ice
61 ! function result Real(krealfp) saturation vapor pressure in Pascals
62 ! t Real(krealfp) temperature in Kelvin
63 !
64 ! gpvs Compute saturation vapor pressure table
65 !
66 ! fpvs Elementally compute saturation vapor pressure
67 ! function result Real(krealfp) saturation vapor pressure in Pascals
68 ! t Real(krealfp) temperature in Kelvin
69 !
70 ! fpvsq Elementally compute saturation vapor pressure
71 ! function result Real(krealfp) saturation vapor pressure in Pascals
72 ! t Real(krealfp) temperature in Kelvin
73 !
74 ! fpvsx Elementally compute saturation vapor pressure
75 ! function result Real(krealfp) saturation vapor pressure in Pascals
76 ! t Real(krealfp) temperature in Kelvin
77 !
78 ! gtdpl Compute dewpoint temperature over liquid table
79 !
80 ! ftdpl Elementally compute dewpoint temperature over liquid
81 ! function result Real(krealfp) dewpoint temperature in Kelvin
82 ! pv Real(krealfp) vapor pressure in Pascals
83 !
84 ! ftdplq Elementally compute dewpoint temperature over liquid
85 ! function result Real(krealfp) dewpoint temperature in Kelvin
86 ! pv Real(krealfp) vapor pressure in Pascals
87 !
88 ! ftdplx Elementally compute dewpoint temperature over liquid
89 ! function result Real(krealfp) dewpoint temperature in Kelvin
90 ! pv Real(krealfp) vapor pressure in Pascals
91 !
92 ! ftdplxg Elementally compute dewpoint temperature over liquid
93 ! function result Real(krealfp) dewpoint temperature in Kelvin
94 ! t Real(krealfp) guess dewpoint temperature in Kelvin
95 ! pv Real(krealfp) vapor pressure in Pascals
96 !
97 ! gtdpi Compute dewpoint temperature table over ice
98 !
99 ! ftdpi Elementally compute dewpoint temperature over ice
100 ! function result Real(krealfp) dewpoint temperature in Kelvin
101 ! pv Real(krealfp) vapor pressure in Pascals
102 !
103 ! ftdpiq Elementally compute dewpoint temperature over ice
104 ! function result Real(krealfp) dewpoint temperature in Kelvin
105 ! pv Real(krealfp) vapor pressure in Pascals
106 !
107 ! ftdpix Elementally compute dewpoint temperature over ice
108 ! function result Real(krealfp) dewpoint temperature in Kelvin
109 ! pv Real(krealfp) vapor pressure in Pascals
110 !
111 ! ftdpixg Elementally compute dewpoint temperature over ice
112 ! function result Real(krealfp) dewpoint temperature in Kelvin
113 ! t Real(krealfp) guess dewpoint temperature in Kelvin
114 ! pv Real(krealfp) vapor pressure in Pascals
115 !
116 ! gtdp Compute dewpoint temperature table
117 !
118 ! ftdp Elementally compute dewpoint temperature
119 ! function result Real(krealfp) dewpoint temperature in Kelvin
120 ! pv Real(krealfp) vapor pressure in Pascals
121 !
122 ! ftdpq Elementally compute dewpoint temperature
123 ! function result Real(krealfp) dewpoint temperature in Kelvin
124 ! pv Real(krealfp) vapor pressure in Pascals
125 !
126 ! ftdpx Elementally compute dewpoint temperature
127 ! function result Real(krealfp) dewpoint temperature in Kelvin
128 ! pv Real(krealfp) vapor pressure in Pascals
129 !
130 ! ftdpxg Elementally compute dewpoint temperature
131 ! function result Real(krealfp) dewpoint temperature in Kelvin
132 ! t Real(krealfp) guess dewpoint temperature in Kelvin
133 ! pv Real(krealfp) vapor pressure in Pascals
134 !
135 ! gthe Compute equivalent potential temperature table
136 !
137 ! fthe Elementally compute equivalent potential temperature
138 ! function result Real(krealfp) equivalent potential temperature in Kelvin
139 ! t Real(krealfp) LCL temperature in Kelvin
140 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
141 !
142 ! ftheq Elementally compute equivalent potential temperature
143 ! function result Real(krealfp) equivalent potential temperature in Kelvin
144 ! t Real(krealfp) LCL temperature in Kelvin
145 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
146 !
147 ! fthex Elementally compute equivalent potential temperature
148 ! function result Real(krealfp) equivalent potential temperature in Kelvin
149 ! t Real(krealfp) LCL temperature in Kelvin
150 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
151 !
152 ! gtma Compute moist adiabat tables
153 !
154 ! stma Elementally compute moist adiabat temperature and moisture
155 ! the Real(krealfp) equivalent potential temperature in Kelvin
156 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
157 ! tma Real(krealfp) parcel temperature in Kelvin
158 ! qma Real(krealfp) parcel specific humidity in kg/kg
159 !
160 ! stmaq Elementally compute moist adiabat temperature and moisture
161 ! the Real(krealfp) equivalent potential temperature in Kelvin
162 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
163 ! tma Real(krealfp) parcel temperature in Kelvin
164 ! qma Real(krealfp) parcel specific humidity in kg/kg
165 !
166 ! stmax Elementally compute moist adiabat temperature and moisture
167 ! the Real(krealfp) equivalent potential temperature in Kelvin
168 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
169 ! tma Real(krealfp) parcel temperature in Kelvin
170 ! qma Real(krealfp) parcel specific humidity in kg/kg
171 !
172 ! stmaxg Elementally compute moist adiabat temperature and moisture
173 ! tg Real(krealfp) guess parcel temperature in Kelvin
174 ! the Real(krealfp) equivalent potential temperature in Kelvin
175 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
176 ! tma Real(krealfp) parcel temperature in Kelvin
177 ! qma Real(krealfp) parcel specific humidity in kg/kg
178 !
179 ! gpkap Compute pressure to the kappa table
180 !
181 ! fpkap Elementally raise pressure to the kappa power.
182 ! function result Real(krealfp) p over 1e5 Pa to the kappa power
183 ! p Real(krealfp) pressure in Pascals
184 !
185 ! fpkapq Elementally raise pressure to the kappa power.
186 ! function result Real(krealfp) p over 1e5 Pa to the kappa power
187 ! p Real(krealfp) pressure in Pascals
188 !
189 ! fpkapo Elementally raise pressure to the kappa power.
190 ! function result Real(krealfp) p over 1e5 Pa to the kappa power
191 ! p Real(krealfp) surface pressure in Pascals
192 !
193 ! fpkapx Elementally raise pressure to the kappa power.
194 ! function result Real(krealfp) p over 1e5 Pa to the kappa power
195 ! p Real(krealfp) pressure in Pascals
196 !
197 ! grkap Compute pressure to the 1/kappa table
198 !
199 ! frkap Elementally raise pressure to the 1/kappa power.
200 ! function result Real(krealfp) pressure in Pascals
201 ! pkap Real(krealfp) p over 1e5 Pa to the 1/kappa power
202 !
203 ! frkapq Elementally raise pressure to the kappa power.
204 ! function result Real(krealfp) pressure in Pascals
205 ! pkap Real(krealfp) p over 1e5 Pa to the kappa power
206 !
207 ! frkapx Elementally raise pressure to the kappa power.
208 ! function result Real(krealfp) pressure in Pascals
209 ! pkap Real(krealfp) p over 1e5 Pa to the kappa power
210 !
211 ! gtlcl Compute LCL temperature table
212 !
213 ! ftlcl Elementally compute LCL temperature.
214 ! function result Real(krealfp) temperature at the LCL in Kelvin
215 ! t Real(krealfp) temperature in Kelvin
216 ! tdpd Real(krealfp) dewpoint depression in Kelvin
217 !
218 ! ftlclq Elementally compute LCL temperature.
219 ! function result Real(krealfp) temperature at the LCL in Kelvin
220 ! t Real(krealfp) temperature in Kelvin
221 ! tdpd Real(krealfp) dewpoint depression in Kelvin
222 !
223 ! ftlclo Elementally compute LCL temperature.
224 ! function result Real(krealfp) temperature at the LCL in Kelvin
225 ! t Real(krealfp) temperature in Kelvin
226 ! tdpd Real(krealfp) dewpoint depression in Kelvin
227 !
228 ! ftlclx Elementally compute LCL temperature.
229 ! function result Real(krealfp) temperature at the LCL in Kelvin
230 ! t Real(krealfp) temperature in Kelvin
231 ! tdpd Real(krealfp) dewpoint depression in Kelvin
232 !
233 ! gfuncphys Compute all physics function tables
234 !
235 ! Attributes:
236 ! Language: Fortran 90
237 !
238 !$$$
239 use module_gfs_machine,only:kind_phys
240 use module_gfs_physcons
241 implicit none
242 private
243 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244 ! Public Variables
245 ! integer,public,parameter:: krealfp=selected_real_kind(15,45)
246 integer,public,parameter:: krealfp=kind_phys
247 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248 ! Private Variables
249 real(krealfp),parameter:: psatb=con_psat*1.e-5
250 integer,parameter:: nxpvsl=7501
251 real(krealfp) c1xpvsl,c2xpvsl,tbpvsl(nxpvsl)
252 integer,parameter:: nxpvsi=7501
253 real(krealfp) c1xpvsi,c2xpvsi,tbpvsi(nxpvsi)
254 integer,parameter:: nxpvs=7501
255 real(krealfp) c1xpvs,c2xpvs,tbpvs(nxpvs)
256 integer,parameter:: nxtdpl=5001
257 real(krealfp) c1xtdpl,c2xtdpl,tbtdpl(nxtdpl)
258 integer,parameter:: nxtdpi=5001
259 real(krealfp) c1xtdpi,c2xtdpi,tbtdpi(nxtdpi)
260 integer,parameter:: nxtdp=5001
261 real(krealfp) c1xtdp,c2xtdp,tbtdp(nxtdp)
262 integer,parameter:: nxthe=241,nythe=151
263 real(krealfp) c1xthe,c2xthe,c1ythe,c2ythe,tbthe(nxthe,nythe)
264 integer,parameter:: nxma=151,nyma=121
265 real(krealfp) c1xma,c2xma,c1yma,c2yma,tbtma(nxma,nyma),tbqma(nxma,nyma)
266 ! integer,parameter:: nxpkap=5501
267 integer,parameter:: nxpkap=11001
268 real(krealfp) c1xpkap,c2xpkap,tbpkap(nxpkap)
269 integer,parameter:: nxrkap=5501
270 real(krealfp) c1xrkap,c2xrkap,tbrkap(nxrkap)
271 integer,parameter:: nxtlcl=151,nytlcl=61
272 real(krealfp) c1xtlcl,c2xtlcl,c1ytlcl,c2ytlcl,tbtlcl(nxtlcl,nytlcl)
273 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274 ! Public Subprograms
275 public gpvsl,fpvsl,fpvslq,fpvslx
276 public gpvsi,fpvsi,fpvsiq,fpvsix
277 public gpvs,fpvs,fpvsq,fpvsx
278 public gtdpl,ftdpl,ftdplq,ftdplx,ftdplxg
279 public gtdpi,ftdpi,ftdpiq,ftdpix,ftdpixg
280 public gtdp,ftdp,ftdpq,ftdpx,ftdpxg
281 public gthe,fthe,ftheq,fthex
282 public gtma,stma,stmaq,stmax,stmaxg
283 public gpkap,fpkap,fpkapq,fpkapo,fpkapx
284 public grkap,frkap,frkapq,frkapx
285 public gtlcl,ftlcl,ftlclq,ftlclo,ftlclx
286 public gfuncphys
287 contains
288 !-------------------------------------------------------------------------------
289 subroutine gpvsl
290 !$$$ Subprogram Documentation Block
291 !
292 ! Subprogram: gpvsl Compute saturation vapor pressure table over liquid
293 ! Author: N Phillips W/NMC2X2 Date: 30 dec 82
294 !
295 ! Abstract: Computes saturation vapor pressure table as a function of
296 ! temperature for the table lookup function fpvsl.
297 ! Exact saturation vapor pressures are calculated in subprogram fpvslx.
298 ! The current implementation computes a table with a length
299 ! of 7501 for temperatures ranging from 180. to 330. Kelvin.
300 !
301 ! Program History Log:
302 ! 91-05-07 Iredell
303 ! 94-12-30 Iredell expand table
304 ! 1999-03-01 Iredell f90 module
305 !
306 ! Usage: call gpvsl
307 !
308 ! Subprograms called:
309 ! (fpvslx) inlinable function to compute saturation vapor pressure
310 !
311 ! Attributes:
312 ! Language: Fortran 90.
313 !
314 !$$$
315 implicit none
316 integer jx
317 real(krealfp) xmin,xmax,xinc,x,t
318 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
319 xmin=180.0_krealfp
320 xmax=330.0_krealfp
321 xinc=(xmax-xmin)/(nxpvsl-1)
322 ! c1xpvsl=1.-xmin/xinc
323 c2xpvsl=1./xinc
324 c1xpvsl=1.-xmin*c2xpvsl
325 do jx=1,nxpvsl
326 x=xmin+(jx-1)*xinc
327 t=x
328 tbpvsl(jx)=fpvslx(t)
329 enddo
330 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
331 end subroutine
332 !-------------------------------------------------------------------------------
333 ! elemental function fpvsl(t)
334 function fpvsl(t)
335 !$$$ Subprogram Documentation Block
336 !
337 ! Subprogram: fpvsl Compute saturation vapor pressure over liquid
338 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
339 !
340 ! Abstract: Compute saturation vapor pressure from the temperature.
341 ! A linear interpolation is done between values in a lookup table
342 ! computed in gpvsl. See documentation for fpvslx for details.
343 ! Input values outside table range are reset to table extrema.
344 ! The interpolation accuracy is almost 6 decimal places.
345 ! On the Cray, fpvsl is about 4 times faster than exact calculation.
346 ! This function should be expanded inline in the calling routine.
347 !
348 ! Program History Log:
349 ! 91-05-07 Iredell made into inlinable function
350 ! 94-12-30 Iredell expand table
351 ! 1999-03-01 Iredell f90 module
352 !
353 ! Usage: pvsl=fpvsl(t)
354 !
355 ! Input argument list:
356 ! t Real(krealfp) temperature in Kelvin
357 !
358 ! Output argument list:
359 ! fpvsl Real(krealfp) saturation vapor pressure in Pascals
360 !
361 ! Attributes:
362 ! Language: Fortran 90.
363 !
364 !$$$
365 implicit none
366 real(krealfp) fpvsl
367 real(krealfp),intent(in):: t
368 integer jx
369 real(krealfp) xj
370 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371 xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp))
372 jx=min(xj,nxpvsl-1._krealfp)
373 fpvsl=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx))
374 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
375 end function
376 !-------------------------------------------------------------------------------
377 ! elemental function fpvslq(t)
378 function fpvslq(t)
379 !$$$ Subprogram Documentation Block
380 !
381 ! Subprogram: fpvslq Compute saturation vapor pressure over liquid
382 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
383 !
384 ! Abstract: Compute saturation vapor pressure from the temperature.
385 ! A quadratic interpolation is done between values in a lookup table
386 ! computed in gpvsl. See documentation for fpvslx for details.
387 ! Input values outside table range are reset to table extrema.
388 ! The interpolation accuracy is almost 9 decimal places.
389 ! On the Cray, fpvslq is about 3 times faster than exact calculation.
390 ! This function should be expanded inline in the calling routine.
391 !
392 ! Program History Log:
393 ! 91-05-07 Iredell made into inlinable function
394 ! 94-12-30 Iredell quadratic interpolation
395 ! 1999-03-01 Iredell f90 module
396 !
397 ! Usage: pvsl=fpvslq(t)
398 !
399 ! Input argument list:
400 ! t Real(krealfp) temperature in Kelvin
401 !
402 ! Output argument list:
403 ! fpvslq Real(krealfp) saturation vapor pressure in Pascals
404 !
405 ! Attributes:
406 ! Language: Fortran 90.
407 !
408 !$$$
409 implicit none
410 real(krealfp) fpvslq
411 real(krealfp),intent(in):: t
412 integer jx
413 real(krealfp) xj,dxj,fj1,fj2,fj3
414 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
415 xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp))
416 jx=min(max(nint(xj),2),nxpvsl-1)
417 dxj=xj-jx
418 fj1=tbpvsl(jx-1)
419 fj2=tbpvsl(jx)
420 fj3=tbpvsl(jx+1)
421 fpvslq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
422 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
423 end function
424 !-------------------------------------------------------------------------------
425 ! elemental function fpvslx(t)
426 function fpvslx(t)
427 !$$$ Subprogram Documentation Block
428 !
429 ! Subprogram: fpvslx Compute saturation vapor pressure over liquid
430 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
431 !
432 ! Abstract: Exactly compute saturation vapor pressure from temperature.
433 ! The water model assumes a perfect gas, constant specific heats
434 ! for gas and liquid, and neglects the volume of the liquid.
435 ! The model does account for the variation of the latent heat
436 ! of condensation with temperature. The ice option is not included.
437 ! The Clausius-Clapeyron equation is integrated from the triple point
438 ! to get the formula
439 ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
440 ! where tr is ttp/t and other values are physical constants.
441 ! This function should be expanded inline in the calling routine.
442 !
443 ! Program History Log:
444 ! 91-05-07 Iredell made into inlinable function
445 ! 94-12-30 Iredell exact computation
446 ! 1999-03-01 Iredell f90 module
447 !
448 ! Usage: pvsl=fpvslx(t)
449 !
450 ! Input argument list:
451 ! t Real(krealfp) temperature in Kelvin
452 !
453 ! Output argument list:
454 ! fpvslx Real(krealfp) saturation vapor pressure in Pascals
455 !
456 ! Attributes:
457 ! Language: Fortran 90.
458 !
459 !$$$
460 implicit none
461 real(krealfp) fpvslx
462 real(krealfp),intent(in):: t
463 real(krealfp),parameter:: dldt=con_cvap-con_cliq
464 real(krealfp),parameter:: heat=con_hvap
465 real(krealfp),parameter:: xpona=-dldt/con_rv
466 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
467 real(krealfp) tr
468 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469 tr=con_ttp/t
470 fpvslx=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
471 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
472 end function
473 !-------------------------------------------------------------------------------
474 subroutine gpvsi
475 !$$$ Subprogram Documentation Block
476 !
477 ! Subprogram: gpvsi Compute saturation vapor pressure table over ice
478 ! Author: N Phillips W/NMC2X2 Date: 30 dec 82
479 !
480 ! Abstract: Computes saturation vapor pressure table as a function of
481 ! temperature for the table lookup function fpvsi.
482 ! Exact saturation vapor pressures are calculated in subprogram fpvsix.
483 ! The current implementation computes a table with a length
484 ! of 7501 for temperatures ranging from 180. to 330. Kelvin.
485 !
486 ! Program History Log:
487 ! 91-05-07 Iredell
488 ! 94-12-30 Iredell expand table
489 ! 1999-03-01 Iredell f90 module
490 ! 2001-02-26 Iredell ice phase
491 !
492 ! Usage: call gpvsi
493 !
494 ! Subprograms called:
495 ! (fpvsix) inlinable function to compute saturation vapor pressure
496 !
497 ! Attributes:
498 ! Language: Fortran 90.
499 !
500 !$$$
501 implicit none
502 integer jx
503 real(krealfp) xmin,xmax,xinc,x,t
504 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
505 xmin=180.0_krealfp
506 xmax=330.0_krealfp
507 xinc=(xmax-xmin)/(nxpvsi-1)
508 ! c1xpvsi=1.-xmin/xinc
509 c2xpvsi=1./xinc
510 c1xpvsi=1.-xmin*c2xpvsi
511 do jx=1,nxpvsi
512 x=xmin+(jx-1)*xinc
513 t=x
514 tbpvsi(jx)=fpvsix(t)
515 enddo
516 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
517 end subroutine
518 !-------------------------------------------------------------------------------
519 ! elemental function fpvsi(t)
520 function fpvsi(t)
521 !$$$ Subprogram Documentation Block
522 !
523 ! Subprogram: fpvsi Compute saturation vapor pressure over ice
524 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
525 !
526 ! Abstract: Compute saturation vapor pressure from the temperature.
527 ! A linear interpolation is done between values in a lookup table
528 ! computed in gpvsi. See documentation for fpvsix for details.
529 ! Input values outside table range are reset to table extrema.
530 ! The interpolation accuracy is almost 6 decimal places.
531 ! On the Cray, fpvsi is about 4 times faster than exact calculation.
532 ! This function should be expanded inline in the calling routine.
533 !
534 ! Program History Log:
535 ! 91-05-07 Iredell made into inlinable function
536 ! 94-12-30 Iredell expand table
537 ! 1999-03-01 Iredell f90 module
538 ! 2001-02-26 Iredell ice phase
539 !
540 ! Usage: pvsi=fpvsi(t)
541 !
542 ! Input argument list:
543 ! t Real(krealfp) temperature in Kelvin
544 !
545 ! Output argument list:
546 ! fpvsi Real(krealfp) saturation vapor pressure in Pascals
547 !
548 ! Attributes:
549 ! Language: Fortran 90.
550 !
551 !$$$
552 implicit none
553 real(krealfp) fpvsi
554 real(krealfp),intent(in):: t
555 integer jx
556 real(krealfp) xj
557 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
558 xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp))
559 jx=min(xj,nxpvsi-1._krealfp)
560 fpvsi=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx))
561 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
562 end function
563 !-------------------------------------------------------------------------------
564 ! elemental function fpvsiq(t)
565 function fpvsiq(t)
566 !$$$ Subprogram Documentation Block
567 !
568 ! Subprogram: fpvsiq Compute saturation vapor pressure over ice
569 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
570 !
571 ! Abstract: Compute saturation vapor pressure from the temperature.
572 ! A quadratic interpolation is done between values in a lookup table
573 ! computed in gpvsi. See documentation for fpvsix for details.
574 ! Input values outside table range are reset to table extrema.
575 ! The interpolation accuracy is almost 9 decimal places.
576 ! On the Cray, fpvsiq is about 3 times faster than exact calculation.
577 ! This function should be expanded inline in the calling routine.
578 !
579 ! Program History Log:
580 ! 91-05-07 Iredell made into inlinable function
581 ! 94-12-30 Iredell quadratic interpolation
582 ! 1999-03-01 Iredell f90 module
583 ! 2001-02-26 Iredell ice phase
584 !
585 ! Usage: pvsi=fpvsiq(t)
586 !
587 ! Input argument list:
588 ! t Real(krealfp) temperature in Kelvin
589 !
590 ! Output argument list:
591 ! fpvsiq Real(krealfp) saturation vapor pressure in Pascals
592 !
593 ! Attributes:
594 ! Language: Fortran 90.
595 !
596 !$$$
597 implicit none
598 real(krealfp) fpvsiq
599 real(krealfp),intent(in):: t
600 integer jx
601 real(krealfp) xj,dxj,fj1,fj2,fj3
602 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
603 xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp))
604 jx=min(max(nint(xj),2),nxpvsi-1)
605 dxj=xj-jx
606 fj1=tbpvsi(jx-1)
607 fj2=tbpvsi(jx)
608 fj3=tbpvsi(jx+1)
609 fpvsiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
610 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
611 end function
612 !-------------------------------------------------------------------------------
613 ! elemental function fpvsix(t)
614 function fpvsix(t)
615 !$$$ Subprogram Documentation Block
616 !
617 ! Subprogram: fpvsix Compute saturation vapor pressure over ice
618 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
619 !
620 ! Abstract: Exactly compute saturation vapor pressure from temperature.
621 ! The water model assumes a perfect gas, constant specific heats
622 ! for gas and ice, and neglects the volume of the ice.
623 ! The model does account for the variation of the latent heat
624 ! of condensation with temperature. The liquid option is not included.
625 ! The Clausius-Clapeyron equation is integrated from the triple point
626 ! to get the formula
627 ! pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr))
628 ! where tr is ttp/t and other values are physical constants.
629 ! This function should be expanded inline in the calling routine.
630 !
631 ! Program History Log:
632 ! 91-05-07 Iredell made into inlinable function
633 ! 94-12-30 Iredell exact computation
634 ! 1999-03-01 Iredell f90 module
635 ! 2001-02-26 Iredell ice phase
636 !
637 ! Usage: pvsi=fpvsix(t)
638 !
639 ! Input argument list:
640 ! t Real(krealfp) temperature in Kelvin
641 !
642 ! Output argument list:
643 ! fpvsix Real(krealfp) saturation vapor pressure in Pascals
644 !
645 ! Attributes:
646 ! Language: Fortran 90.
647 !
648 !$$$
649 implicit none
650 real(krealfp) fpvsix
651 real(krealfp),intent(in):: t
652 real(krealfp),parameter:: dldt=con_cvap-con_csol
653 real(krealfp),parameter:: heat=con_hvap+con_hfus
654 real(krealfp),parameter:: xpona=-dldt/con_rv
655 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
656 real(krealfp) tr
657 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
658 tr=con_ttp/t
659 fpvsix=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
660 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
661 end function
662 !-------------------------------------------------------------------------------
663 subroutine gpvs
664 !$$$ Subprogram Documentation Block
665 !
666 ! Subprogram: gpvs Compute saturation vapor pressure table
667 ! Author: N Phillips W/NMC2X2 Date: 30 dec 82
668 !
669 ! Abstract: Computes saturation vapor pressure table as a function of
670 ! temperature for the table lookup function fpvs.
671 ! Exact saturation vapor pressures are calculated in subprogram fpvsx.
672 ! The current implementation computes a table with a length
673 ! of 7501 for temperatures ranging from 180. to 330. Kelvin.
674 !
675 ! Program History Log:
676 ! 91-05-07 Iredell
677 ! 94-12-30 Iredell expand table
678 ! 1999-03-01 Iredell f90 module
679 ! 2001-02-26 Iredell ice phase
680 !
681 ! Usage: call gpvs
682 !
683 ! Subprograms called:
684 ! (fpvsx) inlinable function to compute saturation vapor pressure
685 !
686 ! Attributes:
687 ! Language: Fortran 90.
688 !
689 !$$$
690 implicit none
691 integer jx
692 real(krealfp) xmin,xmax,xinc,x,t
693 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
694 xmin=180.0_krealfp
695 xmax=330.0_krealfp
696 xinc=(xmax-xmin)/(nxpvs-1)
697 ! c1xpvs=1.-xmin/xinc
698 c2xpvs=1./xinc
699 c1xpvs=1.-xmin*c2xpvs
700 do jx=1,nxpvs
701 x=xmin+(jx-1)*xinc
702 t=x
703 tbpvs(jx)=fpvsx(t)
704 enddo
705 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
706 end subroutine
707 !-------------------------------------------------------------------------------
708 ! elemental function fpvs(t)
709 function fpvs(t)
710 !$$$ Subprogram Documentation Block
711 !
712 ! Subprogram: fpvs Compute saturation vapor pressure
713 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
714 !
715 ! Abstract: Compute saturation vapor pressure from the temperature.
716 ! A linear interpolation is done between values in a lookup table
717 ! computed in gpvs. See documentation for fpvsx for details.
718 ! Input values outside table range are reset to table extrema.
719 ! The interpolation accuracy is almost 6 decimal places.
720 ! On the Cray, fpvs is about 4 times faster than exact calculation.
721 ! This function should be expanded inline in the calling routine.
722 !
723 ! Program History Log:
724 ! 91-05-07 Iredell made into inlinable function
725 ! 94-12-30 Iredell expand table
726 ! 1999-03-01 Iredell f90 module
727 ! 2001-02-26 Iredell ice phase
728 !
729 ! Usage: pvs=fpvs(t)
730 !
731 ! Input argument list:
732 ! t Real(krealfp) temperature in Kelvin
733 !
734 ! Output argument list:
735 ! fpvs Real(krealfp) saturation vapor pressure in Pascals
736 !
737 ! Attributes:
738 ! Language: Fortran 90.
739 !
740 !$$$
741 implicit none
742 real(krealfp) fpvs
743 real(krealfp),intent(in):: t
744 integer jx
745 real(krealfp) xj
746 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
747 xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
748 jx=min(xj,nxpvs-1._krealfp)
749 fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
750 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
751 end function
752 !-------------------------------------------------------------------------------
753 ! elemental function fpvsq(t)
754 function fpvsq(t)
755 !$$$ Subprogram Documentation Block
756 !
757 ! Subprogram: fpvsq Compute saturation vapor pressure
758 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
759 !
760 ! Abstract: Compute saturation vapor pressure from the temperature.
761 ! A quadratic interpolation is done between values in a lookup table
762 ! computed in gpvs. See documentation for fpvsx for details.
763 ! Input values outside table range are reset to table extrema.
764 ! The interpolation accuracy is almost 9 decimal places.
765 ! On the Cray, fpvsq is about 3 times faster than exact calculation.
766 ! This function should be expanded inline in the calling routine.
767 !
768 ! Program History Log:
769 ! 91-05-07 Iredell made into inlinable function
770 ! 94-12-30 Iredell quadratic interpolation
771 ! 1999-03-01 Iredell f90 module
772 ! 2001-02-26 Iredell ice phase
773 !
774 ! Usage: pvs=fpvsq(t)
775 !
776 ! Input argument list:
777 ! t Real(krealfp) temperature in Kelvin
778 !
779 ! Output argument list:
780 ! fpvsq Real(krealfp) saturation vapor pressure in Pascals
781 !
782 ! Attributes:
783 ! Language: Fortran 90.
784 !
785 !$$$
786 implicit none
787 real(krealfp) fpvsq
788 real(krealfp),intent(in):: t
789 integer jx
790 real(krealfp) xj,dxj,fj1,fj2,fj3
791 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
792 xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
793 jx=min(max(nint(xj),2),nxpvs-1)
794 dxj=xj-jx
795 fj1=tbpvs(jx-1)
796 fj2=tbpvs(jx)
797 fj3=tbpvs(jx+1)
798 fpvsq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
799 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
800 end function
801 !-------------------------------------------------------------------------------
802 ! elemental function fpvsx(t)
803 function fpvsx(t)
804 !$$$ Subprogram Documentation Block
805 !
806 ! Subprogram: fpvsx Compute saturation vapor pressure
807 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
808 !
809 ! Abstract: Exactly compute saturation vapor pressure from temperature.
810 ! The saturation vapor pressure over either liquid and ice is computed
811 ! over liquid for temperatures above the triple point,
812 ! over ice for temperatures 20 degress below the triple point,
813 ! and a linear combination of the two for temperatures in between.
814 ! The water model assumes a perfect gas, constant specific heats
815 ! for gas, liquid and ice, and neglects the volume of the condensate.
816 ! The model does account for the variation of the latent heat
817 ! of condensation and sublimation with temperature.
818 ! The Clausius-Clapeyron equation is integrated from the triple point
819 ! to get the formula
820 ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
821 ! where tr is ttp/t and other values are physical constants.
822 ! The reference for this computation is Emanuel(1994), pages 116-117.
823 ! This function should be expanded inline in the calling routine.
824 !
825 ! Program History Log:
826 ! 91-05-07 Iredell made into inlinable function
827 ! 94-12-30 Iredell exact computation
828 ! 1999-03-01 Iredell f90 module
829 ! 2001-02-26 Iredell ice phase
830 !
831 ! Usage: pvs=fpvsx(t)
832 !
833 ! Input argument list:
834 ! t Real(krealfp) temperature in Kelvin
835 !
836 ! Output argument list:
837 ! fpvsx Real(krealfp) saturation vapor pressure in Pascals
838 !
839 ! Attributes:
840 ! Language: Fortran 90.
841 !
842 !$$$
843 implicit none
844 real(krealfp) fpvsx
845 real(krealfp),intent(in):: t
846 real(krealfp),parameter:: tliq=con_ttp
847 real(krealfp),parameter:: tice=con_ttp-20.0
848 real(krealfp),parameter:: dldtl=con_cvap-con_cliq
849 real(krealfp),parameter:: heatl=con_hvap
850 real(krealfp),parameter:: xponal=-dldtl/con_rv
851 real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
852 real(krealfp),parameter:: dldti=con_cvap-con_csol
853 real(krealfp),parameter:: heati=con_hvap+con_hfus
854 real(krealfp),parameter:: xponai=-dldti/con_rv
855 real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
856 real(krealfp) tr,w,pvl,pvi
857 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
858 tr=con_ttp/t
859 if(t.ge.tliq) then
860 fpvsx=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
861 elseif(t.lt.tice) then
862 fpvsx=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
863 else
864 w=(t-tice)/(tliq-tice)
865 pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
866 pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
867 fpvsx=w*pvl+(1.-w)*pvi
868 endif
869 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
870 end function
871 !-------------------------------------------------------------------------------
872 subroutine gtdpl
873 !$$$ Subprogram Documentation Block
874 !
875 ! Subprogram: gtdpl Compute dewpoint temperature over liquid table
876 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
877 !
878 ! Abstract: Compute dewpoint temperature table as a function of
879 ! vapor pressure for inlinable function ftdpl.
880 ! Exact dewpoint temperatures are calculated in subprogram ftdplxg.
881 ! The current implementation computes a table with a length
882 ! of 5001 for vapor pressures ranging from 1 to 10001 Pascals
883 ! giving a dewpoint temperature range of 208 to 319 Kelvin.
884 !
885 ! Program History Log:
886 ! 91-05-07 Iredell
887 ! 94-12-30 Iredell expand table
888 ! 1999-03-01 Iredell f90 module
889 !
890 ! Usage: call gtdpl
891 !
892 ! Subprograms called:
893 ! (ftdplxg) inlinable function to compute dewpoint temperature over liquid
894 !
895 ! Attributes:
896 ! Language: Fortran 90.
897 !
898 !$$$
899 implicit none
900 integer jx
901 real(krealfp) xmin,xmax,xinc,t,x,pv
902 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
903 xmin=1
904 xmax=10001
905 xinc=(xmax-xmin)/(nxtdpl-1)
906 c1xtdpl=1.-xmin/xinc
907 c2xtdpl=1./xinc
908 t=208.0
909 do jx=1,nxtdpl
910 x=xmin+(jx-1)*xinc
911 pv=x
912 t=ftdplxg(t,pv)
913 tbtdpl(jx)=t
914 enddo
915 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
916 end subroutine
917 !-------------------------------------------------------------------------------
918 ! elemental function ftdpl(pv)
919 function ftdpl(pv)
920 !$$$ Subprogram Documentation Block
921 !
922 ! Subprogram: ftdpl Compute dewpoint temperature over liquid
923 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
924 !
925 ! Abstract: Compute dewpoint temperature from vapor pressure.
926 ! A linear interpolation is done between values in a lookup table
927 ! computed in gtdpl. See documentation for ftdplxg for details.
928 ! Input values outside table range are reset to table extrema.
929 ! The interpolation accuracy is better than 0.0005 Kelvin
930 ! for dewpoint temperatures greater than 250 Kelvin,
931 ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
932 ! On the Cray, ftdpl is about 75 times faster than exact calculation.
933 ! This function should be expanded inline in the calling routine.
934 !
935 ! Program History Log:
936 ! 91-05-07 Iredell made into inlinable function
937 ! 94-12-30 Iredell expand table
938 ! 1999-03-01 Iredell f90 module
939 !
940 ! Usage: tdpl=ftdpl(pv)
941 !
942 ! Input argument list:
943 ! pv Real(krealfp) vapor pressure in Pascals
944 !
945 ! Output argument list:
946 ! ftdpl Real(krealfp) dewpoint temperature in Kelvin
947 !
948 ! Attributes:
949 ! Language: Fortran 90.
950 !
951 !$$$
952 implicit none
953 real(krealfp) ftdpl
954 real(krealfp),intent(in):: pv
955 integer jx
956 real(krealfp) xj
957 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
958 xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
959 jx=min(xj,nxtdpl-1._krealfp)
960 ftdpl=tbtdpl(jx)+(xj-jx)*(tbtdpl(jx+1)-tbtdpl(jx))
961 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
962 end function
963 !-------------------------------------------------------------------------------
964 ! elemental function ftdplq(pv)
965 function ftdplq(pv)
966 !$$$ Subprogram Documentation Block
967 !
968 ! Subprogram: ftdplq Compute dewpoint temperature over liquid
969 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
970 !
971 ! Abstract: Compute dewpoint temperature from vapor pressure.
972 ! A quadratic interpolation is done between values in a lookup table
973 ! computed in gtdpl. see documentation for ftdplxg for details.
974 ! Input values outside table range are reset to table extrema.
975 ! the interpolation accuracy is better than 0.00001 Kelvin
976 ! for dewpoint temperatures greater than 250 Kelvin,
977 ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
978 ! On the Cray, ftdplq is about 60 times faster than exact calculation.
979 ! This function should be expanded inline in the calling routine.
980 !
981 ! Program History Log:
982 ! 91-05-07 Iredell made into inlinable function
983 ! 94-12-30 Iredell quadratic interpolation
984 ! 1999-03-01 Iredell f90 module
985 !
986 ! Usage: tdpl=ftdplq(pv)
987 !
988 ! Input argument list:
989 ! pv Real(krealfp) vapor pressure in Pascals
990 !
991 ! Output argument list:
992 ! ftdplq Real(krealfp) dewpoint temperature in Kelvin
993 !
994 ! Attributes:
995 ! Language: Fortran 90.
996 !
997 !$$$
998 implicit none
999 real(krealfp) ftdplq
1000 real(krealfp),intent(in):: pv
1001 integer jx
1002 real(krealfp) xj,dxj,fj1,fj2,fj3
1003 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1004 xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
1005 jx=min(max(nint(xj),2),nxtdpl-1)
1006 dxj=xj-jx
1007 fj1=tbtdpl(jx-1)
1008 fj2=tbtdpl(jx)
1009 fj3=tbtdpl(jx+1)
1010 ftdplq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1011 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1012 end function
1013 !-------------------------------------------------------------------------------
1014 ! elemental function ftdplx(pv)
1015 function ftdplx(pv)
1016 !$$$ Subprogram Documentation Block
1017 !
1018 ! Subprogram: ftdplx Compute dewpoint temperature over liquid
1019 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1020 !
1021 ! Abstract: exactly compute dewpoint temperature from vapor pressure.
1022 ! An approximate dewpoint temperature for function ftdplxg
1023 ! is obtained using ftdpl so gtdpl must be already called.
1024 ! See documentation for ftdplxg for details.
1025 !
1026 ! Program History Log:
1027 ! 91-05-07 Iredell made into inlinable function
1028 ! 94-12-30 Iredell exact computation
1029 ! 1999-03-01 Iredell f90 module
1030 !
1031 ! Usage: tdpl=ftdplx(pv)
1032 !
1033 ! Input argument list:
1034 ! pv Real(krealfp) vapor pressure in Pascals
1035 !
1036 ! Output argument list:
1037 ! ftdplx Real(krealfp) dewpoint temperature in Kelvin
1038 !
1039 ! Subprograms called:
1040 ! (ftdpl) inlinable function to compute dewpoint temperature over liquid
1041 ! (ftdplxg) inlinable function to compute dewpoint temperature over liquid
1042 !
1043 ! Attributes:
1044 ! Language: Fortran 90.
1045 !
1046 !$$$
1047 implicit none
1048 real(krealfp) ftdplx
1049 real(krealfp),intent(in):: pv
1050 real(krealfp) tg
1051 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1052 tg=ftdpl(pv)
1053 ftdplx=ftdplxg(tg,pv)
1054 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1055 end function
1056 !-------------------------------------------------------------------------------
1057 ! elemental function ftdplxg(tg,pv)
1058 function ftdplxg(tg,pv)
1059 !$$$ Subprogram Documentation Block
1060 !
1061 ! Subprogram: ftdplxg Compute dewpoint temperature over liquid
1062 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1063 !
1064 ! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1065 ! A guess dewpoint temperature must be provided.
1066 ! The water model assumes a perfect gas, constant specific heats
1067 ! for gas and liquid, and neglects the volume of the liquid.
1068 ! The model does account for the variation of the latent heat
1069 ! of condensation with temperature. The ice option is not included.
1070 ! The Clausius-Clapeyron equation is integrated from the triple point
1071 ! to get the formula
1072 ! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
1073 ! where tr is ttp/t and other values are physical constants.
1074 ! The formula is inverted by iterating Newtonian approximations
1075 ! for each pvs until t is found to within 1.e-6 Kelvin.
1076 ! This function can be expanded inline in the calling routine.
1077 !
1078 ! Program History Log:
1079 ! 91-05-07 Iredell made into inlinable function
1080 ! 94-12-30 Iredell exact computation
1081 ! 1999-03-01 Iredell f90 module
1082 !
1083 ! Usage: tdpl=ftdplxg(tg,pv)
1084 !
1085 ! Input argument list:
1086 ! tg Real(krealfp) guess dewpoint temperature in Kelvin
1087 ! pv Real(krealfp) vapor pressure in Pascals
1088 !
1089 ! Output argument list:
1090 ! ftdplxg Real(krealfp) dewpoint temperature in Kelvin
1091 !
1092 ! Attributes:
1093 ! Language: Fortran 90.
1094 !
1095 !$$$
1096 implicit none
1097 real(krealfp) ftdplxg
1098 real(krealfp),intent(in):: tg,pv
1099 real(krealfp),parameter:: terrm=1.e-6
1100 real(krealfp),parameter:: dldt=con_cvap-con_cliq
1101 real(krealfp),parameter:: heat=con_hvap
1102 real(krealfp),parameter:: xpona=-dldt/con_rv
1103 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
1104 real(krealfp) t,tr,pvt,el,dpvt,terr
1105 integer i
1106 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1107 t=tg
1108 do i=1,100
1109 tr=con_ttp/t
1110 pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
1111 el=heat+dldt*(t-con_ttp)
1112 dpvt=el*pvt/(con_rv*t**2)
1113 terr=(pvt-pv)/dpvt
1114 t=t-terr
1115 if(abs(terr).le.terrm) exit
1116 enddo
1117 ftdplxg=t
1118 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1119 end function
1120 !-------------------------------------------------------------------------------
1121 subroutine gtdpi
1122 !$$$ Subprogram Documentation Block
1123 !
1124 ! Subprogram: gtdpi Compute dewpoint temperature over ice table
1125 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1126 !
1127 ! Abstract: Compute dewpoint temperature table as a function of
1128 ! vapor pressure for inlinable function ftdpi.
1129 ! Exact dewpoint temperatures are calculated in subprogram ftdpixg.
1130 ! The current implementation computes a table with a length
1131 ! of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals
1132 ! giving a dewpoint temperature range of 197 to 279 Kelvin.
1133 !
1134 ! Program History Log:
1135 ! 91-05-07 Iredell
1136 ! 94-12-30 Iredell expand table
1137 ! 1999-03-01 Iredell f90 module
1138 ! 2001-02-26 Iredell ice phase
1139 !
1140 ! Usage: call gtdpi
1141 !
1142 ! Subprograms called:
1143 ! (ftdpixg) inlinable function to compute dewpoint temperature over ice
1144 !
1145 ! Attributes:
1146 ! Language: Fortran 90.
1147 !
1148 !$$$
1149 implicit none
1150 integer jx
1151 real(krealfp) xmin,xmax,xinc,t,x,pv
1152 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1153 xmin=0.1
1154 xmax=1000.1
1155 xinc=(xmax-xmin)/(nxtdpi-1)
1156 c1xtdpi=1.-xmin/xinc
1157 c2xtdpi=1./xinc
1158 t=197.0
1159 do jx=1,nxtdpi
1160 x=xmin+(jx-1)*xinc
1161 pv=x
1162 t=ftdpixg(t,pv)
1163 tbtdpi(jx)=t
1164 enddo
1165 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1166 end subroutine
1167 !-------------------------------------------------------------------------------
1168 ! elemental function ftdpi(pv)
1169 function ftdpi(pv)
1170 !$$$ Subprogram Documentation Block
1171 !
1172 ! Subprogram: ftdpi Compute dewpoint temperature over ice
1173 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1174 !
1175 ! Abstract: Compute dewpoint temperature from vapor pressure.
1176 ! A linear interpolation is done between values in a lookup table
1177 ! computed in gtdpi. See documentation for ftdpixg for details.
1178 ! Input values outside table range are reset to table extrema.
1179 ! The interpolation accuracy is better than 0.0005 Kelvin
1180 ! for dewpoint temperatures greater than 250 Kelvin,
1181 ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
1182 ! On the Cray, ftdpi is about 75 times faster than exact calculation.
1183 ! This function should be expanded inline in the calling routine.
1184 !
1185 ! Program History Log:
1186 ! 91-05-07 Iredell made into inlinable function
1187 ! 94-12-30 Iredell expand table
1188 ! 1999-03-01 Iredell f90 module
1189 ! 2001-02-26 Iredell ice phase
1190 !
1191 ! Usage: tdpi=ftdpi(pv)
1192 !
1193 ! Input argument list:
1194 ! pv Real(krealfp) vapor pressure in Pascals
1195 !
1196 ! Output argument list:
1197 ! ftdpi Real(krealfp) dewpoint temperature in Kelvin
1198 !
1199 ! Attributes:
1200 ! Language: Fortran 90.
1201 !
1202 !$$$
1203 implicit none
1204 real(krealfp) ftdpi
1205 real(krealfp),intent(in):: pv
1206 integer jx
1207 real(krealfp) xj
1208 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1209 xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
1210 jx=min(xj,nxtdpi-1._krealfp)
1211 ftdpi=tbtdpi(jx)+(xj-jx)*(tbtdpi(jx+1)-tbtdpi(jx))
1212 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1213 end function
1214 !-------------------------------------------------------------------------------
1215 ! elemental function ftdpiq(pv)
1216 function ftdpiq(pv)
1217 !$$$ Subprogram Documentation Block
1218 !
1219 ! Subprogram: ftdpiq Compute dewpoint temperature over ice
1220 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1221 !
1222 ! Abstract: Compute dewpoint temperature from vapor pressure.
1223 ! A quadratic interpolation is done between values in a lookup table
1224 ! computed in gtdpi. see documentation for ftdpixg for details.
1225 ! Input values outside table range are reset to table extrema.
1226 ! the interpolation accuracy is better than 0.00001 Kelvin
1227 ! for dewpoint temperatures greater than 250 Kelvin,
1228 ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
1229 ! On the Cray, ftdpiq is about 60 times faster than exact calculation.
1230 ! This function should be expanded inline in the calling routine.
1231 !
1232 ! Program History Log:
1233 ! 91-05-07 Iredell made into inlinable function
1234 ! 94-12-30 Iredell quadratic interpolation
1235 ! 1999-03-01 Iredell f90 module
1236 ! 2001-02-26 Iredell ice phase
1237 !
1238 ! Usage: tdpi=ftdpiq(pv)
1239 !
1240 ! Input argument list:
1241 ! pv Real(krealfp) vapor pressure in Pascals
1242 !
1243 ! Output argument list:
1244 ! ftdpiq Real(krealfp) dewpoint temperature in Kelvin
1245 !
1246 ! Attributes:
1247 ! Language: Fortran 90.
1248 !
1249 !$$$
1250 implicit none
1251 real(krealfp) ftdpiq
1252 real(krealfp),intent(in):: pv
1253 integer jx
1254 real(krealfp) xj,dxj,fj1,fj2,fj3
1255 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1256 xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
1257 jx=min(max(nint(xj),2),nxtdpi-1)
1258 dxj=xj-jx
1259 fj1=tbtdpi(jx-1)
1260 fj2=tbtdpi(jx)
1261 fj3=tbtdpi(jx+1)
1262 ftdpiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1263 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1264 end function
1265 !-------------------------------------------------------------------------------
1266 ! elemental function ftdpix(pv)
1267 function ftdpix(pv)
1268 !$$$ Subprogram Documentation Block
1269 !
1270 ! Subprogram: ftdpix Compute dewpoint temperature over ice
1271 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1272 !
1273 ! Abstract: exactly compute dewpoint temperature from vapor pressure.
1274 ! An approximate dewpoint temperature for function ftdpixg
1275 ! is obtained using ftdpi so gtdpi must be already called.
1276 ! See documentation for ftdpixg for details.
1277 !
1278 ! Program History Log:
1279 ! 91-05-07 Iredell made into inlinable function
1280 ! 94-12-30 Iredell exact computation
1281 ! 1999-03-01 Iredell f90 module
1282 ! 2001-02-26 Iredell ice phase
1283 !
1284 ! Usage: tdpi=ftdpix(pv)
1285 !
1286 ! Input argument list:
1287 ! pv Real(krealfp) vapor pressure in Pascals
1288 !
1289 ! Output argument list:
1290 ! ftdpix Real(krealfp) dewpoint temperature in Kelvin
1291 !
1292 ! Subprograms called:
1293 ! (ftdpi) inlinable function to compute dewpoint temperature over ice
1294 ! (ftdpixg) inlinable function to compute dewpoint temperature over ice
1295 !
1296 ! Attributes:
1297 ! Language: Fortran 90.
1298 !
1299 !$$$
1300 implicit none
1301 real(krealfp) ftdpix
1302 real(krealfp),intent(in):: pv
1303 real(krealfp) tg
1304 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1305 tg=ftdpi(pv)
1306 ftdpix=ftdpixg(tg,pv)
1307 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1308 end function
1309 !-------------------------------------------------------------------------------
1310 ! elemental function ftdpixg(tg,pv)
1311 function ftdpixg(tg,pv)
1312 !$$$ Subprogram Documentation Block
1313 !
1314 ! Subprogram: ftdpixg Compute dewpoint temperature over ice
1315 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1316 !
1317 ! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1318 ! A guess dewpoint temperature must be provided.
1319 ! The water model assumes a perfect gas, constant specific heats
1320 ! for gas and ice, and neglects the volume of the ice.
1321 ! The model does account for the variation of the latent heat
1322 ! of sublimation with temperature. The liquid option is not included.
1323 ! The Clausius-Clapeyron equation is integrated from the triple point
1324 ! to get the formula
1325 ! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
1326 ! where tr is ttp/t and other values are physical constants.
1327 ! The formula is inverted by iterating Newtonian approximations
1328 ! for each pvs until t is found to within 1.e-6 Kelvin.
1329 ! This function can be expanded inline in the calling routine.
1330 !
1331 ! Program History Log:
1332 ! 91-05-07 Iredell made into inlinable function
1333 ! 94-12-30 Iredell exact computation
1334 ! 1999-03-01 Iredell f90 module
1335 ! 2001-02-26 Iredell ice phase
1336 !
1337 ! Usage: tdpi=ftdpixg(tg,pv)
1338 !
1339 ! Input argument list:
1340 ! tg Real(krealfp) guess dewpoint temperature in Kelvin
1341 ! pv Real(krealfp) vapor pressure in Pascals
1342 !
1343 ! Output argument list:
1344 ! ftdpixg Real(krealfp) dewpoint temperature in Kelvin
1345 !
1346 ! Attributes:
1347 ! Language: Fortran 90.
1348 !
1349 !$$$
1350 implicit none
1351 real(krealfp) ftdpixg
1352 real(krealfp),intent(in):: tg,pv
1353 real(krealfp),parameter:: terrm=1.e-6
1354 real(krealfp),parameter:: dldt=con_cvap-con_csol
1355 real(krealfp),parameter:: heat=con_hvap+con_hfus
1356 real(krealfp),parameter:: xpona=-dldt/con_rv
1357 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
1358 real(krealfp) t,tr,pvt,el,dpvt,terr
1359 integer i
1360 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1361 t=tg
1362 do i=1,100
1363 tr=con_ttp/t
1364 pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
1365 el=heat+dldt*(t-con_ttp)
1366 dpvt=el*pvt/(con_rv*t**2)
1367 terr=(pvt-pv)/dpvt
1368 t=t-terr
1369 if(abs(terr).le.terrm) exit
1370 enddo
1371 ftdpixg=t
1372 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1373 end function
1374 !-------------------------------------------------------------------------------
1375 subroutine gtdp
1376 !$$$ Subprogram Documentation Block
1377 !
1378 ! Subprogram: gtdp Compute dewpoint temperature table
1379 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1380 !
1381 ! Abstract: Compute dewpoint temperature table as a function of
1382 ! vapor pressure for inlinable function ftdp.
1383 ! Exact dewpoint temperatures are calculated in subprogram ftdpxg.
1384 ! The current implementation computes a table with a length
1385 ! of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals
1386 ! giving a dewpoint temperature range of 208 to 319 Kelvin.
1387 !
1388 ! Program History Log:
1389 ! 91-05-07 Iredell
1390 ! 94-12-30 Iredell expand table
1391 ! 1999-03-01 Iredell f90 module
1392 ! 2001-02-26 Iredell ice phase
1393 !
1394 ! Usage: call gtdp
1395 !
1396 ! Subprograms called:
1397 ! (ftdpxg) inlinable function to compute dewpoint temperature
1398 !
1399 ! Attributes:
1400 ! Language: Fortran 90.
1401 !
1402 !$$$
1403 implicit none
1404 integer jx
1405 real(krealfp) xmin,xmax,xinc,t,x,pv
1406 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1407 xmin=0.5
1408 xmax=10000.5
1409 xinc=(xmax-xmin)/(nxtdp-1)
1410 c1xtdp=1.-xmin/xinc
1411 c2xtdp=1./xinc
1412 t=208.0
1413 do jx=1,nxtdp
1414 x=xmin+(jx-1)*xinc
1415 pv=x
1416 t=ftdpxg(t,pv)
1417 tbtdp(jx)=t
1418 enddo
1419 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1420 end subroutine
1421 !-------------------------------------------------------------------------------
1422 ! elemental function ftdp(pv)
1423 function ftdp(pv)
1424 !$$$ Subprogram Documentation Block
1425 !
1426 ! Subprogram: ftdp Compute dewpoint temperature
1427 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1428 !
1429 ! Abstract: Compute dewpoint temperature from vapor pressure.
1430 ! A linear interpolation is done between values in a lookup table
1431 ! computed in gtdp. See documentation for ftdpxg for details.
1432 ! Input values outside table range are reset to table extrema.
1433 ! The interpolation accuracy is better than 0.0005 Kelvin
1434 ! for dewpoint temperatures greater than 250 Kelvin,
1435 ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
1436 ! On the Cray, ftdp is about 75 times faster than exact calculation.
1437 ! This function should be expanded inline in the calling routine.
1438 !
1439 ! Program History Log:
1440 ! 91-05-07 Iredell made into inlinable function
1441 ! 94-12-30 Iredell expand table
1442 ! 1999-03-01 Iredell f90 module
1443 ! 2001-02-26 Iredell ice phase
1444 !
1445 ! Usage: tdp=ftdp(pv)
1446 !
1447 ! Input argument list:
1448 ! pv Real(krealfp) vapor pressure in Pascals
1449 !
1450 ! Output argument list:
1451 ! ftdp Real(krealfp) dewpoint temperature in Kelvin
1452 !
1453 ! Attributes:
1454 ! Language: Fortran 90.
1455 !
1456 !$$$
1457 implicit none
1458 real(krealfp) ftdp
1459 real(krealfp),intent(in):: pv
1460 integer jx
1461 real(krealfp) xj
1462 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1463 xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
1464 jx=min(xj,nxtdp-1._krealfp)
1465 ftdp=tbtdp(jx)+(xj-jx)*(tbtdp(jx+1)-tbtdp(jx))
1466 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1467 end function
1468 !-------------------------------------------------------------------------------
1469 ! elemental function ftdpq(pv)
1470 function ftdpq(pv)
1471 !$$$ Subprogram Documentation Block
1472 !
1473 ! Subprogram: ftdpq Compute dewpoint temperature
1474 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1475 !
1476 ! Abstract: Compute dewpoint temperature from vapor pressure.
1477 ! A quadratic interpolation is done between values in a lookup table
1478 ! computed in gtdp. see documentation for ftdpxg for details.
1479 ! Input values outside table range are reset to table extrema.
1480 ! the interpolation accuracy is better than 0.00001 Kelvin
1481 ! for dewpoint temperatures greater than 250 Kelvin,
1482 ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
1483 ! On the Cray, ftdpq is about 60 times faster than exact calculation.
1484 ! This function should be expanded inline in the calling routine.
1485 !
1486 ! Program History Log:
1487 ! 91-05-07 Iredell made into inlinable function
1488 ! 94-12-30 Iredell quadratic interpolation
1489 ! 1999-03-01 Iredell f90 module
1490 ! 2001-02-26 Iredell ice phase
1491 !
1492 ! Usage: tdp=ftdpq(pv)
1493 !
1494 ! Input argument list:
1495 ! pv Real(krealfp) vapor pressure in Pascals
1496 !
1497 ! Output argument list:
1498 ! ftdpq Real(krealfp) dewpoint temperature in Kelvin
1499 !
1500 ! Attributes:
1501 ! Language: Fortran 90.
1502 !
1503 !$$$
1504 implicit none
1505 real(krealfp) ftdpq
1506 real(krealfp),intent(in):: pv
1507 integer jx
1508 real(krealfp) xj,dxj,fj1,fj2,fj3
1509 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1510 xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
1511 jx=min(max(nint(xj),2),nxtdp-1)
1512 dxj=xj-jx
1513 fj1=tbtdp(jx-1)
1514 fj2=tbtdp(jx)
1515 fj3=tbtdp(jx+1)
1516 ftdpq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1517 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1518 end function
1519 !-------------------------------------------------------------------------------
1520 ! elemental function ftdpx(pv)
1521 function ftdpx(pv)
1522 !$$$ Subprogram Documentation Block
1523 !
1524 ! Subprogram: ftdpx Compute dewpoint temperature
1525 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1526 !
1527 ! Abstract: exactly compute dewpoint temperature from vapor pressure.
1528 ! An approximate dewpoint temperature for function ftdpxg
1529 ! is obtained using ftdp so gtdp must be already called.
1530 ! See documentation for ftdpxg for details.
1531 !
1532 ! Program History Log:
1533 ! 91-05-07 Iredell made into inlinable function
1534 ! 94-12-30 Iredell exact computation
1535 ! 1999-03-01 Iredell f90 module
1536 ! 2001-02-26 Iredell ice phase
1537 !
1538 ! Usage: tdp=ftdpx(pv)
1539 !
1540 ! Input argument list:
1541 ! pv Real(krealfp) vapor pressure in Pascals
1542 !
1543 ! Output argument list:
1544 ! ftdpx Real(krealfp) dewpoint temperature in Kelvin
1545 !
1546 ! Subprograms called:
1547 ! (ftdp) inlinable function to compute dewpoint temperature
1548 ! (ftdpxg) inlinable function to compute dewpoint temperature
1549 !
1550 ! Attributes:
1551 ! Language: Fortran 90.
1552 !
1553 !$$$
1554 implicit none
1555 real(krealfp) ftdpx
1556 real(krealfp),intent(in):: pv
1557 real(krealfp) tg
1558 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1559 tg=ftdp(pv)
1560 ftdpx=ftdpxg(tg,pv)
1561 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1562 end function
1563 !-------------------------------------------------------------------------------
1564 ! elemental function ftdpxg(tg,pv)
1565 function ftdpxg(tg,pv)
1566 !$$$ Subprogram Documentation Block
1567 !
1568 ! Subprogram: ftdpxg Compute dewpoint temperature
1569 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1570 !
1571 ! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1572 ! A guess dewpoint temperature must be provided.
1573 ! The saturation vapor pressure over either liquid and ice is computed
1574 ! over liquid for temperatures above the triple point,
1575 ! over ice for temperatures 20 degress below the triple point,
1576 ! and a linear combination of the two for temperatures in between.
1577 ! The water model assumes a perfect gas, constant specific heats
1578 ! for gas, liquid and ice, and neglects the volume of the condensate.
1579 ! The model does account for the variation of the latent heat
1580 ! of condensation and sublimation with temperature.
1581 ! The Clausius-Clapeyron equation is integrated from the triple point
1582 ! to get the formula
1583 ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
1584 ! where tr is ttp/t and other values are physical constants.
1585 ! The reference for this decision is Emanuel(1994), pages 116-117.
1586 ! The formula is inverted by iterating Newtonian approximations
1587 ! for each pvs until t is found to within 1.e-6 Kelvin.
1588 ! This function can be expanded inline in the calling routine.
1589 !
1590 ! Program History Log:
1591 ! 91-05-07 Iredell made into inlinable function
1592 ! 94-12-30 Iredell exact computation
1593 ! 1999-03-01 Iredell f90 module
1594 ! 2001-02-26 Iredell ice phase
1595 !
1596 ! Usage: tdp=ftdpxg(tg,pv)
1597 !
1598 ! Input argument list:
1599 ! tg Real(krealfp) guess dewpoint temperature in Kelvin
1600 ! pv Real(krealfp) vapor pressure in Pascals
1601 !
1602 ! Output argument list:
1603 ! ftdpxg Real(krealfp) dewpoint temperature in Kelvin
1604 !
1605 ! Attributes:
1606 ! Language: Fortran 90.
1607 !
1608 !$$$
1609 implicit none
1610 real(krealfp) ftdpxg
1611 real(krealfp),intent(in):: tg,pv
1612 real(krealfp),parameter:: terrm=1.e-6
1613 real(krealfp),parameter:: tliq=con_ttp
1614 real(krealfp),parameter:: tice=con_ttp-20.0
1615 real(krealfp),parameter:: dldtl=con_cvap-con_cliq
1616 real(krealfp),parameter:: heatl=con_hvap
1617 real(krealfp),parameter:: xponal=-dldtl/con_rv
1618 real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
1619 real(krealfp),parameter:: dldti=con_cvap-con_csol
1620 real(krealfp),parameter:: heati=con_hvap+con_hfus
1621 real(krealfp),parameter:: xponai=-dldti/con_rv
1622 real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
1623 real(krealfp) t,tr,w,pvtl,pvti,pvt,ell,eli,el,dpvt,terr
1624 integer i
1625 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1626 t=tg
1627 do i=1,100
1628 tr=con_ttp/t
1629 if(t.ge.tliq) then
1630 pvt=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1631 el=heatl+dldtl*(t-con_ttp)
1632 dpvt=el*pvt/(con_rv*t**2)
1633 elseif(t.lt.tice) then
1634 pvt=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1635 el=heati+dldti*(t-con_ttp)
1636 dpvt=el*pvt/(con_rv*t**2)
1637 else
1638 w=(t-tice)/(tliq-tice)
1639 pvtl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1640 pvti=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1641 pvt=w*pvtl+(1.-w)*pvti
1642 ell=heatl+dldtl*(t-con_ttp)
1643 eli=heati+dldti*(t-con_ttp)
1644 dpvt=(w*ell*pvtl+(1.-w)*eli*pvti)/(con_rv*t**2)
1645 endif
1646 terr=(pvt-pv)/dpvt
1647 t=t-terr
1648 if(abs(terr).le.terrm) exit
1649 enddo
1650 ftdpxg=t
1651 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1652 end function
1653 !-------------------------------------------------------------------------------
1654 subroutine gthe
1655 !$$$ Subprogram Documentation Block
1656 !
1657 ! Subprogram: gthe Compute equivalent potential temperature table
1658 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1659 !
1660 ! Abstract: Compute equivalent potential temperature table
1661 ! as a function of LCL temperature and pressure over 1e5 Pa
1662 ! to the kappa power for function fthe.
1663 ! Equivalent potential temperatures are calculated in subprogram fthex
1664 ! the current implementation computes a table with a first dimension
1665 ! of 241 for temperatures ranging from 183.16 to 303.16 Kelvin
1666 ! and a second dimension of 151 for pressure over 1e5 Pa
1667 ! to the kappa power ranging from 0.04**rocp to 1.10**rocp.
1668 !
1669 ! Program History Log:
1670 ! 91-05-07 Iredell
1671 ! 94-12-30 Iredell expand table
1672 ! 1999-03-01 Iredell f90 module
1673 !
1674 ! Usage: call gthe
1675 !
1676 ! Subprograms called:
1677 ! (fthex) inlinable function to compute equiv. pot. temperature
1678 !
1679 ! Attributes:
1680 ! Language: Fortran 90.
1681 !
1682 !$$$
1683 implicit none
1684 integer jx,jy
1685 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,t
1686 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1687 xmin=con_ttp-90._krealfp
1688 xmax=con_ttp+30._krealfp
1689 ymin=0.04_krealfp**con_rocp
1690 ymax=1.10_krealfp**con_rocp
1691 xinc=(xmax-xmin)/(nxthe-1)
1692 c1xthe=1.-xmin/xinc
1693 c2xthe=1./xinc
1694 yinc=(ymax-ymin)/(nythe-1)
1695 c1ythe=1.-ymin/yinc
1696 c2ythe=1./yinc
1697 do jy=1,nythe
1698 y=ymin+(jy-1)*yinc
1699 pk=y
1700 do jx=1,nxthe
1701 x=xmin+(jx-1)*xinc
1702 t=x
1703 tbthe(jx,jy)=fthex(t,pk)
1704 enddo
1705 enddo
1706 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1707 end subroutine
1708 !-------------------------------------------------------------------------------
1709 ! elemental function fthe(t,pk)
1710 function fthe(t,pk)
1711 !$$$ Subprogram Documentation Block
1712 !
1713 ! Subprogram: fthe Compute equivalent potential temperature
1714 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1715 !
1716 ! Abstract: Compute equivalent potential temperature at the LCL
1717 ! from temperature and pressure over 1e5 Pa to the kappa power.
1718 ! A bilinear interpolation is done between values in a lookup table
1719 ! computed in gthe. see documentation for fthex for details.
1720 ! Input values outside table range are reset to table extrema,
1721 ! except zero is returned for too cold or high LCLs.
1722 ! The interpolation accuracy is better than 0.01 Kelvin.
1723 ! On the Cray, fthe is almost 6 times faster than exact calculation.
1724 ! This function should be expanded inline in the calling routine.
1725 !
1726 ! Program History Log:
1727 ! 91-05-07 Iredell made into inlinable function
1728 ! 94-12-30 Iredell expand table
1729 ! 1999-03-01 Iredell f90 module
1730 !
1731 ! Usage: the=fthe(t,pk)
1732 !
1733 ! Input argument list:
1734 ! t Real(krealfp) LCL temperature in Kelvin
1735 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
1736 !
1737 ! Output argument list:
1738 ! fthe Real(krealfp) equivalent potential temperature in Kelvin
1739 !
1740 ! Attributes:
1741 ! Language: Fortran 90.
1742 !
1743 !$$$
1744 implicit none
1745 real(krealfp) fthe
1746 real(krealfp),intent(in):: t,pk
1747 integer jx,jy
1748 real(krealfp) xj,yj,ftx1,ftx2
1749 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1750 xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
1751 yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
1752 if(xj.ge.1..and.yj.ge.1.) then
1753 jx=min(xj,nxthe-1._krealfp)
1754 jy=min(yj,nythe-1._krealfp)
1755 ftx1=tbthe(jx,jy)+(xj-jx)*(tbthe(jx+1,jy)-tbthe(jx,jy))
1756 ftx2=tbthe(jx,jy+1)+(xj-jx)*(tbthe(jx+1,jy+1)-tbthe(jx,jy+1))
1757 fthe=ftx1+(yj-jy)*(ftx2-ftx1)
1758 else
1759 fthe=0.
1760 endif
1761 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1762 end function
1763 !-------------------------------------------------------------------------------
1764 ! elemental function ftheq(t,pk)
1765 function ftheq(t,pk)
1766 !$$$ Subprogram Documentation Block
1767 !
1768 ! Subprogram: ftheq Compute equivalent potential temperature
1769 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1770 !
1771 ! Abstract: Compute equivalent potential temperature at the LCL
1772 ! from temperature and pressure over 1e5 Pa to the kappa power.
1773 ! A biquadratic interpolation is done between values in a lookup table
1774 ! computed in gthe. see documentation for fthex for details.
1775 ! Input values outside table range are reset to table extrema,
1776 ! except zero is returned for too cold or high LCLs.
1777 ! The interpolation accuracy is better than 0.0002 Kelvin.
1778 ! On the Cray, ftheq is almost 3 times faster than exact calculation.
1779 ! This function should be expanded inline in the calling routine.
1780 !
1781 ! Program History Log:
1782 ! 91-05-07 Iredell made into inlinable function
1783 ! 94-12-30 Iredell quadratic interpolation
1784 ! 1999-03-01 Iredell f90 module
1785 !
1786 ! Usage: the=ftheq(t,pk)
1787 !
1788 ! Input argument list:
1789 ! t Real(krealfp) LCL temperature in Kelvin
1790 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
1791 !
1792 ! Output argument list:
1793 ! ftheq Real(krealfp) equivalent potential temperature in Kelvin
1794 !
1795 ! Attributes:
1796 ! Language: Fortran 90.
1797 !
1798 !$$$
1799 implicit none
1800 real(krealfp) ftheq
1801 real(krealfp),intent(in):: t,pk
1802 integer jx,jy
1803 real(krealfp) xj,yj,dxj,dyj
1804 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
1805 real(krealfp) ftx1,ftx2,ftx3
1806 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1807 xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
1808 yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
1809 if(xj.ge.1..and.yj.ge.1.) then
1810 jx=min(max(nint(xj),2),nxthe-1)
1811 jy=min(max(nint(yj),2),nythe-1)
1812 dxj=xj-jx
1813 dyj=yj-jy
1814 ft11=tbthe(jx-1,jy-1)
1815 ft12=tbthe(jx-1,jy)
1816 ft13=tbthe(jx-1,jy+1)
1817 ft21=tbthe(jx,jy-1)
1818 ft22=tbthe(jx,jy)
1819 ft23=tbthe(jx,jy+1)
1820 ft31=tbthe(jx+1,jy-1)
1821 ft32=tbthe(jx+1,jy)
1822 ft33=tbthe(jx+1,jy+1)
1823 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
1824 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
1825 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
1826 ftheq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
1827 else
1828 ftheq=0.
1829 endif
1830 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1831 end function
1832 !-------------------------------------------------------------------------------
1833 ! elemental function fthex(t,pk)
1834 function fthex(t,pk)
1835 !$$$ Subprogram Documentation Block
1836 !
1837 ! Subprogram: fthex Compute equivalent potential temperature
1838 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1839 !
1840 ! Abstract: Exactly compute equivalent potential temperature at the LCL
1841 ! from temperature and pressure over 1e5 Pa to the kappa power.
1842 ! Equivalent potential temperature is constant for a saturated parcel
1843 ! rising adiabatically up a moist adiabat when the heat and mass
1844 ! of the condensed water are neglected. Ice is also neglected.
1845 ! The formula for equivalent potential temperature (Holton) is
1846 ! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
1847 ! where t is the temperature, pv is the saturated vapor pressure,
1848 ! pd is the dry pressure p-pv, el is the temperature dependent
1849 ! latent heat of condensation hvap+dldt*(t-ttp), and other values
1850 ! are physical constants defined in parameter statements in the code.
1851 ! Zero is returned if the input values make saturation impossible.
1852 ! This function should be expanded inline in the calling routine.
1853 !
1854 ! Program History Log:
1855 ! 91-05-07 Iredell made into inlinable function
1856 ! 94-12-30 Iredell exact computation
1857 ! 1999-03-01 Iredell f90 module
1858 !
1859 ! Usage: the=fthex(t,pk)
1860 !
1861 ! Input argument list:
1862 ! t Real(krealfp) LCL temperature in Kelvin
1863 ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
1864 !
1865 ! Output argument list:
1866 ! fthex Real(krealfp) equivalent potential temperature in Kelvin
1867 !
1868 ! Attributes:
1869 ! Language: Fortran 90.
1870 !
1871 !$$$
1872 implicit none
1873 real(krealfp) fthex
1874 real(krealfp),intent(in):: t,pk
1875 real(krealfp) p,tr,pv,pd,el,expo,expmax
1876 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1877 p=pk**con_cpor
1878 tr=con_ttp/t
1879 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
1880 pd=p-pv
1881 if(pd.gt.pv) then
1882 el=con_hvap+con_dldt*(t-con_ttp)
1883 expo=el*con_eps*pv/(con_cp*t*pd)
1884 fthex=t*pd**(-con_rocp)*exp(expo)
1885 else
1886 fthex=0.
1887 endif
1888 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1889 end function
1890 !-------------------------------------------------------------------------------
1891 subroutine gtma
1892 !$$$ Subprogram Documentation Block
1893 !
1894 ! Subprogram: gtma Compute moist adiabat tables
1895 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1896 !
1897 ! Abstract: Compute temperature and specific humidity tables
1898 ! as a function of equivalent potential temperature and
1899 ! pressure over 1e5 Pa to the kappa power for subprogram stma.
1900 ! Exact parcel temperatures are calculated in subprogram stmaxg.
1901 ! The current implementation computes a table with a first dimension
1902 ! of 151 for equivalent potential temperatures ranging from 200 to 500
1903 ! Kelvin and a second dimension of 121 for pressure over 1e5 Pa
1904 ! to the kappa power ranging from 0.01**rocp to 1.10**rocp.
1905 !
1906 ! Program History Log:
1907 ! 91-05-07 Iredell
1908 ! 94-12-30 Iredell expand table
1909 ! 1999-03-01 Iredell f90 module
1910 !
1911 ! Usage: call gtma
1912 !
1913 ! Subprograms called:
1914 ! (stmaxg) inlinable subprogram to compute parcel temperature
1915 !
1916 ! Attributes:
1917 ! Language: Fortran 90.
1918 !
1919 !$$$
1920 implicit none
1921 integer jx,jy
1922 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,the,t,q,tg
1923 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1924 xmin=200._krealfp
1925 xmax=500._krealfp
1926 ymin=0.01_krealfp**con_rocp
1927 ymax=1.10_krealfp**con_rocp
1928 xinc=(xmax-xmin)/(nxma-1)
1929 c1xma=1.-xmin/xinc
1930 c2xma=1./xinc
1931 yinc=(ymax-ymin)/(nyma-1)
1932 c1yma=1.-ymin/yinc
1933 c2yma=1./yinc
1934 do jy=1,nyma
1935 y=ymin+(jy-1)*yinc
1936 pk=y
1937 tg=xmin*y
1938 do jx=1,nxma
1939 x=xmin+(jx-1)*xinc
1940 the=x
1941 call stmaxg(tg,the,pk,t,q)
1942 tbtma(jx,jy)=t
1943 tbqma(jx,jy)=q
1944 tg=t
1945 enddo
1946 enddo
1947 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1948 end subroutine
1949 !-------------------------------------------------------------------------------
1950 ! elemental subroutine stma(the,pk,tma,qma)
1951 subroutine stma(the,pk,tma,qma)
1952 !$$$ Subprogram Documentation Block
1953 !
1954 ! Subprogram: stma Compute moist adiabat temperature
1955 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1956 !
1957 ! Abstract: Compute temperature and specific humidity of a parcel
1958 ! lifted up a moist adiabat from equivalent potential temperature
1959 ! at the LCL and pressure over 1e5 Pa to the kappa power.
1960 ! Bilinear interpolations are done between values in a lookup table
1961 ! computed in gtma. See documentation for stmaxg for details.
1962 ! Input values outside table range are reset to table extrema.
1963 ! The interpolation accuracy is better than 0.01 Kelvin
1964 ! and 5.e-6 kg/kg for temperature and humidity, respectively.
1965 ! On the Cray, stma is about 35 times faster than exact calculation.
1966 ! This subprogram should be expanded inline in the calling routine.
1967 !
1968 ! Program History Log:
1969 ! 91-05-07 Iredell made into inlinable function
1970 ! 94-12-30 Iredell expand table
1971 ! 1999-03-01 Iredell f90 module
1972 !
1973 ! Usage: call stma(the,pk,tma,qma)
1974 !
1975 ! Input argument list:
1976 ! the Real(krealfp) equivalent potential temperature in Kelvin
1977 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
1978 !
1979 ! Output argument list:
1980 ! tma Real(krealfp) parcel temperature in Kelvin
1981 ! qma Real(krealfp) parcel specific humidity in kg/kg
1982 !
1983 ! Attributes:
1984 ! Language: Fortran 90.
1985 !
1986 !$$$
1987 implicit none
1988 real(krealfp),intent(in):: the,pk
1989 real(krealfp),intent(out):: tma,qma
1990 integer jx,jy
1991 real(krealfp) xj,yj,ftx1,ftx2,qx1,qx2
1992 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1993 xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
1994 yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
1995 jx=min(xj,nxma-1._krealfp)
1996 jy=min(yj,nyma-1._krealfp)
1997 ftx1=tbtma(jx,jy)+(xj-jx)*(tbtma(jx+1,jy)-tbtma(jx,jy))
1998 ftx2=tbtma(jx,jy+1)+(xj-jx)*(tbtma(jx+1,jy+1)-tbtma(jx,jy+1))
1999 tma=ftx1+(yj-jy)*(ftx2-ftx1)
2000 qx1=tbqma(jx,jy)+(xj-jx)*(tbqma(jx+1,jy)-tbqma(jx,jy))
2001 qx2=tbqma(jx,jy+1)+(xj-jx)*(tbqma(jx+1,jy+1)-tbqma(jx,jy+1))
2002 qma=qx1+(yj-jy)*(qx2-qx1)
2003 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2004 end subroutine
2005 !-------------------------------------------------------------------------------
2006 ! elemental subroutine stmaq(the,pk,tma,qma)
2007 subroutine stmaq(the,pk,tma,qma)
2008 !$$$ Subprogram Documentation Block
2009 !
2010 ! Subprogram: stmaq Compute moist adiabat temperature
2011 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2012 !
2013 ! Abstract: Compute temperature and specific humidity of a parcel
2014 ! lifted up a moist adiabat from equivalent potential temperature
2015 ! at the LCL and pressure over 1e5 Pa to the kappa power.
2016 ! Biquadratic interpolations are done between values in a lookup table
2017 ! computed in gtma. See documentation for stmaxg for details.
2018 ! Input values outside table range are reset to table extrema.
2019 ! the interpolation accuracy is better than 0.0005 Kelvin
2020 ! and 1.e-7 kg/kg for temperature and humidity, respectively.
2021 ! On the Cray, stmaq is about 25 times faster than exact calculation.
2022 ! This subprogram should be expanded inline in the calling routine.
2023 !
2024 ! Program History Log:
2025 ! 91-05-07 Iredell made into inlinable function
2026 ! 94-12-30 Iredell quadratic interpolation
2027 ! 1999-03-01 Iredell f90 module
2028 !
2029 ! Usage: call stmaq(the,pk,tma,qma)
2030 !
2031 ! Input argument list:
2032 ! the Real(krealfp) equivalent potential temperature in Kelvin
2033 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2034 !
2035 ! Output argument list:
2036 ! tmaq Real(krealfp) parcel temperature in Kelvin
2037 ! qma Real(krealfp) parcel specific humidity in kg/kg
2038 !
2039 ! Attributes:
2040 ! Language: Fortran 90.
2041 !
2042 !$$$
2043 implicit none
2044 real(krealfp),intent(in):: the,pk
2045 real(krealfp),intent(out):: tma,qma
2046 integer jx,jy
2047 real(krealfp) xj,yj,dxj,dyj
2048 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
2049 real(krealfp) ftx1,ftx2,ftx3
2050 real(krealfp) q11,q12,q13,q21,q22,q23,q31,q32,q33,qx1,qx2,qx3
2051 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2052 xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
2053 yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
2054 jx=min(max(nint(xj),2),nxma-1)
2055 jy=min(max(nint(yj),2),nyma-1)
2056 dxj=xj-jx
2057 dyj=yj-jy
2058 ft11=tbtma(jx-1,jy-1)
2059 ft12=tbtma(jx-1,jy)
2060 ft13=tbtma(jx-1,jy+1)
2061 ft21=tbtma(jx,jy-1)
2062 ft22=tbtma(jx,jy)
2063 ft23=tbtma(jx,jy+1)
2064 ft31=tbtma(jx+1,jy-1)
2065 ft32=tbtma(jx+1,jy)
2066 ft33=tbtma(jx+1,jy+1)
2067 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
2068 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
2069 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
2070 tma=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
2071 q11=tbqma(jx-1,jy-1)
2072 q12=tbqma(jx-1,jy)
2073 q13=tbqma(jx-1,jy+1)
2074 q21=tbqma(jx,jy-1)
2075 q22=tbqma(jx,jy)
2076 q23=tbqma(jx,jy+1)
2077 q31=tbqma(jx+1,jy-1)
2078 q32=tbqma(jx+1,jy)
2079 q33=tbqma(jx+1,jy+1)
2080 qx1=(((q31+q11)/2-q21)*dxj+(q31-q11)/2)*dxj+q21
2081 qx2=(((q32+q12)/2-q22)*dxj+(q32-q12)/2)*dxj+q22
2082 qx3=(((q33+q13)/2-q23)*dxj+(q33-q13)/2)*dxj+q23
2083 qma=(((qx3+qx1)/2-qx2)*dyj+(qx3-qx1)/2)*dyj+qx2
2084 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2085 end subroutine
2086 !-------------------------------------------------------------------------------
2087 ! elemental subroutine stmax(the,pk,tma,qma)
2088 subroutine stmax(the,pk,tma,qma)
2089 !$$$ Subprogram Documentation Block
2090 !
2091 ! Subprogram: stmax Compute moist adiabat temperature
2092 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2093 !
2094 ! Abstract: Exactly compute temperature and humidity of a parcel
2095 ! lifted up a moist adiabat from equivalent potential temperature
2096 ! at the LCL and pressure over 1e5 Pa to the kappa power.
2097 ! An approximate parcel temperature for subprogram stmaxg
2098 ! is obtained using stma so gtma must be already called.
2099 ! See documentation for stmaxg for details.
2100 !
2101 ! Program History Log:
2102 ! 91-05-07 Iredell made into inlinable function
2103 ! 94-12-30 Iredell exact computation
2104 ! 1999-03-01 Iredell f90 module
2105 !
2106 ! Usage: call stmax(the,pk,tma,qma)
2107 !
2108 ! Input argument list:
2109 ! the Real(krealfp) equivalent potential temperature in Kelvin
2110 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2111 !
2112 ! Output argument list:
2113 ! tma Real(krealfp) parcel temperature in Kelvin
2114 ! qma Real(krealfp) parcel specific humidity in kg/kg
2115 !
2116 ! Subprograms called:
2117 ! (stma) inlinable subprogram to compute parcel temperature
2118 ! (stmaxg) inlinable subprogram to compute parcel temperature
2119 !
2120 ! Attributes:
2121 ! Language: Fortran 90.
2122 !
2123 !$$$
2124 implicit none
2125 real(krealfp),intent(in):: the,pk
2126 real(krealfp),intent(out):: tma,qma
2127 real(krealfp) tg,qg
2128 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2129 call stma(the,pk,tg,qg)
2130 call stmaxg(tg,the,pk,tma,qma)
2131 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2132 end subroutine
2133 !-------------------------------------------------------------------------------
2134 ! elemental subroutine stmaxg(tg,the,pk,tma,qma)
2135 subroutine stmaxg(tg,the,pk,tma,qma)
2136 !$$$ Subprogram Documentation Block
2137 !
2138 ! Subprogram: stmaxg Compute moist adiabat temperature
2139 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2140 !
2141 ! Abstract: exactly compute temperature and humidity of a parcel
2142 ! lifted up a moist adiabat from equivalent potential temperature
2143 ! at the LCL and pressure over 1e5 Pa to the kappa power.
2144 ! A guess parcel temperature must be provided.
2145 ! Equivalent potential temperature is constant for a saturated parcel
2146 ! rising adiabatically up a moist adiabat when the heat and mass
2147 ! of the condensed water are neglected. Ice is also neglected.
2148 ! The formula for equivalent potential temperature (Holton) is
2149 ! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
2150 ! where t is the temperature, pv is the saturated vapor pressure,
2151 ! pd is the dry pressure p-pv, el is the temperature dependent
2152 ! latent heat of condensation hvap+dldt*(t-ttp), and other values
2153 ! are physical constants defined in parameter statements in the code.
2154 ! The formula is inverted by iterating Newtonian approximations
2155 ! for each the and p until t is found to within 1.e-4 Kelvin.
2156 ! The specific humidity is then computed from pv and pd.
2157 ! This subprogram can be expanded inline in the calling routine.
2158 !
2159 ! Program History Log:
2160 ! 91-05-07 Iredell made into inlinable function
2161 ! 94-12-30 Iredell exact computation
2162 ! 1999-03-01 Iredell f90 module
2163 !
2164 ! Usage: call stmaxg(tg,the,pk,tma,qma)
2165 !
2166 ! Input argument list:
2167 ! tg Real(krealfp) guess parcel temperature in Kelvin
2168 ! the Real(krealfp) equivalent potential temperature in Kelvin
2169 ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2170 !
2171 ! Output argument list:
2172 ! tma Real(krealfp) parcel temperature in Kelvin
2173 ! qma Real(krealfp) parcel specific humidity in kg/kg
2174 !
2175 ! Attributes:
2176 ! Language: Fortran 90.
2177 !
2178 !$$$
2179 implicit none
2180 real(krealfp),intent(in):: tg,the,pk
2181 real(krealfp),intent(out):: tma,qma
2182 real(krealfp),parameter:: terrm=1.e-4
2183 real(krealfp) t,p,tr,pv,pd,el,expo,thet,dthet,terr
2184 integer i
2185 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2186 t=tg
2187 p=pk**con_cpor
2188 do i=1,100
2189 tr=con_ttp/t
2190 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2191 pd=p-pv
2192 el=con_hvap+con_dldt*(t-con_ttp)
2193 expo=el*con_eps*pv/(con_cp*t*pd)
2194 thet=t*pd**(-con_rocp)*exp(expo)
2195 dthet=thet/t*(1.+expo*(con_dldt*t/el+el*p/(con_rv*t*pd)))
2196 terr=(thet-the)/dthet
2197 t=t-terr
2198 if(abs(terr).le.terrm) exit
2199 enddo
2200 tma=t
2201 tr=con_ttp/t
2202 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2203 pd=p-pv
2204 qma=con_eps*pv/(pd+con_eps*pv)
2205 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2206 end subroutine
2207 !-------------------------------------------------------------------------------
2208 subroutine gpkap
2209 !$$$ Subprogram documentation block
2210 !
2211 ! Subprogram: gpkap Compute coefficients for p**kappa
2212 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2213 !
2214 ! Abstract: Computes pressure to the kappa table as a function of pressure
2215 ! for the table lookup function fpkap.
2216 ! Exact pressure to the kappa values are calculated in subprogram fpkapx.
2217 ! The current implementation computes a table with a length
2218 ! of 5501 for pressures ranging up to 110000 Pascals.
2219 !
2220 ! Program History Log:
2221 ! 94-12-30 Iredell
2222 ! 1999-03-01 Iredell f90 module
2223 ! 1999-03-24 Iredell table lookup
2224 !
2225 ! Usage: call gpkap
2226 !
2227 ! Subprograms called:
2228 ! fpkapx function to compute exact pressure to the kappa
2229 !
2230 ! Attributes:
2231 ! Language: Fortran 90.
2232 !
2233 !$$$
2234 implicit none
2235 integer jx
2236 real(krealfp) xmin,xmax,xinc,x,p
2237 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2238 xmin=0._krealfp
2239 xmax=110000._krealfp
2240 xinc=(xmax-xmin)/(nxpkap-1)
2241 c1xpkap=1.-xmin/xinc
2242 c2xpkap=1./xinc
2243 do jx=1,nxpkap
2244 x=xmin+(jx-1)*xinc
2245 p=x
2246 tbpkap(jx)=fpkapx(p)
2247 enddo
2248 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2249 end subroutine
2250 !-------------------------------------------------------------------------------
2251 ! elemental function fpkap(p)
2252 function fpkap(p)
2253 !$$$ Subprogram Documentation Block
2254 !
2255 ! Subprogram: fpkap raise pressure to the kappa power.
2256 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2257 !
2258 ! Abstract: Raise pressure over 1e5 Pa to the kappa power.
2259 ! A linear interpolation is done between values in a lookup table
2260 ! computed in gpkap. See documentation for fpkapx for details.
2261 ! Input values outside table range are reset to table extrema.
2262 ! The interpolation accuracy ranges from 9 decimal places
2263 ! at 100000 Pascals to 5 decimal places at 1000 Pascals.
2264 ! On the Cray, fpkap is over 5 times faster than exact calculation.
2265 ! This function should be expanded inline in the calling routine.
2266 !
2267 ! Program History Log:
2268 ! 91-05-07 Iredell made into inlinable function
2269 ! 94-12-30 Iredell standardized kappa,
2270 ! increased range and accuracy
2271 ! 1999-03-01 Iredell f90 module
2272 ! 1999-03-24 Iredell table lookup
2273 !
2274 ! Usage: pkap=fpkap(p)
2275 !
2276 ! Input argument list:
2277 ! p Real(krealfp) pressure in Pascals
2278 !
2279 ! Output argument list:
2280 ! fpkap Real(krealfp) p over 1e5 Pa to the kappa power
2281 !
2282 ! Attributes:
2283 ! Language: Fortran 90.
2284 !
2285 !$$$
2286 implicit none
2287 real(krealfp) fpkap
2288 real(krealfp),intent(in):: p
2289 integer jx
2290 real(krealfp) xj
2291 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2292 xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
2293 jx=min(xj,nxpkap-1._krealfp)
2294 fpkap=tbpkap(jx)+(xj-jx)*(tbpkap(jx+1)-tbpkap(jx))
2295 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2296 end function
2297 !-------------------------------------------------------------------------------
2298 ! elemental function fpkapq(p)
2299 function fpkapq(p)
2300 !$$$ Subprogram Documentation Block
2301 !
2302 ! Subprogram: fpkapq raise pressure to the kappa power.
2303 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2304 !
2305 ! Abstract: Raise pressure over 1e5 Pa to the kappa power.
2306 ! A quadratic interpolation is done between values in a lookup table
2307 ! computed in gpkap. see documentation for fpkapx for details.
2308 ! Input values outside table range are reset to table extrema.
2309 ! The interpolation accuracy ranges from 12 decimal places
2310 ! at 100000 Pascals to 7 decimal places at 1000 Pascals.
2311 ! On the Cray, fpkap is over 4 times faster than exact calculation.
2312 ! This function should be expanded inline in the calling routine.
2313 !
2314 ! Program History Log:
2315 ! 91-05-07 Iredell made into inlinable function
2316 ! 94-12-30 Iredell standardized kappa,
2317 ! increased range and accuracy
2318 ! 1999-03-01 Iredell f90 module
2319 ! 1999-03-24 Iredell table lookup
2320 !
2321 ! Usage: pkap=fpkapq(p)
2322 !
2323 ! Input argument list:
2324 ! p Real(krealfp) pressure in Pascals
2325 !
2326 ! Output argument list:
2327 ! fpkapq Real(krealfp) p over 1e5 Pa to the kappa power
2328 !
2329 ! Attributes:
2330 ! Language: Fortran 90.
2331 !
2332 !$$$
2333 implicit none
2334 real(krealfp) fpkapq
2335 real(krealfp),intent(in):: p
2336 integer jx
2337 real(krealfp) xj,dxj,fj1,fj2,fj3
2338 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2339 xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
2340 jx=min(max(nint(xj),2),nxpkap-1)
2341 dxj=xj-jx
2342 fj1=tbpkap(jx-1)
2343 fj2=tbpkap(jx)
2344 fj3=tbpkap(jx+1)
2345 fpkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
2346 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2347 end function
2348 !-------------------------------------------------------------------------------
2349 function fpkapo(p)
2350 !$$$ Subprogram documentation block
2351 !
2352 ! Subprogram: fpkapo raise surface pressure to the kappa power.
2353 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2354 !
2355 ! Abstract: Raise surface pressure over 1e5 Pa to the kappa power
2356 ! using a rational weighted chebyshev approximation.
2357 ! The numerator is of order 2 and the denominator is of order 4.
2358 ! The pressure range is 40000-110000 Pa and kappa is defined in fpkapx.
2359 ! The accuracy of this approximation is almost 8 decimal places.
2360 ! On the Cray, fpkap is over 10 times faster than exact calculation.
2361 ! This function should be expanded inline in the calling routine.
2362 !
2363 ! Program History Log:
2364 ! 91-05-07 Iredell made into inlinable function
2365 ! 94-12-30 Iredell standardized kappa,
2366 ! increased range and accuracy
2367 ! 1999-03-01 Iredell f90 module
2368 !
2369 ! Usage: pkap=fpkapo(p)
2370 !
2371 ! Input argument list:
2372 ! p Real(krealfp) surface pressure in Pascals
2373 ! p should be in the range 40000 to 110000
2374 !
2375 ! Output argument list:
2376 ! fpkapo Real(krealfp) p over 1e5 Pa to the kappa power
2377 !
2378 ! Attributes:
2379 ! Language: Fortran 90.
2380 !
2381 !$$$
2382 implicit none
2383 real(krealfp) fpkapo
2384 real(krealfp),intent(in):: p
2385 integer,parameter:: nnpk=2,ndpk=4
2386 real(krealfp):: cnpk(0:nnpk)=(/3.13198449e-1,5.78544829e-2,&
2387 8.35491871e-4/)
2388 real(krealfp):: cdpk(0:ndpk)=(/1.,8.15968401e-2,5.72839518e-4,&
2389 -4.86959812e-7,5.24459889e-10/)
2390 integer n
2391 real(krealfp) pkpa,fnpk,fdpk
2392 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2393 pkpa=p*1.e-3_krealfp
2394 fnpk=cnpk(nnpk)
2395 do n=nnpk-1,0,-1
2396 fnpk=pkpa*fnpk+cnpk(n)
2397 enddo
2398 fdpk=cdpk(ndpk)
2399 do n=ndpk-1,0,-1
2400 fdpk=pkpa*fdpk+cdpk(n)
2401 enddo
2402 fpkapo=fnpk/fdpk
2403 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2404 end function
2405 !-------------------------------------------------------------------------------
2406 ! elemental function fpkapx(p)
2407 function fpkapx(p)
2408 !$$$ Subprogram documentation block
2409 !
2410 ! Subprogram: fpkapx raise pressure to the kappa power.
2411 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2412 !
2413 ! Abstract: raise pressure over 1e5 Pa to the kappa power.
2414 ! Kappa is equal to rd/cp where rd and cp are physical constants.
2415 ! This function should be expanded inline in the calling routine.
2416 !
2417 ! Program History Log:
2418 ! 94-12-30 Iredell made into inlinable function
2419 ! 1999-03-01 Iredell f90 module
2420 !
2421 ! Usage: pkap=fpkapx(p)
2422 !
2423 ! Input argument list:
2424 ! p Real(krealfp) pressure in Pascals
2425 !
2426 ! Output argument list:
2427 ! fpkapx Real(krealfp) p over 1e5 Pa to the kappa power
2428 !
2429 ! Attributes:
2430 ! Language: Fortran 90.
2431 !
2432 !$$$
2433 implicit none
2434 real(krealfp) fpkapx
2435 real(krealfp),intent(in):: p
2436 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2437 fpkapx=(p/1.e5_krealfp)**con_rocp
2438 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2439 end function
2440 !-------------------------------------------------------------------------------
2441 subroutine grkap
2442 !$$$ Subprogram documentation block
2443 !
2444 ! Subprogram: grkap Compute coefficients for p**(1/kappa)
2445 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2446 !
2447 ! Abstract: Computes pressure to the 1/kappa table as a function of pressure
2448 ! for the table lookup function frkap.
2449 ! Exact pressure to the 1/kappa values are calculated in subprogram frkapx.
2450 ! The current implementation computes a table with a length
2451 ! of 5501 for pressures ranging up to 110000 Pascals.
2452 !
2453 ! Program History Log:
2454 ! 94-12-30 Iredell
2455 ! 1999-03-01 Iredell f90 module
2456 ! 1999-03-24 Iredell table lookup
2457 !
2458 ! Usage: call grkap
2459 !
2460 ! Subprograms called:
2461 ! frkapx function to compute exact pressure to the 1/kappa
2462 !
2463 ! Attributes:
2464 ! Language: Fortran 90.
2465 !
2466 !$$$
2467 implicit none
2468 integer jx
2469 real(krealfp) xmin,xmax,xinc,x,p
2470 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2471 xmin=0._krealfp
2472 xmax=fpkapx(110000._krealfp)
2473 xinc=(xmax-xmin)/(nxrkap-1)
2474 c1xrkap=1.-xmin/xinc
2475 c2xrkap=1./xinc
2476 do jx=1,nxrkap
2477 x=xmin+(jx-1)*xinc
2478 p=x
2479 tbrkap(jx)=frkapx(p)
2480 enddo
2481 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2482 end subroutine
2483 !-------------------------------------------------------------------------------
2484 ! elemental function frkap(pkap)
2485 function frkap(pkap)
2486 !$$$ Subprogram Documentation Block
2487 !
2488 ! Subprogram: frkap raise pressure to the 1/kappa power.
2489 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2490 !
2491 ! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
2492 ! A linear interpolation is done between values in a lookup table
2493 ! computed in grkap. See documentation for frkapx for details.
2494 ! Input values outside table range are reset to table extrema.
2495 ! The interpolation accuracy is better than 7 decimal places.
2496 ! On the IBM, fpkap is about 4 times faster than exact calculation.
2497 ! This function should be expanded inline in the calling routine.
2498 !
2499 ! Program History Log:
2500 ! 91-05-07 Iredell made into inlinable function
2501 ! 94-12-30 Iredell standardized kappa,
2502 ! increased range and accuracy
2503 ! 1999-03-01 Iredell f90 module
2504 ! 1999-03-24 Iredell table lookup
2505 !
2506 ! Usage: p=frkap(pkap)
2507 !
2508 ! Input argument list:
2509 ! pkap Real(krealfp) p over 1e5 Pa to the kappa power
2510 !
2511 ! Output argument list:
2512 ! frkap Real(krealfp) pressure in Pascals
2513 !
2514 ! Attributes:
2515 ! Language: Fortran 90.
2516 !
2517 !$$$
2518 implicit none
2519 real(krealfp) frkap
2520 real(krealfp),intent(in):: pkap
2521 integer jx
2522 real(krealfp) xj
2523 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2524 xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
2525 jx=min(xj,nxrkap-1._krealfp)
2526 frkap=tbrkap(jx)+(xj-jx)*(tbrkap(jx+1)-tbrkap(jx))
2527 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2528 end function
2529 !-------------------------------------------------------------------------------
2530 ! elemental function frkapq(pkap)
2531 function frkapq(pkap)
2532 !$$$ Subprogram Documentation Block
2533 !
2534 ! Subprogram: frkapq raise pressure to the 1/kappa power.
2535 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2536 !
2537 ! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
2538 ! A quadratic interpolation is done between values in a lookup table
2539 ! computed in grkap. see documentation for frkapx for details.
2540 ! Input values outside table range are reset to table extrema.
2541 ! The interpolation accuracy is better than 11 decimal places.
2542 ! On the IBM, fpkap is almost 4 times faster than exact calculation.
2543 ! This function should be expanded inline in the calling routine.
2544 !
2545 ! Program History Log:
2546 ! 91-05-07 Iredell made into inlinable function
2547 ! 94-12-30 Iredell standardized kappa,
2548 ! increased range and accuracy
2549 ! 1999-03-01 Iredell f90 module
2550 ! 1999-03-24 Iredell table lookup
2551 !
2552 ! Usage: p=frkapq(pkap)
2553 !
2554 ! Input argument list:
2555 ! pkap Real(krealfp) p over 1e5 Pa to the kappa power
2556 !
2557 ! Output argument list:
2558 ! frkapq Real(krealfp) pressure in Pascals
2559 !
2560 ! Attributes:
2561 ! Language: Fortran 90.
2562 !
2563 !$$$
2564 implicit none
2565 real(krealfp) frkapq
2566 real(krealfp),intent(in):: pkap
2567 integer jx
2568 real(krealfp) xj,dxj,fj1,fj2,fj3
2569 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2570 xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
2571 jx=min(max(nint(xj),2),nxrkap-1)
2572 dxj=xj-jx
2573 fj1=tbrkap(jx-1)
2574 fj2=tbrkap(jx)
2575 fj3=tbrkap(jx+1)
2576 frkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
2577 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2578 end function
2579 !-------------------------------------------------------------------------------
2580 ! elemental function frkapx(pkap)
2581 function frkapx(pkap)
2582 !$$$ Subprogram documentation block
2583 !
2584 ! Subprogram: frkapx raise pressure to the 1/kappa power.
2585 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2586 !
2587 ! Abstract: raise pressure over 1e5 Pa to the 1/kappa power.
2588 ! Kappa is equal to rd/cp where rd and cp are physical constants.
2589 ! This function should be expanded inline in the calling routine.
2590 !
2591 ! Program History Log:
2592 ! 94-12-30 Iredell made into inlinable function
2593 ! 1999-03-01 Iredell f90 module
2594 !
2595 ! Usage: p=frkapx(pkap)
2596 !
2597 ! Input argument list:
2598 ! pkap Real(krealfp) p over 1e5 Pa to the kappa power
2599 !
2600 ! Output argument list:
2601 ! frkapx Real(krealfp) pressure in Pascals
2602 !
2603 ! Attributes:
2604 ! Language: Fortran 90.
2605 !
2606 !$$$
2607 implicit none
2608 real(krealfp) frkapx
2609 real(krealfp),intent(in):: pkap
2610 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2611 frkapx=pkap**(1/con_rocp)*1.e5_krealfp
2612 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2613 end function
2614 !-------------------------------------------------------------------------------
2615 subroutine gtlcl
2616 !$$$ Subprogram Documentation Block
2617 !
2618 ! Subprogram: gtlcl Compute equivalent potential temperature table
2619 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2620 !
2621 ! Abstract: Compute lifting condensation level temperature table
2622 ! as a function of temperature and dewpoint depression for function ftlcl.
2623 ! Lifting condensation level temperature is calculated in subprogram ftlclx
2624 ! The current implementation computes a table with a first dimension
2625 ! of 151 for temperatures ranging from 180.0 to 330.0 Kelvin
2626 ! and a second dimension of 61 for dewpoint depression ranging from
2627 ! 0 to 60 Kelvin.
2628 !
2629 ! Program History Log:
2630 ! 1999-03-01 Iredell f90 module
2631 !
2632 ! Usage: call gtlcl
2633 !
2634 ! Subprograms called:
2635 ! (ftlclx) inlinable function to compute LCL temperature
2636 !
2637 ! Attributes:
2638 ! Language: Fortran 90.
2639 !
2640 !$$$
2641 implicit none
2642 integer jx,jy
2643 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,tdpd,t
2644 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2645 xmin=180._krealfp
2646 xmax=330._krealfp
2647 ymin=0._krealfp
2648 ymax=60._krealfp
2649 xinc=(xmax-xmin)/(nxtlcl-1)
2650 c1xtlcl=1.-xmin/xinc
2651 c2xtlcl=1./xinc
2652 yinc=(ymax-ymin)/(nytlcl-1)
2653 c1ytlcl=1.-ymin/yinc
2654 c2ytlcl=1./yinc
2655 do jy=1,nytlcl
2656 y=ymin+(jy-1)*yinc
2657 tdpd=y
2658 do jx=1,nxtlcl
2659 x=xmin+(jx-1)*xinc
2660 t=x
2661 tbtlcl(jx,jy)=ftlclx(t,tdpd)
2662 enddo
2663 enddo
2664 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2665 end subroutine
2666 !-------------------------------------------------------------------------------
2667 ! elemental function ftlcl(t,tdpd)
2668 function ftlcl(t,tdpd)
2669 !$$$ Subprogram Documentation Block
2670 !
2671 ! Subprogram: ftlcl Compute LCL temperature
2672 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2673 !
2674 ! Abstract: Compute temperature at the lifting condensation level
2675 ! from temperature and dewpoint depression.
2676 ! A bilinear interpolation is done between values in a lookup table
2677 ! computed in gtlcl. See documentation for ftlclx for details.
2678 ! Input values outside table range are reset to table extrema.
2679 ! The interpolation accuracy is better than 0.0005 Kelvin.
2680 ! On the Cray, ftlcl is ? times faster than exact calculation.
2681 ! This function should be expanded inline in the calling routine.
2682 !
2683 ! Program History Log:
2684 ! 1999-03-01 Iredell f90 module
2685 !
2686 ! Usage: tlcl=ftlcl(t,tdpd)
2687 !
2688 ! Input argument list:
2689 ! t Real(krealfp) LCL temperature in Kelvin
2690 ! tdpd Real(krealfp) dewpoint depression in Kelvin
2691 !
2692 ! Output argument list:
2693 ! ftlcl Real(krealfp) temperature at the LCL in Kelvin
2694 !
2695 ! Attributes:
2696 ! Language: Fortran 90.
2697 !
2698 !$$$
2699 implicit none
2700 real(krealfp) ftlcl
2701 real(krealfp),intent(in):: t,tdpd
2702 integer jx,jy
2703 real(krealfp) xj,yj,ftx1,ftx2
2704 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2705 xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
2706 yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
2707 jx=min(xj,nxtlcl-1._krealfp)
2708 jy=min(yj,nytlcl-1._krealfp)
2709 ftx1=tbtlcl(jx,jy)+(xj-jx)*(tbtlcl(jx+1,jy)-tbtlcl(jx,jy))
2710 ftx2=tbtlcl(jx,jy+1)+(xj-jx)*(tbtlcl(jx+1,jy+1)-tbtlcl(jx,jy+1))
2711 ftlcl=ftx1+(yj-jy)*(ftx2-ftx1)
2712 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2713 end function
2714 !-------------------------------------------------------------------------------
2715 ! elemental function ftlclq(t,tdpd)
2716 function ftlclq(t,tdpd)
2717 !$$$ Subprogram Documentation Block
2718 !
2719 ! Subprogram: ftlclq Compute LCL temperature
2720 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2721 !
2722 ! Abstract: Compute temperature at the lifting condensation level
2723 ! from temperature and dewpoint depression.
2724 ! A biquadratic interpolation is done between values in a lookup table
2725 ! computed in gtlcl. see documentation for ftlclx for details.
2726 ! Input values outside table range are reset to table extrema.
2727 ! The interpolation accuracy is better than 0.000003 Kelvin.
2728 ! On the Cray, ftlclq is ? times faster than exact calculation.
2729 ! This function should be expanded inline in the calling routine.
2730 !
2731 ! Program History Log:
2732 ! 1999-03-01 Iredell f90 module
2733 !
2734 ! Usage: tlcl=ftlclq(t,tdpd)
2735 !
2736 ! Input argument list:
2737 ! t Real(krealfp) LCL temperature in Kelvin
2738 ! tdpd Real(krealfp) dewpoint depression in Kelvin
2739 !
2740 ! Output argument list:
2741 ! ftlcl Real(krealfp) temperature at the LCL in Kelvin
2742 !
2743 ! Attributes:
2744 ! Language: Fortran 90.
2745 !
2746 !$$$
2747 implicit none
2748 real(krealfp) ftlclq
2749 real(krealfp),intent(in):: t,tdpd
2750 integer jx,jy
2751 real(krealfp) xj,yj,dxj,dyj
2752 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
2753 real(krealfp) ftx1,ftx2,ftx3
2754 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2755 xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
2756 yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
2757 jx=min(max(nint(xj),2),nxtlcl-1)
2758 jy=min(max(nint(yj),2),nytlcl-1)
2759 dxj=xj-jx
2760 dyj=yj-jy
2761 ft11=tbtlcl(jx-1,jy-1)
2762 ft12=tbtlcl(jx-1,jy)
2763 ft13=tbtlcl(jx-1,jy+1)
2764 ft21=tbtlcl(jx,jy-1)
2765 ft22=tbtlcl(jx,jy)
2766 ft23=tbtlcl(jx,jy+1)
2767 ft31=tbtlcl(jx+1,jy-1)
2768 ft32=tbtlcl(jx+1,jy)
2769 ft33=tbtlcl(jx+1,jy+1)
2770 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
2771 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
2772 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
2773 ftlclq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
2774 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2775 end function
2776 !-------------------------------------------------------------------------------
2777 function ftlclo(t,tdpd)
2778 !$$$ Subprogram documentation block
2779 !
2780 ! Subprogram: ftlclo Compute LCL temperature.
2781 ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2782 !
2783 ! Abstract: Compute temperature at the lifting condensation level
2784 ! from temperature and dewpoint depression. the formula used is
2785 ! a polynomial taken from Phillips mstadb routine which empirically
2786 ! approximates the original exact implicit relationship.
2787 ! (This kind of approximation is customary (inman, 1969), but
2788 ! the original source for this particular one is not yet known. -MI)
2789 ! Its accuracy is about 0.03 Kelvin for a dewpoint depression of 30.
2790 ! This function should be expanded inline in the calling routine.
2791 !
2792 ! Program History Log:
2793 ! 91-05-07 Iredell made into inlinable function
2794 ! 1999-03-01 Iredell f90 module
2795 !
2796 ! Usage: tlcl=ftlclo(t,tdpd)
2797 !
2798 ! Input argument list:
2799 ! t Real(krealfp) temperature in Kelvin
2800 ! tdpd Real(krealfp) dewpoint depression in Kelvin
2801 !
2802 ! Output argument list:
2803 ! ftlclo Real(krealfp) temperature at the LCL in Kelvin
2804 !
2805 ! Attributes:
2806 ! Language: Fortran 90.
2807 !
2808 !$$$
2809 implicit none
2810 real(krealfp) ftlclo
2811 real(krealfp),intent(in):: t,tdpd
2812 real(krealfp),parameter:: clcl1= 0.954442e+0,clcl2= 0.967772e-3,&
2813 clcl3=-0.710321e-3,clcl4=-0.270742e-5
2814 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2815 ftlclo=t-tdpd*(clcl1+clcl2*t+tdpd*(clcl3+clcl4*t))
2816 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2817 end function
2818 !-------------------------------------------------------------------------------
2819 ! elemental function ftlclx(t,tdpd)
2820 function ftlclx(t,tdpd)
2821 !$$$ Subprogram documentation block
2822 !
2823 ! Subprogram: ftlclx Compute LCL temperature.
2824 ! Author: Iredell org: w/NMC2X2 Date: 25 March 1999
2825 !
2826 ! Abstract: Compute temperature at the lifting condensation level
2827 ! from temperature and dewpoint depression. A parcel lifted
2828 ! adiabatically becomes saturated at the lifting condensation level.
2829 ! The water model assumes a perfect gas, constant specific heats
2830 ! for gas and liquid, and neglects the volume of the liquid.
2831 ! The model does account for the variation of the latent heat
2832 ! of condensation with temperature. The ice option is not included.
2833 ! The Clausius-Clapeyron equation is integrated from the triple point
2834 ! to get the formulas
2835 ! pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl))
2836 ! pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew))
2837 ! where pvlcl is the saturated parcel vapor pressure at the LCL,
2838 ! pvdew is the unsaturated parcel vapor pressure initially,
2839 ! trlcl is ttp/tlcl and trdew is ttp/tdew. The adiabatic lifting
2840 ! of the parcel is represented by the following formula
2841 ! pvdew=pvlcl*(t/tlcl)**(1/kappa)
2842 ! This formula is inverted by iterating Newtonian approximations
2843 ! until tlcl is found to within 1.e-6 Kelvin. Note that the minimum
2844 ! returned temperature is 180 Kelvin.
2845 !
2846 ! Program History Log:
2847 ! 1999-03-25 Iredell
2848 !
2849 ! Usage: tlcl=ftlclx(t,tdpd)
2850 !
2851 ! Input argument list:
2852 ! t Real(krealfp) temperature in Kelvin
2853 ! tdpd Real(krealfp) dewpoint depression in Kelvin
2854 !
2855 ! Output argument list:
2856 ! ftlclx Real(krealfp) temperature at the LCL in Kelvin
2857 !
2858 ! Attributes:
2859 ! Language: Fortran 90.
2860 !
2861 !$$$
2862 implicit none
2863 real(krealfp) ftlclx
2864 real(krealfp),intent(in):: t,tdpd
2865 real(krealfp),parameter:: terrm=1.e-4,tlmin=180.,tlminx=tlmin-5.
2866 real(krealfp) tr,pvdew,tlcl,ta,pvlcl,el,dpvlcl,terr,terrp
2867 integer i
2868 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2869 tr=con_ttp/(t-tdpd)
2870 pvdew=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2871 tlcl=t-tdpd
2872 do i=1,100
2873 tr=con_ttp/tlcl
2874 ta=t/tlcl
2875 pvlcl=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))*ta**(1/con_rocp)
2876 el=con_hvap+con_dldt*(tlcl-con_ttp)
2877 dpvlcl=(el/(con_rv*t**2)+1/(con_rocp*tlcl))*pvlcl
2878 terr=(pvlcl-pvdew)/dpvlcl
2879 tlcl=tlcl-terr
2880 if(abs(terr).le.terrm.or.tlcl.lt.tlminx) exit
2881 enddo
2882 ftlclx=max(tlcl,tlmin)
2883 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2884 end function
2885 !-------------------------------------------------------------------------------
2886 subroutine gfuncphys
2887 !$$$ Subprogram Documentation Block
2888 !
2889 ! Subprogram: gfuncphys Compute all physics function tables
2890 ! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2891 !
2892 ! Abstract: Compute all physics function tables. Lookup tables are
2893 ! set up for computing saturation vapor pressure, dewpoint temperature,
2894 ! equivalent potential temperature, moist adiabatic temperature and humidity,
2895 ! pressure to the kappa, and lifting condensation level temperature.
2896 !
2897 ! Program History Log:
2898 ! 1999-03-01 Iredell f90 module
2899 !
2900 ! Usage: call gfuncphys
2901 !
2902 ! Subprograms called:
2903 ! gpvsl compute saturation vapor pressure over liquid table
2904 ! gpvsi compute saturation vapor pressure over ice table
2905 ! gpvs compute saturation vapor pressure table
2906 ! gtdpl compute dewpoint temperature over liquid table
2907 ! gtdpi compute dewpoint temperature over ice table
2908 ! gtdp compute dewpoint temperature table
2909 ! gthe compute equivalent potential temperature table
2910 ! gtma compute moist adiabat tables
2911 ! gpkap compute pressure to the kappa table
2912 ! grkap compute pressure to the 1/kappa table
2913 ! gtlcl compute LCL temperature table
2914 !
2915 ! Attributes:
2916 ! Language: Fortran 90.
2917 !
2918 !$$$
2919 implicit none
2920 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2921 call gpvsl
2922 call gpvsi
2923 call gpvs
2924 call gtdpl
2925 call gtdpi
2926 call gtdp
2927 call gthe
2928 call gtma
2929 call gpkap
2930 call grkap
2931 call gtlcl
2932 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2933 end subroutine
2934 !-------------------------------------------------------------------------------
2935 end module module_gfs_funcphys