CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
funcphys.f90
1
3
7
27!$$$ Module Documentation Block
28!
29! Module: funcphys API for basic thermodynamic physics
30! Author: Iredell Org: W/NX23 Date: 1999-03-01
31!
32! Abstract: This module provides an Application Program Interface
33! for computing basic thermodynamic physics functions, in particular
34! (1) saturation vapor pressure as a function of temperature,
35! (2) dewpoint temperature as a function of vapor pressure,
36! (3) equivalent potential temperature as a function of temperature
37! and scaled pressure to the kappa power,
38! (4) temperature and specific humidity along a moist adiabat
39! as functions of equivalent potential temperature and
40! scaled pressure to the kappa power,
41! (5) scaled pressure to the kappa power as a function of pressure, and
42! (6) temperature at the lifting condensation level as a function
43! of temperature and dewpoint depression.
44! The entry points required to set up lookup tables start with a "g".
45! All the other entry points are functions starting with an "f" or
46! are subroutines starting with an "s". These other functions and
47! subroutines are elemental; that is, they return a scalar if they
48! are passed only scalars, but they return an array if they are passed
49! an array. These other functions and subroutines can be inlined, too.
50!
51! Program History Log:
52! 1999-03-01 Mark Iredell
53! 1999-10-15 Mark Iredell SI unit for pressure (Pascals)
54! 2001-02-26 Mark Iredell Ice phase changes of Hong and Moorthi
55!
56! Public Variables:
57! krealfp Integer parameter kind or length of reals (=kind_phys)
58!
59! Public Subprograms:
60! gpvsl Compute saturation vapor pressure over liquid table
61!
62! fpvsl Elementally compute saturation vapor pressure over liquid
63! function result Real(krealfp) saturation vapor pressure in Pascals
64! t Real(krealfp) temperature in Kelvin
65!
66! fpvslq Elementally compute saturation vapor pressure over liquid
67! function result Real(krealfp) saturation vapor pressure in Pascals
68! t Real(krealfp) temperature in Kelvin
69!
70! fpvslx Elementally compute saturation vapor pressure over liquid
71! function result Real(krealfp) saturation vapor pressure in Pascals
72! t Real(krealfp) temperature in Kelvin
73!
74! gpvsi Compute saturation vapor pressure over ice table
75!
76! fpvsi Elementally compute saturation vapor pressure over ice
77! function result Real(krealfp) saturation vapor pressure in Pascals
78! t Real(krealfp) temperature in Kelvin
79!
80! fpvsiq Elementally compute saturation vapor pressure over ice
81! function result Real(krealfp) saturation vapor pressure in Pascals
82! t Real(krealfp) temperature in Kelvin
83!
84! fpvsix Elementally compute saturation vapor pressure over ice
85! function result Real(krealfp) saturation vapor pressure in Pascals
86! t Real(krealfp) temperature in Kelvin
87!
88! gpvs Compute saturation vapor pressure table
89!
90! fpvs Elementally compute saturation vapor pressure
91! function result Real(krealfp) saturation vapor pressure in Pascals
92! t Real(krealfp) temperature in Kelvin
93!
94! fpvsq Elementally compute saturation vapor pressure
95! function result Real(krealfp) saturation vapor pressure in Pascals
96! t Real(krealfp) temperature in Kelvin
97!
98! fpvsx Elementally compute saturation vapor pressure
99! function result Real(krealfp) saturation vapor pressure in Pascals
100! t Real(krealfp) temperature in Kelvin
101!
102! gtdpl Compute dewpoint temperature over liquid table
103!
104! ftdpl Elementally compute dewpoint temperature over liquid
105! function result Real(krealfp) dewpoint temperature in Kelvin
106! pv Real(krealfp) vapor pressure in Pascals
107!
108! ftdplq Elementally compute dewpoint temperature over liquid
109! function result Real(krealfp) dewpoint temperature in Kelvin
110! pv Real(krealfp) vapor pressure in Pascals
111!
112! ftdplx Elementally compute dewpoint temperature over liquid
113! function result Real(krealfp) dewpoint temperature in Kelvin
114! pv Real(krealfp) vapor pressure in Pascals
115!
116! ftdplxg Elementally compute dewpoint temperature over liquid
117! function result Real(krealfp) dewpoint temperature in Kelvin
118! t Real(krealfp) guess dewpoint temperature in Kelvin
119! pv Real(krealfp) vapor pressure in Pascals
120!
121! gtdpi Compute dewpoint temperature table over ice
122!
123! ftdpi Elementally compute dewpoint temperature over ice
124! function result Real(krealfp) dewpoint temperature in Kelvin
125! pv Real(krealfp) vapor pressure in Pascals
126!
127! ftdpiq Elementally compute dewpoint temperature over ice
128! function result Real(krealfp) dewpoint temperature in Kelvin
129! pv Real(krealfp) vapor pressure in Pascals
130!
131! ftdpix Elementally compute dewpoint temperature over ice
132! function result Real(krealfp) dewpoint temperature in Kelvin
133! pv Real(krealfp) vapor pressure in Pascals
134!
135! ftdpixg Elementally compute dewpoint temperature over ice
136! function result Real(krealfp) dewpoint temperature in Kelvin
137! t Real(krealfp) guess dewpoint temperature in Kelvin
138! pv Real(krealfp) vapor pressure in Pascals
139!
140! gtdp Compute dewpoint temperature table
141!
142! ftdp Elementally compute dewpoint temperature
143! function result Real(krealfp) dewpoint temperature in Kelvin
144! pv Real(krealfp) vapor pressure in Pascals
145!
146! ftdpq Elementally compute dewpoint temperature
147! function result Real(krealfp) dewpoint temperature in Kelvin
148! pv Real(krealfp) vapor pressure in Pascals
149!
150! ftdpx Elementally compute dewpoint temperature
151! function result Real(krealfp) dewpoint temperature in Kelvin
152! pv Real(krealfp) vapor pressure in Pascals
153!
154! ftdpxg Elementally compute dewpoint temperature
155! function result Real(krealfp) dewpoint temperature in Kelvin
156! t Real(krealfp) guess dewpoint temperature in Kelvin
157! pv Real(krealfp) vapor pressure in Pascals
158!
159! gthe Compute equivalent potential temperature table
160!
161! fthe Elementally compute equivalent potential temperature
162! function result Real(krealfp) equivalent potential temperature in Kelvin
163! t Real(krealfp) LCL temperature in Kelvin
164! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
165!
166! ftheq Elementally compute equivalent potential temperature
167! function result Real(krealfp) equivalent potential temperature in Kelvin
168! t Real(krealfp) LCL temperature in Kelvin
169! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
170!
171! fthex Elementally compute equivalent potential temperature
172! function result Real(krealfp) equivalent potential temperature in Kelvin
173! t Real(krealfp) LCL temperature in Kelvin
174! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
175!
176! gtma Compute moist adiabat tables
177!
178! stma Elementally compute moist adiabat temperature and moisture
179! the Real(krealfp) equivalent potential temperature in Kelvin
180! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
181! tma Real(krealfp) parcel temperature in Kelvin
182! qma Real(krealfp) parcel specific humidity in kg/kg
183!
184! stmaq Elementally compute moist adiabat temperature and moisture
185! the Real(krealfp) equivalent potential temperature in Kelvin
186! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
187! tma Real(krealfp) parcel temperature in Kelvin
188! qma Real(krealfp) parcel specific humidity in kg/kg
189!
190! stmax Elementally compute moist adiabat temperature and moisture
191! the Real(krealfp) equivalent potential temperature in Kelvin
192! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
193! tma Real(krealfp) parcel temperature in Kelvin
194! qma Real(krealfp) parcel specific humidity in kg/kg
195!
196! stmaxg Elementally compute moist adiabat temperature and moisture
197! tg Real(krealfp) guess parcel temperature in Kelvin
198! the Real(krealfp) equivalent potential temperature in Kelvin
199! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
200! tma Real(krealfp) parcel temperature in Kelvin
201! qma Real(krealfp) parcel specific humidity in kg/kg
202!
203! gpkap Compute pressure to the kappa table
204!
205! fpkap Elementally raise pressure to the kappa power.
206! function result Real(krealfp) p over 1e5 Pa to the kappa power
207! p Real(krealfp) pressure in Pascals
208!
209! fpkapq Elementally raise pressure to the kappa power.
210! function result Real(krealfp) p over 1e5 Pa to the kappa power
211! p Real(krealfp) pressure in Pascals
212!
213! fpkapo Elementally raise pressure to the kappa power.
214! function result Real(krealfp) p over 1e5 Pa to the kappa power
215! p Real(krealfp) surface pressure in Pascals
216!
217! fpkapx Elementally raise pressure to the kappa power.
218! function result Real(krealfp) p over 1e5 Pa to the kappa power
219! p Real(krealfp) pressure in Pascals
220!
221! grkap Compute pressure to the 1/kappa table
222!
223! frkap Elementally raise pressure to the 1/kappa power.
224! function result Real(krealfp) pressure in Pascals
225! pkap Real(krealfp) p over 1e5 Pa to the 1/kappa power
226!
227! frkapq Elementally raise pressure to the kappa power.
228! function result Real(krealfp) pressure in Pascals
229! pkap Real(krealfp) p over 1e5 Pa to the kappa power
230!
231! frkapx Elementally raise pressure to the kappa power.
232! function result Real(krealfp) pressure in Pascals
233! pkap Real(krealfp) p over 1e5 Pa to the kappa power
234!
235! gtlcl Compute LCL temperature table
236!
237! ftlcl Elementally compute LCL temperature.
238! function result Real(krealfp) temperature at the LCL in Kelvin
239! t Real(krealfp) temperature in Kelvin
240! tdpd Real(krealfp) dewpoint depression in Kelvin
241!
242! ftlclq Elementally compute LCL temperature.
243! function result Real(krealfp) temperature at the LCL in Kelvin
244! t Real(krealfp) temperature in Kelvin
245! tdpd Real(krealfp) dewpoint depression in Kelvin
246!
247! ftlclo Elementally compute LCL temperature.
248! function result Real(krealfp) temperature at the LCL in Kelvin
249! t Real(krealfp) temperature in Kelvin
250! tdpd Real(krealfp) dewpoint depression in Kelvin
251!
252! ftlclx Elementally compute LCL temperature.
253! function result Real(krealfp) temperature at the LCL in Kelvin
254! t Real(krealfp) temperature in Kelvin
255! tdpd Real(krealfp) dewpoint depression in Kelvin
256!
257! gfuncphys Compute all physics function tables
258!
259! Attributes:
260! Language: Fortran 90
261!
262!$$$
263 use machine,only:kind_phys,r8=>kind_dbl_prec,r4=>kind_sngl_prec
264 use physcons
265 implicit none
266 private
267! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268! Public Variables
269! integer,public,parameter:: krealfp=selected_real_kind(15,45)
270 integer,public,parameter:: krealfp=kind_phys
271! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272! Private Variables
273 real(krealfp),parameter:: psatb=con_psat*1.e-5
274 integer,parameter:: nxpvsl=7501
275 real(krealfp) c1xpvsl,c2xpvsl,tbpvsl(nxpvsl)
276 integer,parameter:: nxpvsi=7501
277 real(krealfp) c1xpvsi,c2xpvsi,tbpvsi(nxpvsi)
278 integer,parameter:: nxpvs=7501
279 real(krealfp) c1xpvs,c2xpvs,tbpvs(nxpvs)
280 integer,parameter:: nxtdpl=5001
281 real(krealfp) c1xtdpl,c2xtdpl,tbtdpl(nxtdpl)
282 integer,parameter:: nxtdpi=5001
283 real(krealfp) c1xtdpi,c2xtdpi,tbtdpi(nxtdpi)
284 integer,parameter:: nxtdp=5001
285 real(krealfp) c1xtdp,c2xtdp,tbtdp(nxtdp)
286 integer,parameter:: nxthe=241,nythe=151
287 real(krealfp) c1xthe,c2xthe,c1ythe,c2ythe,tbthe(nxthe,nythe)
288 integer,parameter:: nxma=151,nyma=121
289 real(krealfp) c1xma,c2xma,c1yma,c2yma,tbtma(nxma,nyma),tbqma(nxma,nyma)
290! integer,parameter:: nxpkap=5501
291 integer,parameter:: nxpkap=11001
292 real(krealfp) c1xpkap,c2xpkap,tbpkap(nxpkap)
293 integer,parameter:: nxrkap=5501
294 real(krealfp) c1xrkap,c2xrkap,tbrkap(nxrkap)
295 integer,parameter:: nxtlcl=151,nytlcl=61
296 real(krealfp) c1xtlcl,c2xtlcl,c1ytlcl,c2ytlcl,tbtlcl(nxtlcl,nytlcl)
297! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
298! Public Subprograms
299 public gpvsl,fpvsl,fpvslq,fpvslx
300 public gpvsi,fpvsi,fpvsiq,fpvsix
301 public gpvs,fpvs,fpvsq,fpvsx
302 public gtdpl,ftdpl,ftdplq,ftdplx,ftdplxg
303 public gtdpi,ftdpi,ftdpiq,ftdpix,ftdpixg
304 public gtdp,ftdp,ftdpq,ftdpx,ftdpxg
305 public gthe,fthe,ftheq,fthex
306 public gtma,stma,stmaq,stmax,stmaxg
307 public gpkap,fpkap,fpkapq,fpkapo,fpkapx
308 public grkap,frkap,frkapq,frkapx
309 public gtlcl,ftlcl,ftlclq,ftlclo,ftlclx
310 public gfuncphys
311
312 interface fpvsl
313 module procedure fpvsl_r4, fpvsl_r8
314 end interface fpvsl
315 interface fpvsi
316 module procedure fpvsi_r4, fpvsi_r8
317 end interface fpvsi
318contains
319!-------------------------------------------------------------------------------
325 subroutine gpvsl
326!$$$ Subprogram Documentation Block
327!
328! Subprogram: gpvsl Compute saturation vapor pressure table over liquid
329! Author: N Phillips W/NMC2X2 Date: 30 dec 82
330!
331! Abstract: Computes saturation vapor pressure table as a function of
332! temperature for the table lookup function fpvsl.
333! Exact saturation vapor pressures are calculated in subprogram fpvslx.
334! The current implementation computes a table with a length
335! of 7501 for temperatures ranging from 180. to 330. Kelvin.
336!
337! Program History Log:
338! 91-05-07 Iredell
339! 94-12-30 Iredell expand table
340! 1999-03-01 Iredell f90 module
341!
342! Usage: call gpvsl
343!
344! Subprograms called:
345! (fpvslx) inlinable function to compute saturation vapor pressure
346!
347! Attributes:
348! Language: Fortran 90.
349!
350!$$$
351 implicit none
352 integer jx
353 real(krealfp) xmin,xmax,xinc,x,t
354! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355 xmin=180.0_krealfp
356 xmax=330.0_krealfp
357 xinc=(xmax-xmin)/(nxpvsl-1)
358! c1xpvsl=1.-xmin/xinc
359 c2xpvsl=1./xinc
360 c1xpvsl=1.-xmin*c2xpvsl
361 do jx=1,nxpvsl
362 x=xmin+(jx-1)*xinc
363 t=x
364 tbpvsl(jx)=fpvslx(t)
365 enddo
366! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
367 end subroutine
368!-------------------------------------------------------------------------------
374
375 elemental function fpvsl_r4(t)
376!$$$ Subprogram Documentation Block
377!
378! Subprogram: fpvsl Compute saturation vapor pressure over liquid
379! Author: N Phillips w/NMC2X2 Date: 30 dec 82
380!
381! Abstract: Compute saturation vapor pressure from the temperature.
382! A linear interpolation is done between values in a lookup table
383! computed in gpvsl. See documentation for fpvslx for details.
384! Input values outside table range are reset to table extrema.
385! The interpolation accuracy is almost 6 decimal places.
386! On the Cray, fpvsl is about 4 times faster than exact calculation.
387! This function should be expanded inline in the calling routine.
388!
389! Program History Log:
390! 91-05-07 Iredell made into inlinable function
391! 94-12-30 Iredell expand table
392! 1999-03-01 Iredell f90 module
393!
394! Usage: pvsl=fpvsl(t)
395!
396! Input argument list:
397! t Real(krealfp) temperature in Kelvin
398!
399! Output argument list:
400! fpvsl Real(krealfp) saturation vapor pressure in Pascals
401!
402! Attributes:
403! Language: Fortran 90.
404!
405!$$$
406 implicit none
407 real(r4) fpvsl_r4
408 real(r4),intent(in):: t
409 integer jx
410 real(r4) xj
411! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
412 xj=min(max(c1xpvsl+c2xpvsl*t,1._r4),real(nxpvsl,r4))
413 jx=min(xj,nxpvsl-1._r4)
414 fpvsl_r4=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx))
415! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
416 end function fpvsl_r4
417
418 elemental function fpvsl_r8(t)
419!$$$ Subprogram Documentation Block
420!
421! Subprogram: fpvsl Compute saturation vapor pressure over liquid
422! Author: N Phillips w/NMC2X2 Date: 30 dec 82
423!
424! Abstract: Compute saturation vapor pressure from the temperature.
425! A linear interpolation is done between values in a lookup table
426! computed in gpvsl. See documentation for fpvslx for details.
427! Input values outside table range are reset to table extrema.
428! The interpolation accuracy is almost 6 decimal places.
429! On the Cray, fpvsl is about 4 times faster than exact calculation.
430! This function should be expanded inline in the calling routine.
431!
432! Program History Log:
433! 91-05-07 Iredell made into inlinable function
434! 94-12-30 Iredell expand table
435! 1999-03-01 Iredell f90 module
436!
437! Usage: pvsl=fpvsl(t)
438!
439! Input argument list:
440! t Real(krealfp) temperature in Kelvin
441!
442! Output argument list:
443! fpvsl Real(krealfp) saturation vapor pressure in Pascals
444!
445! Attributes:
446! Language: Fortran 90.
447!
448!$$$
449 implicit none
450 real(r8) fpvsl_r8
451 real(r8),intent(in):: t
452 integer jx
453 real(r8) xj
454! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
455 xj=min(max(c1xpvsl+c2xpvsl*t,1._r8),real(nxpvsl,r8))
456 jx=min(xj,nxpvsl-1._r8)
457 fpvsl_r8=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx))
458! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
459 end function fpvsl_r8
460
461
462
463!-------------------------------------------------------------------------------
468 elemental function fpvslq(t)
469
470!$$$ Subprogram Documentation Block
471!
472! Subprogram: fpvslq Compute saturation vapor pressure over liquid
473! Author: N Phillips w/NMC2X2 Date: 30 dec 82
474!
475! Abstract: Compute saturation vapor pressure from the temperature.
476! A quadratic interpolation is done between values in a lookup table
477! computed in gpvsl. See documentation for fpvslx for details.
478! Input values outside table range are reset to table extrema.
479! The interpolation accuracy is almost 9 decimal places.
480! On the Cray, fpvslq is about 3 times faster than exact calculation.
481! This function should be expanded inline in the calling routine.
482!
483! Program History Log:
484! 91-05-07 Iredell made into inlinable function
485! 94-12-30 Iredell quadratic interpolation
486! 1999-03-01 Iredell f90 module
487!
488! Usage: pvsl=fpvslq(t)
489!
490! Input argument list:
491! t Real(krealfp) temperature in Kelvin
492!
493! Output argument list:
494! fpvslq Real(krealfp) saturation vapor pressure in Pascals
495!
496! Attributes:
497! Language: Fortran 90.
498!
499!$$$
500 implicit none
501 real(krealfp) fpvslq
502 real(krealfp),intent(in):: t
503 integer jx
504 real(krealfp) xj,dxj,fj1,fj2,fj3
505! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506 xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp))
507 jx=min(max(nint(xj),2),nxpvsl-1)
508 dxj=xj-jx
509 fj1=tbpvsl(jx-1)
510 fj2=tbpvsl(jx)
511 fj3=tbpvsl(jx+1)
512 fpvslq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
513! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
514 end function
515!-------------------------------------------------------------------------------
528!\param[out] fpvslx real, saturation vapor pressure in Pascals
529 elemental function fpvslx(t)
530!$$$ Subprogram Documentation Block
531!
532! Subprogram: fpvslx Compute saturation vapor pressure over liquid
533! Author: N Phillips w/NMC2X2 Date: 30 dec 82
534!
535! Abstract: Exactly compute saturation vapor pressure from temperature.
536! The water model assumes a perfect gas, constant specific heats
537! for gas and liquid, and neglects the volume of the liquid.
538! The model does account for the variation of the latent heat
539! of condensation with temperature. The ice option is not included.
540! The Clausius-Clapeyron equation is integrated from the triple point
541! to get the formula
542! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
543! where tr is ttp/t and other values are physical constants.
544! This function should be expanded inline in the calling routine.
545!
546! Program History Log:
547! 91-05-07 Iredell made into inlinable function
548! 94-12-30 Iredell exact computation
549! 1999-03-01 Iredell f90 module
550!
551! Usage: pvsl=fpvslx(t)
552!
553! Input argument list:
554! t Real(krealfp) temperature in Kelvin
555!
556! Output argument list:
557! fpvslx Real(krealfp) saturation vapor pressure in Pascals
558!
559! Attributes:
560! Language: Fortran 90.
561!
562!$$$
563 implicit none
564 real(krealfp) fpvslx
565 real(krealfp),intent(in):: t
566 real(krealfp),parameter:: dldt=con_cvap-con_cliq
567 real(krealfp),parameter:: heat=con_hvap
568 real(krealfp),parameter:: xpona=-dldt/con_rv
569 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
570 real(krealfp) tr
571! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
572 tr=con_ttp/t
573 fpvslx=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
574! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
575 end function
576!-------------------------------------------------------------------------------
583 subroutine gpvsi
584!$$$ Subprogram Documentation Block
585!
586! Subprogram: gpvsi Compute saturation vapor pressure table over ice
587! Author: N Phillips W/NMC2X2 Date: 30 dec 82
588!
589! Abstract: Computes saturation vapor pressure table as a function of
590! temperature for the table lookup function fpvsi.
591! Exact saturation vapor pressures are calculated in subprogram fpvsix.
592! The current implementation computes a table with a length
593! of 7501 for temperatures ranging from 180. to 330. Kelvin.
594!
595! Program History Log:
596! 91-05-07 Iredell
597! 94-12-30 Iredell expand table
598! 1999-03-01 Iredell f90 module
599! 2001-02-26 Iredell ice phase
600!
601! Usage: call gpvsi
602!
603! Subprograms called:
604! (fpvsix) inlinable function to compute saturation vapor pressure
605!
606! Attributes:
607! Language: Fortran 90.
608!
609!$$$
610 implicit none
611 integer jx
612 real(krealfp) xmin,xmax,xinc,x,t
613! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
614 xmin=180.0_krealfp
615 xmax=330.0_krealfp
616 xinc=(xmax-xmin)/(nxpvsi-1)
617! c1xpvsi=1.-xmin/xinc
618 c2xpvsi=1./xinc
619 c1xpvsi=1.-xmin*c2xpvsi
620 do jx=1,nxpvsi
621 x=xmin+(jx-1)*xinc
622 t=x
623 tbpvsi(jx)=fpvsix(t)
624 enddo
625! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
626 end subroutine
627!-------------------------------------------------------------------------------
633
634 elemental function fpvsi_r4(t)
635!$$$ Subprogram Documentation Block
636!
637! Subprogram: fpvsi Compute saturation vapor pressure over ice
638! Author: N Phillips w/NMC2X2 Date: 30 dec 82
639!
640! Abstract: Compute saturation vapor pressure from the temperature.
641! A linear interpolation is done between values in a lookup table
642! computed in gpvsi. See documentation for fpvsix for details.
643! Input values outside table range are reset to table extrema.
644! The interpolation accuracy is almost 6 decimal places.
645! On the Cray, fpvsi is about 4 times faster than exact calculation.
646! This function should be expanded inline in the calling routine.
647!
648! Program History Log:
649! 91-05-07 Iredell made into inlinable function
650! 94-12-30 Iredell expand table
651! 1999-03-01 Iredell f90 module
652! 2001-02-26 Iredell ice phase
653!
654! Usage: pvsi=fpvsi(t)
655!
656! Input argument list:
657! t Real(krealfp) temperature in Kelvin
658!
659! Output argument list:
660! fpvsi Real(krealfp) saturation vapor pressure in Pascals
661!
662! Attributes:
663! Language: Fortran 90.
664!
665!$$$
666 implicit none
667 real(r4) fpvsi_r4
668 real(r4),intent(in):: t
669 integer jx
670 real(r4) xj
671! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
672 xj=min(max(c1xpvsi+c2xpvsi*t,1._r4),real(nxpvsi,r4))
673 jx=min(xj,nxpvsi-1._r4)
674 fpvsi_r4=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx))
675! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
676 end function fpvsi_r4
677
678 elemental function fpvsi_r8(t)
679!$$$ Subprogram Documentation Block
680!
681! Subprogram: fpvsi Compute saturation vapor pressure over ice
682! Author: N Phillips w/NMC2X2 Date: 30 dec 82
683!
684! Abstract: Compute saturation vapor pressure from the temperature.
685! A linear interpolation is done between values in a lookup table
686! computed in gpvsi. See documentation for fpvsix for details.
687! Input values outside table range are reset to table extrema.
688! The interpolation accuracy is almost 6 decimal places.
689! On the Cray, fpvsi is about 4 times faster than exact calculation.
690! This function should be expanded inline in the calling routine.
691!
692! Program History Log:
693! 91-05-07 Iredell made into inlinable function
694! 94-12-30 Iredell expand table
695! 1999-03-01 Iredell f90 module
696! 2001-02-26 Iredell ice phase
697!
698! Usage: pvsi=fpvsi(t)
699!
700! Input argument list:
701! t Real(krealfp) temperature in Kelvin
702!
703! Output argument list:
704! fpvsi Real(krealfp) saturation vapor pressure in Pascals
705!
706! Attributes:
707! Language: Fortran 90.
708!
709!$$$
710 implicit none
711 real(r8) fpvsi_r8
712 real(r8),intent(in):: t
713 integer jx
714 real(r8) xj
715! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
716 xj=min(max(c1xpvsi+c2xpvsi*t,1._r8),real(nxpvsi,r8))
717 jx=min(xj,nxpvsi-1._r8)
718 fpvsi_r8=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx))
719! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
720 end function fpvsi_r8
721
722!-------------------------------------------------------------------------------
728 elemental function fpvsiq(t)
729!$$$ Subprogram Documentation Block
730!
731! Subprogram: fpvsiq Compute saturation vapor pressure over ice
732! Author: N Phillips w/NMC2X2 Date: 30 dec 82
733!
734! Abstract: Compute saturation vapor pressure from the temperature.
735! A quadratic interpolation is done between values in a lookup table
736! computed in gpvsi. See documentation for fpvsix for details.
737! Input values outside table range are reset to table extrema.
738! The interpolation accuracy is almost 9 decimal places.
739! On the Cray, fpvsiq is about 3 times faster than exact calculation.
740! This function should be expanded inline in the calling routine.
741!
742! Program History Log:
743! 91-05-07 Iredell made into inlinable function
744! 94-12-30 Iredell quadratic interpolation
745! 1999-03-01 Iredell f90 module
746! 2001-02-26 Iredell ice phase
747!
748! Usage: pvsi=fpvsiq(t)
749!
750! Input argument list:
751! t Real(krealfp) temperature in Kelvin
752!
753! Output argument list:
754! fpvsiq Real(krealfp) saturation vapor pressure in Pascals
755!
756! Attributes:
757! Language: Fortran 90.
758!
759!$$$
760 implicit none
761 real(krealfp) fpvsiq
762 real(krealfp),intent(in):: t
763 integer jx
764 real(krealfp) xj,dxj,fj1,fj2,fj3
765! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
766 xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp))
767 jx=min(max(nint(xj),2),nxpvsi-1)
768 dxj=xj-jx
769 fj1=tbpvsi(jx-1)
770 fj2=tbpvsi(jx)
771 fj3=tbpvsi(jx+1)
772 fpvsiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
773! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
774 end function
775!-------------------------------------------------------------------------------
786!\param[out] fpvsix real, saturation vapor pressure in Pascals
787 elemental function fpvsix(t)
788!$$$ Subprogram Documentation Block
789!
790! Subprogram: fpvsix Compute saturation vapor pressure over ice
791! Author: N Phillips w/NMC2X2 Date: 30 dec 82
792!
793! Abstract: Exactly compute saturation vapor pressure from temperature.
794! The water model assumes a perfect gas, constant specific heats
795! for gas and ice, and neglects the volume of the ice.
796! The model does account for the variation of the latent heat
797! of condensation with temperature. The liquid option is not included.
798! The Clausius-Clapeyron equation is integrated from the triple point
799! to get the formula
800! pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr))
801! where tr is ttp/t and other values are physical constants.
802! This function should be expanded inline in the calling routine.
803!
804! Program History Log:
805! 91-05-07 Iredell made into inlinable function
806! 94-12-30 Iredell exact computation
807! 1999-03-01 Iredell f90 module
808! 2001-02-26 Iredell ice phase
809!
810! Usage: pvsi=fpvsix(t)
811!
812! Input argument list:
813! t Real(krealfp) temperature in Kelvin
814!
815! Output argument list:
816! fpvsix Real(krealfp) saturation vapor pressure in Pascals
817!
818! Attributes:
819! Language: Fortran 90.
820!
821!$$$
822 implicit none
823 real(krealfp) fpvsix
824 real(krealfp),intent(in):: t
825 real(krealfp),parameter:: dldt=con_cvap-con_csol
826 real(krealfp),parameter:: heat=con_hvap+con_hfus
827 real(krealfp),parameter:: xpona=-dldt/con_rv
828 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
829 real(krealfp) tr
830! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
831 tr=con_ttp/t
832 fpvsix=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
833! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
834 end function
835!-------------------------------------------------------------------------------
841 subroutine gpvs
842!$$$ Subprogram Documentation Block
843!
844! Subprogram: gpvs Compute saturation vapor pressure table
845! Author: N Phillips W/NMC2X2 Date: 30 dec 82
846!
847! Abstract: Computes saturation vapor pressure table as a function of
848! temperature for the table lookup function fpvs.
849! Exact saturation vapor pressures are calculated in subprogram fpvsx.
850! The current implementation computes a table with a length
851! of 7501 for temperatures ranging from 180. to 330. Kelvin.
852!
853! Program History Log:
854! 91-05-07 Iredell
855! 94-12-30 Iredell expand table
856! 1999-03-01 Iredell f90 module
857! 2001-02-26 Iredell ice phase
858!
859! Usage: call gpvs
860!
861! Subprograms called:
862! (fpvsx) inlinable function to compute saturation vapor pressure
863!
864! Attributes:
865! Language: Fortran 90.
866!
867!$$$
868 implicit none
869 integer jx
870 real(krealfp) xmin,xmax,xinc,x,t
871! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
872 xmin=180.0_krealfp
873 xmax=330.0_krealfp
874 xinc=(xmax-xmin)/(nxpvs-1)
875! c1xpvs=1.-xmin/xinc
876 c2xpvs=1./xinc
877 c1xpvs=1.-xmin*c2xpvs
878 do jx=1,nxpvs
879 x=xmin+(jx-1)*xinc
880 t=x
881 tbpvs(jx)=fpvsx(t)
882 enddo
883! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
884 end subroutine
885!-------------------------------------------------------------------------------
891!\param[out] fpvs real, saturation vapor pressure in Pascals
892 elemental function fpvs(t)
893!$$$ Subprogram Documentation Block
894!
895! Subprogram: fpvs Compute saturation vapor pressure
896! Author: N Phillips w/NMC2X2 Date: 30 dec 82
897!
898! Abstract: Compute saturation vapor pressure from the temperature.
899! A linear interpolation is done between values in a lookup table
900! computed in gpvs. See documentation for fpvsx for details.
901! Input values outside table range are reset to table extrema.
902! The interpolation accuracy is almost 6 decimal places.
903! On the Cray, fpvs is about 4 times faster than exact calculation.
904! This function should be expanded inline in the calling routine.
905!
906! Program History Log:
907! 91-05-07 Iredell made into inlinable function
908! 94-12-30 Iredell expand table
909! 1999-03-01 Iredell f90 module
910! 2001-02-26 Iredell ice phase
911!
912! Usage: pvs=fpvs(t)
913!
914! Input argument list:
915! t Real(krealfp) temperature in Kelvin
916!
917! Output argument list:
918! fpvs Real(krealfp) saturation vapor pressure in Pascals
919!
920! Attributes:
921! Language: Fortran 90.
922!
923!$$$
924 implicit none
925 real(krealfp) fpvs
926 real(krealfp),intent(in):: t
927 integer jx
928 real(krealfp) xj
929! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
930 xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
931 jx=min(xj,nxpvs-1._krealfp)
932 fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
933! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
934 end function
935!-------------------------------------------------------------------------------
941!\param[out] fpvsq real, saturation vapor pressure in Pascals
942 elemental function fpvsq(t)
943!$$$ Subprogram Documentation Block
944!
945! Subprogram: fpvsq Compute saturation vapor pressure
946! Author: N Phillips w/NMC2X2 Date: 30 dec 82
947!
948! Abstract: Compute saturation vapor pressure from the temperature.
949! A quadratic interpolation is done between values in a lookup table
950! computed in gpvs. See documentation for fpvsx for details.
951! Input values outside table range are reset to table extrema.
952! The interpolation accuracy is almost 9 decimal places.
953! On the Cray, fpvsq is about 3 times faster than exact calculation.
954! This function should be expanded inline in the calling routine.
955!
956! Program History Log:
957! 91-05-07 Iredell made into inlinable function
958! 94-12-30 Iredell quadratic interpolation
959! 1999-03-01 Iredell f90 module
960! 2001-02-26 Iredell ice phase
961!
962! Usage: pvs=fpvsq(t)
963!
964! Input argument list:
965! t Real(krealfp) temperature in Kelvin
966!
967! Output argument list:
968! fpvsq Real(krealfp) saturation vapor pressure in Pascals
969!
970! Attributes:
971! Language: Fortran 90.
972!
973!$$$
974 implicit none
975 real(krealfp) fpvsq
976 real(krealfp),intent(in):: t
977 integer jx
978 real(krealfp) xj,dxj,fj1,fj2,fj3
979! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
980 xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
981 jx=min(max(nint(xj),2),nxpvs-1)
982 dxj=xj-jx
983 fj1=tbpvs(jx-1)
984 fj2=tbpvs(jx)
985 fj3=tbpvs(jx+1)
986 fpvsq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
987! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
988 end function
989!-------------------------------------------------------------------------------
1006!\param[out] fpvsx real, saturation vapor pressure in Pascals
1007 elemental function fpvsx(t)
1008!$$$ Subprogram Documentation Block
1009!
1010! Subprogram: fpvsx Compute saturation vapor pressure
1011! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1012!
1013! Abstract: Exactly compute saturation vapor pressure from temperature.
1014! The saturation vapor pressure over either liquid and ice is computed
1015! over liquid for temperatures above the triple point,
1016! over ice for temperatures 20 degress below the triple point,
1017! and a linear combination of the two for temperatures in between.
1018! The water model assumes a perfect gas, constant specific heats
1019! for gas, liquid and ice, and neglects the volume of the condensate.
1020! The model does account for the variation of the latent heat
1021! of condensation and sublimation with temperature.
1022! The Clausius-Clapeyron equation is integrated from the triple point
1023! to get the formula
1024! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
1025! where tr is ttp/t and other values are physical constants.
1026! The reference for this computation is Emanuel(1994), pages 116-117.
1027! This function should be expanded inline in the calling routine.
1028!
1029! Program History Log:
1030! 91-05-07 Iredell made into inlinable function
1031! 94-12-30 Iredell exact computation
1032! 1999-03-01 Iredell f90 module
1033! 2001-02-26 Iredell ice phase
1034!
1035! Usage: pvs=fpvsx(t)
1036!
1037! Input argument list:
1038! t Real(krealfp) temperature in Kelvin
1039!
1040! Output argument list:
1041! fpvsx Real(krealfp) saturation vapor pressure in Pascals
1042!
1043! Attributes:
1044! Language: Fortran 90.
1045!
1046!$$$
1047 implicit none
1048 real(krealfp) fpvsx
1049 real(krealfp),intent(in):: t
1050 real(krealfp),parameter:: tliq=con_ttp
1051 real(krealfp),parameter:: tice=con_ttp-20.0
1052 real(krealfp),parameter:: dldtl=con_cvap-con_cliq
1053 real(krealfp),parameter:: heatl=con_hvap
1054 real(krealfp),parameter:: xponal=-dldtl/con_rv
1055 real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
1056 real(krealfp),parameter:: dldti=con_cvap-con_csol
1057 real(krealfp),parameter:: heati=con_hvap+con_hfus
1058 real(krealfp),parameter:: xponai=-dldti/con_rv
1059 real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
1060 real(krealfp) tr,w,pvl,pvi
1061! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1062 tr=con_ttp/t
1063 if(t.ge.tliq) then
1064 fpvsx=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1065 elseif(t.lt.tice) then
1066 fpvsx=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1067 else
1068 w=(t-tice)/(tliq-tice)
1069 pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1070 pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1071 fpvsx=w*pvl+(1.-w)*pvi
1072 endif
1073! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1074 end function
1075!-------------------------------------------------------------------------------
1082 subroutine gtdpl
1083!$$$ Subprogram Documentation Block
1084!
1085! Subprogram: gtdpl Compute dewpoint temperature over liquid table
1086! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1087!
1088! Abstract: Compute dewpoint temperature table as a function of
1089! vapor pressure for inlinable function ftdpl.
1090! Exact dewpoint temperatures are calculated in subprogram ftdplxg.
1091! The current implementation computes a table with a length
1092! of 5001 for vapor pressures ranging from 1 to 10001 Pascals
1093! giving a dewpoint temperature range of 208 to 319 Kelvin.
1094!
1095! Program History Log:
1096! 91-05-07 Iredell
1097! 94-12-30 Iredell expand table
1098! 1999-03-01 Iredell f90 module
1099!
1100! Usage: call gtdpl
1101!
1102! Subprograms called:
1103! (ftdplxg) inlinable function to compute dewpoint temperature over liquid
1104!
1105! Attributes:
1106! Language: Fortran 90.
1107!
1108!$$$
1109 implicit none
1110 integer jx
1111 real(krealfp) xmin,xmax,xinc,t,x,pv
1112! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1113 xmin=1
1114 xmax=10001
1115 xinc=(xmax-xmin)/(nxtdpl-1)
1116 c1xtdpl=1.-xmin/xinc
1117 c2xtdpl=1./xinc
1118 t=208.0
1119 do jx=1,nxtdpl
1120 x=xmin+(jx-1)*xinc
1121 pv=x
1122 t=ftdplxg(t,pv)
1123 tbtdpl(jx)=t
1124 enddo
1125! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1126 end subroutine
1127!-------------------------------------------------------------------------------
1132 elemental function ftdpl(pv)
1133!$$$ Subprogram Documentation Block
1134!
1135! Subprogram: ftdpl Compute dewpoint temperature over liquid
1136! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1137!
1138! Abstract: Compute dewpoint temperature from vapor pressure.
1139! A linear interpolation is done between values in a lookup table
1140! computed in gtdpl. See documentation for ftdplxg for details.
1141! Input values outside table range are reset to table extrema.
1142! The interpolation accuracy is better than 0.0005 Kelvin
1143! for dewpoint temperatures greater than 250 Kelvin,
1144! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
1145! On the Cray, ftdpl is about 75 times faster than exact calculation.
1146! This function should be expanded inline in the calling routine.
1147!
1148! Program History Log:
1149! 91-05-07 Iredell made into inlinable function
1150! 94-12-30 Iredell expand table
1151! 1999-03-01 Iredell f90 module
1152!
1153! Usage: tdpl=ftdpl(pv)
1154!
1155! Input argument list:
1156! pv Real(krealfp) vapor pressure in Pascals
1157!
1158! Output argument list:
1159! ftdpl Real(krealfp) dewpoint temperature in Kelvin
1160!
1161! Attributes:
1162! Language: Fortran 90.
1163!
1164!$$$
1165 implicit none
1166 real(krealfp) ftdpl
1167 real(krealfp),intent(in):: pv
1168 integer jx
1169 real(krealfp) xj
1170! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1171 xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
1172 jx=min(xj,nxtdpl-1._krealfp)
1173 ftdpl=tbtdpl(jx)+(xj-jx)*(tbtdpl(jx+1)-tbtdpl(jx))
1174! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1175 end function
1176!-------------------------------------------------------------------------------
1182!\param[out] ftdplq real, dewpoint temperature in Kelvin
1183 elemental function ftdplq(pv)
1184!$$$ Subprogram Documentation Block
1185!
1186! Subprogram: ftdplq Compute dewpoint temperature over liquid
1187! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1188!
1189! Abstract: Compute dewpoint temperature from vapor pressure.
1190! A quadratic interpolation is done between values in a lookup table
1191! computed in gtdpl. see documentation for ftdplxg for details.
1192! Input values outside table range are reset to table extrema.
1193! the interpolation accuracy is better than 0.00001 Kelvin
1194! for dewpoint temperatures greater than 250 Kelvin,
1195! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
1196! On the Cray, ftdplq is about 60 times faster than exact calculation.
1197! This function should be expanded inline in the calling routine.
1198!
1199! Program History Log:
1200! 91-05-07 Iredell made into inlinable function
1201! 94-12-30 Iredell quadratic interpolation
1202! 1999-03-01 Iredell f90 module
1203!
1204! Usage: tdpl=ftdplq(pv)
1205!
1206! Input argument list:
1207! pv Real(krealfp) vapor pressure in Pascals
1208!
1209! Output argument list:
1210! ftdplq Real(krealfp) dewpoint temperature in Kelvin
1211!
1212! Attributes:
1213! Language: Fortran 90.
1214!
1215!$$$
1216 implicit none
1217 real(krealfp) ftdplq
1218 real(krealfp),intent(in):: pv
1219 integer jx
1220 real(krealfp) xj,dxj,fj1,fj2,fj3
1221! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1222 xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
1223 jx=min(max(nint(xj),2),nxtdpl-1)
1224 dxj=xj-jx
1225 fj1=tbtdpl(jx-1)
1226 fj2=tbtdpl(jx)
1227 fj3=tbtdpl(jx+1)
1228 ftdplq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1229! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1230 end function
1231!-------------------------------------------------------------------------------
1237!\param[out] ftdplx real, dewpoint temperature in Kelvin
1238 elemental function ftdplx(pv)
1239!$$$ Subprogram Documentation Block
1240!
1241! Subprogram: ftdplx Compute dewpoint temperature over liquid
1242! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1243!
1244! Abstract: exactly compute dewpoint temperature from vapor pressure.
1245! An approximate dewpoint temperature for function ftdplxg
1246! is obtained using ftdpl so gtdpl must be already called.
1247! See documentation for ftdplxg for details.
1248!
1249! Program History Log:
1250! 91-05-07 Iredell made into inlinable function
1251! 94-12-30 Iredell exact computation
1252! 1999-03-01 Iredell f90 module
1253!
1254! Usage: tdpl=ftdplx(pv)
1255!
1256! Input argument list:
1257! pv Real(krealfp) vapor pressure in Pascals
1258!
1259! Output argument list:
1260! ftdplx Real(krealfp) dewpoint temperature in Kelvin
1261!
1262! Subprograms called:
1263! (ftdpl) inlinable function to compute dewpoint temperature over liquid
1264! (ftdplxg) inlinable function to compute dewpoint temperature over liquid
1265!
1266! Attributes:
1267! Language: Fortran 90.
1268!
1269!$$$
1270 implicit none
1271 real(krealfp) ftdplx
1272 real(krealfp),intent(in):: pv
1273 real(krealfp) tg
1274! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1275 tg=ftdpl(pv)
1276 ftdplx=ftdplxg(tg,pv)
1277! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1278 end function
1279!-------------------------------------------------------------------------------
1295!\param[out] ftdplxg real, dewpoint temperature in Kelvin
1296 elemental function ftdplxg(tg,pv)
1297!$$$ Subprogram Documentation Block
1298!
1299! Subprogram: ftdplxg Compute dewpoint temperature over liquid
1300! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1301!
1302! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1303! A guess dewpoint temperature must be provided.
1304! The water model assumes a perfect gas, constant specific heats
1305! for gas and liquid, and neglects the volume of the liquid.
1306! The model does account for the variation of the latent heat
1307! of condensation with temperature. The ice option is not included.
1308! The Clausius-Clapeyron equation is integrated from the triple point
1309! to get the formula
1310! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
1311! where tr is ttp/t and other values are physical constants.
1312! The formula is inverted by iterating Newtonian approximations
1313! for each pvs until t is found to within 1.e-6 Kelvin.
1314! This function can be expanded inline in the calling routine.
1315!
1316! Program History Log:
1317! 91-05-07 Iredell made into inlinable function
1318! 94-12-30 Iredell exact computation
1319! 1999-03-01 Iredell f90 module
1320!
1321! Usage: tdpl=ftdplxg(tg,pv)
1322!
1323! Input argument list:
1324! tg Real(krealfp) guess dewpoint temperature in Kelvin
1325! pv Real(krealfp) vapor pressure in Pascals
1326!
1327! Output argument list:
1328! ftdplxg Real(krealfp) dewpoint temperature in Kelvin
1329!
1330! Attributes:
1331! Language: Fortran 90.
1332!
1333!$$$
1334 implicit none
1335 real(krealfp) ftdplxg
1336 real(krealfp),intent(in):: tg,pv
1337 real(krealfp),parameter:: terrm=1.e-6
1338 real(krealfp),parameter:: dldt=con_cvap-con_cliq
1339 real(krealfp),parameter:: heat=con_hvap
1340 real(krealfp),parameter:: xpona=-dldt/con_rv
1341 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
1342 real(krealfp) t,tr,pvt,el,dpvt,terr
1343 integer i
1344! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1345 t=tg
1346 do i=1,100
1347 tr=con_ttp/t
1348 pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
1349 el=heat+dldt*(t-con_ttp)
1350 dpvt=el*pvt/(con_rv*t**2)
1351 terr=(pvt-pv)/dpvt
1352 t=t-terr
1353 if(abs(terr).le.terrm) exit
1354 enddo
1355 ftdplxg=t
1356! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1357 end function
1358!-------------------------------------------------------------------------------
1365 subroutine gtdpi
1366!$$$ Subprogram Documentation Block
1367!
1368! Subprogram: gtdpi Compute dewpoint temperature over ice table
1369! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1370!
1371! Abstract: Compute dewpoint temperature table as a function of
1372! vapor pressure for inlinable function ftdpi.
1373! Exact dewpoint temperatures are calculated in subprogram ftdpixg.
1374! The current implementation computes a table with a length
1375! of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals
1376! giving a dewpoint temperature range of 197 to 279 Kelvin.
1377!
1378! Program History Log:
1379! 91-05-07 Iredell
1380! 94-12-30 Iredell expand table
1381! 1999-03-01 Iredell f90 module
1382! 2001-02-26 Iredell ice phase
1383!
1384! Usage: call gtdpi
1385!
1386! Subprograms called:
1387! (ftdpixg) inlinable function to compute dewpoint temperature over ice
1388!
1389! Attributes:
1390! Language: Fortran 90.
1391!
1392!$$$
1393 implicit none
1394 integer jx
1395 real(krealfp) xmin,xmax,xinc,t,x,pv
1396! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1397 xmin=0.1
1398 xmax=1000.1
1399 xinc=(xmax-xmin)/(nxtdpi-1)
1400 c1xtdpi=1.-xmin/xinc
1401 c2xtdpi=1./xinc
1402 t=197.0
1403 do jx=1,nxtdpi
1404 x=xmin+(jx-1)*xinc
1405 pv=x
1406 t=ftdpixg(t,pv)
1407 tbtdpi(jx)=t
1408 enddo
1409! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1410 end subroutine
1411!-------------------------------------------------------------------------------
1417!\param[out] ftdpi real, dewpoint temperature in Kelvin
1418 elemental function ftdpi(pv)
1419!$$$ Subprogram Documentation Block
1420!
1421! Subprogram: ftdpi Compute dewpoint temperature over ice
1422! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1423!
1424! Abstract: Compute dewpoint temperature from vapor pressure.
1425! A linear interpolation is done between values in a lookup table
1426! computed in gtdpi. See documentation for ftdpixg for details.
1427! Input values outside table range are reset to table extrema.
1428! The interpolation accuracy is better than 0.0005 Kelvin
1429! for dewpoint temperatures greater than 250 Kelvin,
1430! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
1431! On the Cray, ftdpi is about 75 times faster than exact calculation.
1432! This function should be expanded inline in the calling routine.
1433!
1434! Program History Log:
1435! 91-05-07 Iredell made into inlinable function
1436! 94-12-30 Iredell expand table
1437! 1999-03-01 Iredell f90 module
1438! 2001-02-26 Iredell ice phase
1439!
1440! Usage: tdpi=ftdpi(pv)
1441!
1442! Input argument list:
1443! pv Real(krealfp) vapor pressure in Pascals
1444!
1445! Output argument list:
1446! ftdpi Real(krealfp) dewpoint temperature in Kelvin
1447!
1448! Attributes:
1449! Language: Fortran 90.
1450!
1451!$$$
1452 implicit none
1453 real(krealfp) ftdpi
1454 real(krealfp),intent(in):: pv
1455 integer jx
1456 real(krealfp) xj
1457! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1458 xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
1459 jx=min(xj,nxtdpi-1._krealfp)
1460 ftdpi=tbtdpi(jx)+(xj-jx)*(tbtdpi(jx+1)-tbtdpi(jx))
1461! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1462 end function
1463!-------------------------------------------------------------------------------
1469!\param[out] ftdpiq real, dewpoint temperature in Kelvin
1470 elemental function ftdpiq(pv)
1471!$$$ Subprogram Documentation Block
1472!
1473! Subprogram: ftdpiq Compute dewpoint temperature over ice
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 gtdpi. see documentation for ftdpixg 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, ftdpiq 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: tdpi=ftdpiq(pv)
1493!
1494! Input argument list:
1495! pv Real(krealfp) vapor pressure in Pascals
1496!
1497! Output argument list:
1498! ftdpiq Real(krealfp) dewpoint temperature in Kelvin
1499!
1500! Attributes:
1501! Language: Fortran 90.
1502!
1503!$$$
1504 implicit none
1505 real(krealfp) ftdpiq
1506 real(krealfp),intent(in):: pv
1507 integer jx
1508 real(krealfp) xj,dxj,fj1,fj2,fj3
1509! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1510 xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
1511 jx=min(max(nint(xj),2),nxtdpi-1)
1512 dxj=xj-jx
1513 fj1=tbtdpi(jx-1)
1514 fj2=tbtdpi(jx)
1515 fj3=tbtdpi(jx+1)
1516 ftdpiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1517! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1518 end function
1519!-------------------------------------------------------------------------------
1525!\param[out] ftdpix real, dewpoint temperature in Kelvin
1526 elemental function ftdpix(pv)
1527!$$$ Subprogram Documentation Block
1528!
1529! Subprogram: ftdpix Compute dewpoint temperature over ice
1530! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1531!
1532! Abstract: exactly compute dewpoint temperature from vapor pressure.
1533! An approximate dewpoint temperature for function ftdpixg
1534! is obtained using ftdpi so gtdpi must be already called.
1535! See documentation for ftdpixg for details.
1536!
1537! Program History Log:
1538! 91-05-07 Iredell made into inlinable function
1539! 94-12-30 Iredell exact computation
1540! 1999-03-01 Iredell f90 module
1541! 2001-02-26 Iredell ice phase
1542!
1543! Usage: tdpi=ftdpix(pv)
1544!
1545! Input argument list:
1546! pv Real(krealfp) vapor pressure in Pascals
1547!
1548! Output argument list:
1549! ftdpix Real(krealfp) dewpoint temperature in Kelvin
1550!
1551! Subprograms called:
1552! (ftdpi) inlinable function to compute dewpoint temperature over ice
1553! (ftdpixg) inlinable function to compute dewpoint temperature over ice
1554!
1555! Attributes:
1556! Language: Fortran 90.
1557!
1558!$$$
1559 implicit none
1560 real(krealfp) ftdpix
1561 real(krealfp),intent(in):: pv
1562 real(krealfp) tg
1563! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1564 tg=ftdpi(pv)
1565 ftdpix=ftdpixg(tg,pv)
1566! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1567 end function
1568!-------------------------------------------------------------------------------
1584!\param[out] ftdpixg real, dewpoint temperature in Kelvin
1585 elemental function ftdpixg(tg,pv)
1586!$$$ Subprogram Documentation Block
1587!
1588! Subprogram: ftdpixg Compute dewpoint temperature over ice
1589! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1590!
1591! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1592! A guess dewpoint temperature must be provided.
1593! The water model assumes a perfect gas, constant specific heats
1594! for gas and ice, and neglects the volume of the ice.
1595! The model does account for the variation of the latent heat
1596! of sublimation with temperature. The liquid option is not included.
1597! The Clausius-Clapeyron equation is integrated from the triple point
1598! to get the formula
1599! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
1600! where tr is ttp/t and other values are physical constants.
1601! The formula is inverted by iterating Newtonian approximations
1602! for each pvs until t is found to within 1.e-6 Kelvin.
1603! This function can be expanded inline in the calling routine.
1604!
1605! Program History Log:
1606! 91-05-07 Iredell made into inlinable function
1607! 94-12-30 Iredell exact computation
1608! 1999-03-01 Iredell f90 module
1609! 2001-02-26 Iredell ice phase
1610!
1611! Usage: tdpi=ftdpixg(tg,pv)
1612!
1613! Input argument list:
1614! tg Real(krealfp) guess dewpoint temperature in Kelvin
1615! pv Real(krealfp) vapor pressure in Pascals
1616!
1617! Output argument list:
1618! ftdpixg Real(krealfp) dewpoint temperature in Kelvin
1619!
1620! Attributes:
1621! Language: Fortran 90.
1622!
1623!$$$
1624 implicit none
1625 real(krealfp) ftdpixg
1626 real(krealfp),intent(in):: tg,pv
1627 real(krealfp),parameter:: terrm=1.e-6
1628 real(krealfp),parameter:: dldt=con_cvap-con_csol
1629 real(krealfp),parameter:: heat=con_hvap+con_hfus
1630 real(krealfp),parameter:: xpona=-dldt/con_rv
1631 real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
1632 real(krealfp) t,tr,pvt,el,dpvt,terr
1633 integer i
1634! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1635 t=tg
1636 do i=1,100
1637 tr=con_ttp/t
1638 pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
1639 el=heat+dldt*(t-con_ttp)
1640 dpvt=el*pvt/(con_rv*t**2)
1641 terr=(pvt-pv)/dpvt
1642 t=t-terr
1643 if(abs(terr).le.terrm) exit
1644 enddo
1645 ftdpixg=t
1646! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1647 end function
1648!-------------------------------------------------------------------------------
1655 subroutine gtdp
1656!$$$ Subprogram Documentation Block
1657!
1658! Subprogram: gtdp Compute dewpoint temperature table
1659! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1660!
1661! Abstract: Compute dewpoint temperature table as a function of
1662! vapor pressure for inlinable function ftdp.
1663! Exact dewpoint temperatures are calculated in subprogram ftdpxg.
1664! The current implementation computes a table with a length
1665! of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals
1666! giving a dewpoint temperature range of 208 to 319 Kelvin.
1667!
1668! Program History Log:
1669! 91-05-07 Iredell
1670! 94-12-30 Iredell expand table
1671! 1999-03-01 Iredell f90 module
1672! 2001-02-26 Iredell ice phase
1673!
1674! Usage: call gtdp
1675!
1676! Subprograms called:
1677! (ftdpxg) inlinable function to compute dewpoint temperature
1678!
1679! Attributes:
1680! Language: Fortran 90.
1681!
1682!$$$
1683 implicit none
1684 integer jx
1685 real(krealfp) xmin,xmax,xinc,t,x,pv
1686! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1687 xmin=0.5
1688 xmax=10000.5
1689 xinc=(xmax-xmin)/(nxtdp-1)
1690 c1xtdp=1.-xmin/xinc
1691 c2xtdp=1./xinc
1692 t=208.0
1693 do jx=1,nxtdp
1694 x=xmin+(jx-1)*xinc
1695 pv=x
1696 t=ftdpxg(t,pv)
1697 tbtdp(jx)=t
1698 enddo
1699! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1700 end subroutine
1701!-------------------------------------------------------------------------------
1707!\param[out] ftdp real, dewpoint temperature in Kelvin
1708 elemental function ftdp(pv)
1709!$$$ Subprogram Documentation Block
1710!
1711! Subprogram: ftdp Compute dewpoint temperature
1712! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1713!
1714! Abstract: Compute dewpoint temperature from vapor pressure.
1715! A linear interpolation is done between values in a lookup table
1716! computed in gtdp. See documentation for ftdpxg for details.
1717! Input values outside table range are reset to table extrema.
1718! The interpolation accuracy is better than 0.0005 Kelvin
1719! for dewpoint temperatures greater than 250 Kelvin,
1720! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
1721! On the Cray, ftdp is about 75 times faster than exact calculation.
1722! This function should be expanded inline in the calling routine.
1723!
1724! Program History Log:
1725! 91-05-07 Iredell made into inlinable function
1726! 94-12-30 Iredell expand table
1727! 1999-03-01 Iredell f90 module
1728! 2001-02-26 Iredell ice phase
1729!
1730! Usage: tdp=ftdp(pv)
1731!
1732! Input argument list:
1733! pv Real(krealfp) vapor pressure in Pascals
1734!
1735! Output argument list:
1736! ftdp Real(krealfp) dewpoint temperature in Kelvin
1737!
1738! Attributes:
1739! Language: Fortran 90.
1740!
1741!$$$
1742 implicit none
1743 real(krealfp) ftdp
1744 real(krealfp),intent(in):: pv
1745 integer jx
1746 real(krealfp) xj
1747! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1748 xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
1749 jx=min(xj,nxtdp-1._krealfp)
1750 ftdp=tbtdp(jx)+(xj-jx)*(tbtdp(jx+1)-tbtdp(jx))
1751! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1752 end function
1753!-------------------------------------------------------------------------------
1759!\param[out] ftdpq real, dewpoint temperature in Kelvin
1760 elemental function ftdpq(pv)
1761!$$$ Subprogram Documentation Block
1762!
1763! Subprogram: ftdpq Compute dewpoint temperature
1764! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1765!
1766! Abstract: Compute dewpoint temperature from vapor pressure.
1767! A quadratic interpolation is done between values in a lookup table
1768! computed in gtdp. see documentation for ftdpxg for details.
1769! Input values outside table range are reset to table extrema.
1770! the interpolation accuracy is better than 0.00001 Kelvin
1771! for dewpoint temperatures greater than 250 Kelvin,
1772! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
1773! On the Cray, ftdpq is about 60 times faster than exact calculation.
1774! This function should be expanded inline in the calling routine.
1775!
1776! Program History Log:
1777! 91-05-07 Iredell made into inlinable function
1778! 94-12-30 Iredell quadratic interpolation
1779! 1999-03-01 Iredell f90 module
1780! 2001-02-26 Iredell ice phase
1781!
1782! Usage: tdp=ftdpq(pv)
1783!
1784! Input argument list:
1785! pv Real(krealfp) vapor pressure in Pascals
1786!
1787! Output argument list:
1788! ftdpq Real(krealfp) dewpoint temperature in Kelvin
1789!
1790! Attributes:
1791! Language: Fortran 90.
1792!
1793!$$$
1794 implicit none
1795 real(krealfp) ftdpq
1796 real(krealfp),intent(in):: pv
1797 integer jx
1798 real(krealfp) xj,dxj,fj1,fj2,fj3
1799! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1800 xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
1801 jx=min(max(nint(xj),2),nxtdp-1)
1802 dxj=xj-jx
1803 fj1=tbtdp(jx-1)
1804 fj2=tbtdp(jx)
1805 fj3=tbtdp(jx+1)
1806 ftdpq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
1807! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1808 end function
1809!-------------------------------------------------------------------------------
1815!\param[out] ftdpx real, dewpoint temperature in Kelvin
1816 elemental function ftdpx(pv)
1817!$$$ Subprogram Documentation Block
1818!
1819! Subprogram: ftdpx Compute dewpoint temperature
1820! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1821!
1822! Abstract: exactly compute dewpoint temperature from vapor pressure.
1823! An approximate dewpoint temperature for function ftdpxg
1824! is obtained using ftdp so gtdp must be already called.
1825! See documentation for ftdpxg for details.
1826!
1827! Program History Log:
1828! 91-05-07 Iredell made into inlinable function
1829! 94-12-30 Iredell exact computation
1830! 1999-03-01 Iredell f90 module
1831! 2001-02-26 Iredell ice phase
1832!
1833! Usage: tdp=ftdpx(pv)
1834!
1835! Input argument list:
1836! pv Real(krealfp) vapor pressure in Pascals
1837!
1838! Output argument list:
1839! ftdpx Real(krealfp) dewpoint temperature in Kelvin
1840!
1841! Subprograms called:
1842! (ftdp) inlinable function to compute dewpoint temperature
1843! (ftdpxg) inlinable function to compute dewpoint temperature
1844!
1845! Attributes:
1846! Language: Fortran 90.
1847!
1848!$$$
1849 implicit none
1850 real(krealfp) ftdpx
1851 real(krealfp),intent(in):: pv
1852 real(krealfp) tg
1853! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1854 tg=ftdp(pv)
1855 ftdpx=ftdpxg(tg,pv)
1856! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1857 end function
1858!-------------------------------------------------------------------------------
1879!\param[out] ftdpxg real, dewpoint temperature in Kelvin
1880 elemental function ftdpxg(tg,pv)
1881!$$$ Subprogram Documentation Block
1882!
1883! Subprogram: ftdpxg Compute dewpoint temperature
1884! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1885!
1886! Abstract: Exactly compute dewpoint temperature from vapor pressure.
1887! A guess dewpoint temperature must be provided.
1888! The saturation vapor pressure over either liquid and ice is computed
1889! over liquid for temperatures above the triple point,
1890! over ice for temperatures 20 degress below the triple point,
1891! and a linear combination of the two for temperatures in between.
1892! The water model assumes a perfect gas, constant specific heats
1893! for gas, liquid and ice, and neglects the volume of the condensate.
1894! The model does account for the variation of the latent heat
1895! of condensation and sublimation with temperature.
1896! The Clausius-Clapeyron equation is integrated from the triple point
1897! to get the formula
1898! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
1899! where tr is ttp/t and other values are physical constants.
1900! The reference for this decision is Emanuel(1994), pages 116-117.
1901! The formula is inverted by iterating Newtonian approximations
1902! for each pvs until t is found to within 1.e-6 Kelvin.
1903! This function can be expanded inline in the calling routine.
1904!
1905! Program History Log:
1906! 91-05-07 Iredell made into inlinable function
1907! 94-12-30 Iredell exact computation
1908! 1999-03-01 Iredell f90 module
1909! 2001-02-26 Iredell ice phase
1910!
1911! Usage: tdp=ftdpxg(tg,pv)
1912!
1913! Input argument list:
1914! tg Real(krealfp) guess dewpoint temperature in Kelvin
1915! pv Real(krealfp) vapor pressure in Pascals
1916!
1917! Output argument list:
1918! ftdpxg Real(krealfp) dewpoint temperature in Kelvin
1919!
1920! Attributes:
1921! Language: Fortran 90.
1922!
1923!$$$
1924 implicit none
1925 real(krealfp) ftdpxg
1926 real(krealfp),intent(in):: tg,pv
1927 real(krealfp),parameter:: terrm=1.e-6
1928 real(krealfp),parameter:: tliq=con_ttp
1929 real(krealfp),parameter:: tice=con_ttp-20.0
1930 real(krealfp),parameter:: dldtl=con_cvap-con_cliq
1931 real(krealfp),parameter:: heatl=con_hvap
1932 real(krealfp),parameter:: xponal=-dldtl/con_rv
1933 real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
1934 real(krealfp),parameter:: dldti=con_cvap-con_csol
1935 real(krealfp),parameter:: heati=con_hvap+con_hfus
1936 real(krealfp),parameter:: xponai=-dldti/con_rv
1937 real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
1938 real(krealfp) t,tr,w,pvtl,pvti,pvt,ell,eli,el,dpvt,terr
1939 integer i
1940! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1941 t=tg
1942 do i=1,100
1943 tr=con_ttp/t
1944 if(t.ge.tliq) then
1945 pvt=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1946 el=heatl+dldtl*(t-con_ttp)
1947 dpvt=el*pvt/(con_rv*t**2)
1948 elseif(t.lt.tice) then
1949 pvt=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1950 el=heati+dldti*(t-con_ttp)
1951 dpvt=el*pvt/(con_rv*t**2)
1952 else
1953 w=(t-tice)/(tliq-tice)
1954 pvtl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
1955 pvti=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
1956 pvt=w*pvtl+(1.-w)*pvti
1957 ell=heatl+dldtl*(t-con_ttp)
1958 eli=heati+dldti*(t-con_ttp)
1959 dpvt=(w*ell*pvtl+(1.-w)*eli*pvti)/(con_rv*t**2)
1960 endif
1961 terr=(pvt-pv)/dpvt
1962 t=t-terr
1963 if(abs(terr).le.terrm) exit
1964 enddo
1965 ftdpxg=t
1966! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1967 end function
1968!-------------------------------------------------------------------------------
1977 subroutine gthe
1978!$$$ Subprogram Documentation Block
1979!
1980! Subprogram: gthe Compute equivalent potential temperature table
1981! Author: N Phillips w/NMC2X2 Date: 30 dec 82
1982!
1983! Abstract: Compute equivalent potential temperature table
1984! as a function of LCL temperature and pressure over 1e5 Pa
1985! to the kappa power for function fthe.
1986! Equivalent potential temperatures are calculated in subprogram fthex
1987! the current implementation computes a table with a first dimension
1988! of 241 for temperatures ranging from 183.16 to 303.16 Kelvin
1989! and a second dimension of 151 for pressure over 1e5 Pa
1990! to the kappa power ranging from 0.04**rocp to 1.10**rocp.
1991!
1992! Program History Log:
1993! 91-05-07 Iredell
1994! 94-12-30 Iredell expand table
1995! 1999-03-01 Iredell f90 module
1996!
1997! Usage: call gthe
1998!
1999! Subprograms called:
2000! (fthex) inlinable function to compute equiv. pot. temperature
2001!
2002! Attributes:
2003! Language: Fortran 90.
2004!
2005!$$$
2006 implicit none
2007 integer jx,jy
2008 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,t
2009! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2010 xmin=con_ttp-90._krealfp
2011 xmax=con_ttp+30._krealfp
2012 ymin=0.04_krealfp**con_rocp
2013 ymax=1.10_krealfp**con_rocp
2014 xinc=(xmax-xmin)/(nxthe-1)
2015 c1xthe=1.-xmin/xinc
2016 c2xthe=1./xinc
2017 yinc=(ymax-ymin)/(nythe-1)
2018 c1ythe=1.-ymin/yinc
2019 c2ythe=1./yinc
2020 do jy=1,nythe
2021 y=ymin+(jy-1)*yinc
2022 pk=y
2023 do jx=1,nxthe
2024 x=xmin+(jx-1)*xinc
2025 t=x
2026 tbthe(jx,jy)=fthex(t,pk)
2027 enddo
2028 enddo
2029! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2030 end subroutine
2031!-------------------------------------------------------------------------------
2040!\param[out] fthe real, equivalent potential temperature in Kelvin
2041 elemental function fthe(t,pk)
2042!$$$ Subprogram Documentation Block
2043!
2044! Subprogram: fthe Compute equivalent potential temperature
2045! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2046!
2047! Abstract: Compute equivalent potential temperature at the LCL
2048! from temperature and pressure over 1e5 Pa to the kappa power.
2049! A bilinear interpolation is done between values in a lookup table
2050! computed in gthe. see documentation for fthex for details.
2051! Input values outside table range are reset to table extrema,
2052! except zero is returned for too cold or high LCLs.
2053! The interpolation accuracy is better than 0.01 Kelvin.
2054! On the Cray, fthe is almost 6 times faster than exact calculation.
2055! This function should be expanded inline in the calling routine.
2056!
2057! Program History Log:
2058! 91-05-07 Iredell made into inlinable function
2059! 94-12-30 Iredell expand table
2060! 1999-03-01 Iredell f90 module
2061!
2062! Usage: the=fthe(t,pk)
2063!
2064! Input argument list:
2065! t Real(krealfp) LCL temperature in Kelvin
2066! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
2067!
2068! Output argument list:
2069! fthe Real(krealfp) equivalent potential temperature in Kelvin
2070!
2071! Attributes:
2072! Language: Fortran 90.
2073!
2074!$$$
2075 implicit none
2076 real(krealfp) fthe
2077 real(krealfp),intent(in):: t,pk
2078 integer jx,jy
2079 real(krealfp) xj,yj,ftx1,ftx2
2080! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2081 xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
2082 yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
2083 if(xj.ge.1..and.yj.ge.1.) then
2084 jx=min(xj,nxthe-1._krealfp)
2085 jy=min(yj,nythe-1._krealfp)
2086 ftx1=tbthe(jx,jy)+(xj-jx)*(tbthe(jx+1,jy)-tbthe(jx,jy))
2087 ftx2=tbthe(jx,jy+1)+(xj-jx)*(tbthe(jx+1,jy+1)-tbthe(jx,jy+1))
2088 fthe=ftx1+(yj-jy)*(ftx2-ftx1)
2089 else
2090 fthe=0.
2091 endif
2092! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2093 end function
2094!-------------------------------------------------------------------------------
2103!\param[out] ftheq real, equivalent potential temperature in Kelvin
2104 elemental function ftheq(t,pk)
2105!$$$ Subprogram Documentation Block
2106!
2107! Subprogram: ftheq Compute equivalent potential temperature
2108! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2109!
2110! Abstract: Compute equivalent potential temperature at the LCL
2111! from temperature and pressure over 1e5 Pa to the kappa power.
2112! A biquadratic interpolation is done between values in a lookup table
2113! computed in gthe. see documentation for fthex for details.
2114! Input values outside table range are reset to table extrema,
2115! except zero is returned for too cold or high LCLs.
2116! The interpolation accuracy is better than 0.0002 Kelvin.
2117! On the Cray, ftheq is almost 3 times faster than exact calculation.
2118! This function should be expanded inline in the calling routine.
2119!
2120! Program History Log:
2121! 91-05-07 Iredell made into inlinable function
2122! 94-12-30 Iredell quadratic interpolation
2123! 1999-03-01 Iredell f90 module
2124!
2125! Usage: the=ftheq(t,pk)
2126!
2127! Input argument list:
2128! t Real(krealfp) LCL temperature in Kelvin
2129! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
2130!
2131! Output argument list:
2132! ftheq Real(krealfp) equivalent potential temperature in Kelvin
2133!
2134! Attributes:
2135! Language: Fortran 90.
2136!
2137!$$$
2138 implicit none
2139 real(krealfp) ftheq
2140 real(krealfp),intent(in):: t,pk
2141 integer jx,jy
2142 real(krealfp) xj,yj,dxj,dyj
2143 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
2144 real(krealfp) ftx1,ftx2,ftx3
2145! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2146 xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
2147 yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
2148 if(xj.ge.1..and.yj.ge.1.) then
2149 jx=min(max(nint(xj),2),nxthe-1)
2150 jy=min(max(nint(yj),2),nythe-1)
2151 dxj=xj-jx
2152 dyj=yj-jy
2153 ft11=tbthe(jx-1,jy-1)
2154 ft12=tbthe(jx-1,jy)
2155 ft13=tbthe(jx-1,jy+1)
2156 ft21=tbthe(jx,jy-1)
2157 ft22=tbthe(jx,jy)
2158 ft23=tbthe(jx,jy+1)
2159 ft31=tbthe(jx+1,jy-1)
2160 ft32=tbthe(jx+1,jy)
2161 ft33=tbthe(jx+1,jy+1)
2162 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
2163 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
2164 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
2165 ftheq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
2166 else
2167 ftheq=0.
2168 endif
2169! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2170 end function
2171!-------------------------------------------------------------------------------
2172! elemental function fthex(t,pk)
2188!\param[out] fthex real, equivalent potential temperature in Kelvin
2189 function fthex(t,pk)
2190!$$$ Subprogram Documentation Block
2191!
2192! Subprogram: fthex Compute equivalent potential temperature
2193! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2194!
2195! Abstract: Exactly compute equivalent potential temperature at the LCL
2196! from temperature and pressure over 1e5 Pa to the kappa power.
2197! Equivalent potential temperature is constant for a saturated parcel
2198! rising adiabatically up a moist adiabat when the heat and mass
2199! of the condensed water are neglected. Ice is also neglected.
2200! The formula for equivalent potential temperature (Holton) is
2201! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
2202! where t is the temperature, pv is the saturated vapor pressure,
2203! pd is the dry pressure p-pv, el is the temperature dependent
2204! latent heat of condensation hvap+dldt*(t-ttp), and other values
2205! are physical constants defined in parameter statements in the code.
2206! Zero is returned if the input values make saturation impossible.
2207! This function should be expanded inline in the calling routine.
2208!
2209! Program History Log:
2210! 91-05-07 Iredell made into inlinable function
2211! 94-12-30 Iredell exact computation
2212! 1999-03-01 Iredell f90 module
2213!
2214! Usage: the=fthex(t,pk)
2215!
2216! Input argument list:
2217! t Real(krealfp) LCL temperature in Kelvin
2218! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
2219!
2220! Output argument list:
2221! fthex Real(krealfp) equivalent potential temperature in Kelvin
2222!
2223! Attributes:
2224! Language: Fortran 90.
2225!
2226!$$$
2227 implicit none
2228 real(krealfp) fthex
2229 real(krealfp),intent(in):: t,pk
2230 real(krealfp) p,tr,pv,pd,el,expo,expmax
2231! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2232 p=pk**con_cpor
2233 tr=con_ttp/t
2234 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2235 pd=p-pv
2236 if(pd.gt.pv) then
2237 el=con_hvap+con_dldt*(t-con_ttp)
2238 expo=el*con_eps*pv/(con_cp*t*pd)
2239 fthex=t*pd**(-con_rocp)*exp(expo)
2240 else
2241 fthex=0.
2242 endif
2243! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2244 end function
2245!-------------------------------------------------------------------------------
2254 subroutine gtma
2255!$$$ Subprogram Documentation Block
2256!
2257! Subprogram: gtma Compute moist adiabat tables
2258! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2259!
2260! Abstract: Compute temperature and specific humidity tables
2261! as a function of equivalent potential temperature and
2262! pressure over 1e5 Pa to the kappa power for subprogram stma.
2263! Exact parcel temperatures are calculated in subprogram stmaxg.
2264! The current implementation computes a table with a first dimension
2265! of 151 for equivalent potential temperatures ranging from 200 to 500
2266! Kelvin and a second dimension of 121 for pressure over 1e5 Pa
2267! to the kappa power ranging from 0.01**rocp to 1.10**rocp.
2268!
2269! Program History Log:
2270! 91-05-07 Iredell
2271! 94-12-30 Iredell expand table
2272! 1999-03-01 Iredell f90 module
2273!
2274! Usage: call gtma
2275!
2276! Subprograms called:
2277! (stmaxg) inlinable subprogram to compute parcel temperature
2278!
2279! Attributes:
2280! Language: Fortran 90.
2281!
2282!$$$
2283 implicit none
2284 integer jx,jy
2285 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,the,t,q,tg
2286! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2287 xmin=200._krealfp
2288 xmax=500._krealfp
2289 ymin=0.01_krealfp**con_rocp
2290 ymax=1.10_krealfp**con_rocp
2291 xinc=(xmax-xmin)/(nxma-1)
2292 c1xma=1.-xmin/xinc
2293 c2xma=1./xinc
2294 yinc=(ymax-ymin)/(nyma-1)
2295 c1yma=1.-ymin/yinc
2296 c2yma=1./yinc
2297 do jy=1,nyma
2298 y=ymin+(jy-1)*yinc
2299 pk=y
2300 tg=xmin*y
2301 do jx=1,nxma
2302 x=xmin+(jx-1)*xinc
2303 the=x
2304 call stmaxg(tg,the,pk,t,q)
2305 tbtma(jx,jy)=t
2306 tbqma(jx,jy)=q
2307 tg=t
2308 enddo
2309 enddo
2310! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2311 end subroutine
2312!-------------------------------------------------------------------------------
2323 elemental subroutine stma(the,pk,tma,qma)
2324!$$$ Subprogram Documentation Block
2325!
2326! Subprogram: stma Compute moist adiabat temperature
2327! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2328!
2329! Abstract: Compute temperature and specific humidity of a parcel
2330! lifted up a moist adiabat from equivalent potential temperature
2331! at the LCL and pressure over 1e5 Pa to the kappa power.
2332! Bilinear interpolations are done between values in a lookup table
2333! computed in gtma. See documentation for stmaxg for details.
2334! Input values outside table range are reset to table extrema.
2335! The interpolation accuracy is better than 0.01 Kelvin
2336! and 5.e-6 kg/kg for temperature and humidity, respectively.
2337! On the Cray, stma is about 35 times faster than exact calculation.
2338! This subprogram should be expanded inline in the calling routine.
2339!
2340! Program History Log:
2341! 91-05-07 Iredell made into inlinable function
2342! 94-12-30 Iredell expand table
2343! 1999-03-01 Iredell f90 module
2344!
2345! Usage: call stma(the,pk,tma,qma)
2346!
2347! Input argument list:
2348! the Real(krealfp) equivalent potential temperature in Kelvin
2349! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2350!
2351! Output argument list:
2352! tma Real(krealfp) parcel temperature in Kelvin
2353! qma Real(krealfp) parcel specific humidity in kg/kg
2354!
2355! Attributes:
2356! Language: Fortran 90.
2357!
2358!$$$
2359 implicit none
2360 real(krealfp),intent(in):: the,pk
2361 real(krealfp),intent(out):: tma,qma
2362 integer jx,jy
2363 real(krealfp) xj,yj,ftx1,ftx2,qx1,qx2
2364! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2365 xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
2366 yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
2367 jx=min(xj,nxma-1._krealfp)
2368 jy=min(yj,nyma-1._krealfp)
2369 ftx1=tbtma(jx,jy)+(xj-jx)*(tbtma(jx+1,jy)-tbtma(jx,jy))
2370 ftx2=tbtma(jx,jy+1)+(xj-jx)*(tbtma(jx+1,jy+1)-tbtma(jx,jy+1))
2371 tma=ftx1+(yj-jy)*(ftx2-ftx1)
2372 qx1=tbqma(jx,jy)+(xj-jx)*(tbqma(jx+1,jy)-tbqma(jx,jy))
2373 qx2=tbqma(jx,jy+1)+(xj-jx)*(tbqma(jx+1,jy+1)-tbqma(jx,jy+1))
2374 qma=qx1+(yj-jy)*(qx2-qx1)
2375! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2376 end subroutine
2377!-------------------------------------------------------------------------------
2388 elemental subroutine stmaq(the,pk,tma,qma)
2389!$$$ Subprogram Documentation Block
2390!
2391! Subprogram: stmaq Compute moist adiabat temperature
2392! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2393!
2394! Abstract: Compute temperature and specific humidity of a parcel
2395! lifted up a moist adiabat from equivalent potential temperature
2396! at the LCL and pressure over 1e5 Pa to the kappa power.
2397! Biquadratic interpolations are done between values in a lookup table
2398! computed in gtma. See documentation for stmaxg for details.
2399! Input values outside table range are reset to table extrema.
2400! the interpolation accuracy is better than 0.0005 Kelvin
2401! and 1.e-7 kg/kg for temperature and humidity, respectively.
2402! On the Cray, stmaq is about 25 times faster than exact calculation.
2403! This subprogram should be expanded inline in the calling routine.
2404!
2405! Program History Log:
2406! 91-05-07 Iredell made into inlinable function
2407! 94-12-30 Iredell quadratic interpolation
2408! 1999-03-01 Iredell f90 module
2409!
2410! Usage: call stmaq(the,pk,tma,qma)
2411!
2412! Input argument list:
2413! the Real(krealfp) equivalent potential temperature in Kelvin
2414! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2415!
2416! Output argument list:
2417! tmaq Real(krealfp) parcel temperature in Kelvin
2418! qma Real(krealfp) parcel specific humidity in kg/kg
2419!
2420! Attributes:
2421! Language: Fortran 90.
2422!
2423!$$$
2424 implicit none
2425 real(krealfp),intent(in):: the,pk
2426 real(krealfp),intent(out):: tma,qma
2427 integer jx,jy
2428 real(krealfp) xj,yj,dxj,dyj
2429 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
2430 real(krealfp) ftx1,ftx2,ftx3
2431 real(krealfp) q11,q12,q13,q21,q22,q23,q31,q32,q33,qx1,qx2,qx3
2432! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2433 xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
2434 yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
2435 jx=min(max(nint(xj),2),nxma-1)
2436 jy=min(max(nint(yj),2),nyma-1)
2437 dxj=xj-jx
2438 dyj=yj-jy
2439 ft11=tbtma(jx-1,jy-1)
2440 ft12=tbtma(jx-1,jy)
2441 ft13=tbtma(jx-1,jy+1)
2442 ft21=tbtma(jx,jy-1)
2443 ft22=tbtma(jx,jy)
2444 ft23=tbtma(jx,jy+1)
2445 ft31=tbtma(jx+1,jy-1)
2446 ft32=tbtma(jx+1,jy)
2447 ft33=tbtma(jx+1,jy+1)
2448 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
2449 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
2450 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
2451 tma=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
2452 q11=tbqma(jx-1,jy-1)
2453 q12=tbqma(jx-1,jy)
2454 q13=tbqma(jx-1,jy+1)
2455 q21=tbqma(jx,jy-1)
2456 q22=tbqma(jx,jy)
2457 q23=tbqma(jx,jy+1)
2458 q31=tbqma(jx+1,jy-1)
2459 q32=tbqma(jx+1,jy)
2460 q33=tbqma(jx+1,jy+1)
2461 qx1=(((q31+q11)/2-q21)*dxj+(q31-q11)/2)*dxj+q21
2462 qx2=(((q32+q12)/2-q22)*dxj+(q32-q12)/2)*dxj+q22
2463 qx3=(((q33+q13)/2-q23)*dxj+(q33-q13)/2)*dxj+q23
2464 qma=(((qx3+qx1)/2-qx2)*dyj+(qx3-qx1)/2)*dyj+qx2
2465! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2466 end subroutine
2467!-------------------------------------------------------------------------------
2478 subroutine stmax(the,pk,tma,qma)
2479!$$$ Subprogram Documentation Block
2480!
2481! Subprogram: stmax Compute moist adiabat temperature
2482! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2483!
2484! Abstract: Exactly compute temperature and humidity of a parcel
2485! lifted up a moist adiabat from equivalent potential temperature
2486! at the LCL and pressure over 1e5 Pa to the kappa power.
2487! An approximate parcel temperature for subprogram stmaxg
2488! is obtained using stma so gtma must be already called.
2489! See documentation for stmaxg for details.
2490!
2491! Program History Log:
2492! 91-05-07 Iredell made into inlinable function
2493! 94-12-30 Iredell exact computation
2494! 1999-03-01 Iredell f90 module
2495!
2496! Usage: call stmax(the,pk,tma,qma)
2497!
2498! Input argument list:
2499! the Real(krealfp) equivalent potential temperature in Kelvin
2500! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2501!
2502! Output argument list:
2503! tma Real(krealfp) parcel temperature in Kelvin
2504! qma Real(krealfp) parcel specific humidity in kg/kg
2505!
2506! Subprograms called:
2507! (stma) inlinable subprogram to compute parcel temperature
2508! (stmaxg) inlinable subprogram to compute parcel temperature
2509!
2510! Attributes:
2511! Language: Fortran 90.
2512!
2513!$$$
2514 implicit none
2515 real(krealfp),intent(in):: the,pk
2516 real(krealfp),intent(out):: tma,qma
2517 real(krealfp) tg,qg
2518! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2519 call stma(the,pk,tg,qg)
2520 call stmaxg(tg,the,pk,tma,qma)
2521! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2522 end subroutine
2523!-------------------------------------------------------------------------------
2546 subroutine stmaxg(tg,the,pk,tma,qma)
2547!$$$ Subprogram Documentation Block
2548!
2549! Subprogram: stmaxg Compute moist adiabat temperature
2550! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2551!
2552! Abstract: exactly compute temperature and humidity of a parcel
2553! lifted up a moist adiabat from equivalent potential temperature
2554! at the LCL and pressure over 1e5 Pa to the kappa power.
2555! A guess parcel temperature must be provided.
2556! Equivalent potential temperature is constant for a saturated parcel
2557! rising adiabatically up a moist adiabat when the heat and mass
2558! of the condensed water are neglected. Ice is also neglected.
2559! The formula for equivalent potential temperature (Holton) is
2560! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
2561! where t is the temperature, pv is the saturated vapor pressure,
2562! pd is the dry pressure p-pv, el is the temperature dependent
2563! latent heat of condensation hvap+dldt*(t-ttp), and other values
2564! are physical constants defined in parameter statements in the code.
2565! The formula is inverted by iterating Newtonian approximations
2566! for each the and p until t is found to within 1.e-4 Kelvin.
2567! The specific humidity is then computed from pv and pd.
2568! This subprogram can be expanded inline in the calling routine.
2569!
2570! Program History Log:
2571! 91-05-07 Iredell made into inlinable function
2572! 94-12-30 Iredell exact computation
2573! 1999-03-01 Iredell f90 module
2574!
2575! Usage: call stmaxg(tg,the,pk,tma,qma)
2576!
2577! Input argument list:
2578! tg Real(krealfp) guess parcel temperature in Kelvin
2579! the Real(krealfp) equivalent potential temperature in Kelvin
2580! pk Real(krealfp) pressure over 1e5 Pa to the kappa power
2581!
2582! Output argument list:
2583! tma Real(krealfp) parcel temperature in Kelvin
2584! qma Real(krealfp) parcel specific humidity in kg/kg
2585!
2586! Attributes:
2587! Language: Fortran 90.
2588!
2589!$$$
2590 implicit none
2591 real(krealfp),intent(in):: tg,the,pk
2592 real(krealfp),intent(out):: tma,qma
2593 real(krealfp),parameter:: terrm=1.e-4
2594 real(krealfp) t,p,tr,pv,pd,el,expo,thet,dthet,terr
2595 integer i
2596! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2597 t=tg
2598 p=pk**con_cpor
2599 do i=1,100
2600 tr=con_ttp/t
2601 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2602 pd=p-pv
2603 el=con_hvap+con_dldt*(t-con_ttp)
2604 expo=el*con_eps*pv/(con_cp*t*pd)
2605 thet=t*pd**(-con_rocp)*exp(expo)
2606 dthet=thet/t*(1.+expo*(con_dldt*t/el+el*p/(con_rv*t*pd)))
2607 terr=(thet-the)/dthet
2608 t=t-terr
2609 if(abs(terr).le.terrm) exit
2610 enddo
2611 tma=t
2612 tr=con_ttp/t
2613 pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
2614 pd=p-pv
2615 qma=con_eps*pv/(pd+con_eps*pv)
2616! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2617 end subroutine
2618!-------------------------------------------------------------------------------
2624 subroutine gpkap
2625!$$$ Subprogram documentation block
2626!
2627! Subprogram: gpkap Compute coefficients for p**kappa
2628! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2629!
2630! Abstract: Computes pressure to the kappa table as a function of pressure
2631! for the table lookup function fpkap.
2632! Exact pressure to the kappa values are calculated in subprogram fpkapx.
2633! The current implementation computes a table with a length
2634! of 5501 for pressures ranging up to 110000 Pascals.
2635!
2636! Program History Log:
2637! 94-12-30 Iredell
2638! 1999-03-01 Iredell f90 module
2639! 1999-03-24 Iredell table lookup
2640!
2641! Usage: call gpkap
2642!
2643! Subprograms called:
2644! fpkapx function to compute exact pressure to the kappa
2645!
2646! Attributes:
2647! Language: Fortran 90.
2648!
2649!$$$
2650 implicit none
2651 integer jx
2652 real(krealfp) xmin,xmax,xinc,x,p
2653! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2654 xmin=0._krealfp
2655 xmax=110000._krealfp
2656 xinc=(xmax-xmin)/(nxpkap-1)
2657 c1xpkap=1.-xmin/xinc
2658 c2xpkap=1./xinc
2659 do jx=1,nxpkap
2660 x=xmin+(jx-1)*xinc
2661 p=x
2662 tbpkap(jx)=fpkapx(p)
2663 enddo
2664! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2665 end subroutine
2666!-------------------------------------------------------------------------------
2672!\param[out] fpkap real, p over 1e5 Pa to the kappa power
2673 elemental function fpkap(p)
2674!$$$ Subprogram Documentation Block
2675!
2676! Subprogram: fpkap raise pressure to the kappa power.
2677! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2678!
2679! Abstract: Raise pressure over 1e5 Pa to the kappa power.
2680! A linear interpolation is done between values in a lookup table
2681! computed in gpkap. See documentation for fpkapx for details.
2682! Input values outside table range are reset to table extrema.
2683! The interpolation accuracy ranges from 9 decimal places
2684! at 100000 Pascals to 5 decimal places at 1000 Pascals.
2685! On the Cray, fpkap is over 5 times faster than exact calculation.
2686! This function should be expanded inline in the calling routine.
2687!
2688! Program History Log:
2689! 91-05-07 Iredell made into inlinable function
2690! 94-12-30 Iredell standardized kappa,
2691! increased range and accuracy
2692! 1999-03-01 Iredell f90 module
2693! 1999-03-24 Iredell table lookup
2694!
2695! Usage: pkap=fpkap(p)
2696!
2697! Input argument list:
2698! p Real(krealfp) pressure in Pascals
2699!
2700! Output argument list:
2701! fpkap Real(krealfp) p over 1e5 Pa to the kappa power
2702!
2703! Attributes:
2704! Language: Fortran 90.
2705!
2706!$$$
2707 implicit none
2708 real(krealfp) fpkap
2709 real(krealfp),intent(in):: p
2710 integer jx
2711 real(krealfp) xj
2712! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2713 xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
2714 jx=min(xj,nxpkap-1._krealfp)
2715 fpkap=tbpkap(jx)+(xj-jx)*(tbpkap(jx+1)-tbpkap(jx))
2716! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2717 end function
2718!-------------------------------------------------------------------------------
2724!\param[out] fpkapq real, p over 1e5 Pa to the kappa power
2725 elemental function fpkapq(p)
2726!$$$ Subprogram Documentation Block
2727!
2728! Subprogram: fpkapq raise pressure to the kappa power.
2729! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2730!
2731! Abstract: Raise pressure over 1e5 Pa to the kappa power.
2732! A quadratic interpolation is done between values in a lookup table
2733! computed in gpkap. see documentation for fpkapx for details.
2734! Input values outside table range are reset to table extrema.
2735! The interpolation accuracy ranges from 12 decimal places
2736! at 100000 Pascals to 7 decimal places at 1000 Pascals.
2737! On the Cray, fpkap is over 4 times faster than exact calculation.
2738! This function should be expanded inline in the calling routine.
2739!
2740! Program History Log:
2741! 91-05-07 Iredell made into inlinable function
2742! 94-12-30 Iredell standardized kappa,
2743! increased range and accuracy
2744! 1999-03-01 Iredell f90 module
2745! 1999-03-24 Iredell table lookup
2746!
2747! Usage: pkap=fpkapq(p)
2748!
2749! Input argument list:
2750! p Real(krealfp) pressure in Pascals
2751!
2752! Output argument list:
2753! fpkapq Real(krealfp) p over 1e5 Pa to the kappa power
2754!
2755! Attributes:
2756! Language: Fortran 90.
2757!
2758!$$$
2759 implicit none
2760 real(krealfp) fpkapq
2761 real(krealfp),intent(in):: p
2762 integer jx
2763 real(krealfp) xj,dxj,fj1,fj2,fj3
2764! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2765 xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
2766 jx=min(max(nint(xj),2),nxpkap-1)
2767 dxj=xj-jx
2768 fj1=tbpkap(jx-1)
2769 fj2=tbpkap(jx)
2770 fj3=tbpkap(jx+1)
2771 fpkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
2772! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2773 end function
2774!-------------------------------------------------------------------------------
2781!\param[out] fpkapo real, p over 1e5 Pa to the kappa power
2782 function fpkapo(p)
2783!$$$ Subprogram documentation block
2784!
2785! Subprogram: fpkapo raise surface pressure to the kappa power.
2786! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2787!
2788! Abstract: Raise surface pressure over 1e5 Pa to the kappa power
2789! using a rational weighted chebyshev approximation.
2790! The numerator is of order 2 and the denominator is of order 4.
2791! The pressure range is 40000-110000 Pa and kappa is defined in fpkapx.
2792! The accuracy of this approximation is almost 8 decimal places.
2793! On the Cray, fpkap is over 10 times faster than exact calculation.
2794! This function should be expanded inline in the calling routine.
2795!
2796! Program History Log:
2797! 91-05-07 Iredell made into inlinable function
2798! 94-12-30 Iredell standardized kappa,
2799! increased range and accuracy
2800! 1999-03-01 Iredell f90 module
2801!
2802! Usage: pkap=fpkapo(p)
2803!
2804! Input argument list:
2805! p Real(krealfp) surface pressure in Pascals
2806! p should be in the range 40000 to 110000
2807!
2808! Output argument list:
2809! fpkapo Real(krealfp) p over 1e5 Pa to the kappa power
2810!
2811! Attributes:
2812! Language: Fortran 90.
2813!
2814!$$$
2815 implicit none
2816 real(krealfp) fpkapo
2817 real(krealfp),intent(in):: p
2818 integer,parameter:: nnpk=2,ndpk=4
2819 real(krealfp):: cnpk(0:nnpk)=(/3.13198449e-1,5.78544829e-2,&
2820 8.35491871e-4/)
2821 real(krealfp):: cdpk(0:ndpk)=(/1.,8.15968401e-2,5.72839518e-4,&
2822 -4.86959812e-7,5.24459889e-10/)
2823 integer n
2824 real(krealfp) pkpa,fnpk,fdpk
2825! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2826 pkpa=p*1.e-3_krealfp
2827 fnpk=cnpk(nnpk)
2828 do n=nnpk-1,0,-1
2829 fnpk=pkpa*fnpk+cnpk(n)
2830 enddo
2831 fdpk=cdpk(ndpk)
2832 do n=ndpk-1,0,-1
2833 fdpk=pkpa*fdpk+cdpk(n)
2834 enddo
2835 fpkapo=fnpk/fdpk
2836! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2837 end function
2838!-------------------------------------------------------------------------------
2842!\param[out] fpkapx real, p over 1e5 Pa to the kappa power
2843 elemental function fpkapx(p)
2844!$$$ Subprogram documentation block
2845!
2846! Subprogram: fpkapx raise pressure to the kappa power.
2847! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2848!
2849! Abstract: raise pressure over 1e5 Pa to the kappa power.
2850! Kappa is equal to rd/cp where rd and cp are physical constants.
2851! This function should be expanded inline in the calling routine.
2852!
2853! Program History Log:
2854! 94-12-30 Iredell made into inlinable function
2855! 1999-03-01 Iredell f90 module
2856!
2857! Usage: pkap=fpkapx(p)
2858!
2859! Input argument list:
2860! p Real(krealfp) pressure in Pascals
2861!
2862! Output argument list:
2863! fpkapx Real(krealfp) p over 1e5 Pa to the kappa power
2864!
2865! Attributes:
2866! Language: Fortran 90.
2867!
2868!$$$
2869 implicit none
2870 real(krealfp) fpkapx
2871 real(krealfp),intent(in):: p
2872! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2873 fpkapx=(p/1.e5_krealfp)**con_rocp
2874! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2875 end function
2876!-------------------------------------------------------------------------------
2882 subroutine grkap
2883!$$$ Subprogram documentation block
2884!
2885! Subprogram: grkap Compute coefficients for p**(1/kappa)
2886! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
2887!
2888! Abstract: Computes pressure to the 1/kappa table as a function of pressure
2889! for the table lookup function frkap.
2890! Exact pressure to the 1/kappa values are calculated in subprogram frkapx.
2891! The current implementation computes a table with a length
2892! of 5501 for pressures ranging up to 110000 Pascals.
2893!
2894! Program History Log:
2895! 94-12-30 Iredell
2896! 1999-03-01 Iredell f90 module
2897! 1999-03-24 Iredell table lookup
2898!
2899! Usage: call grkap
2900!
2901! Subprograms called:
2902! frkapx function to compute exact pressure to the 1/kappa
2903!
2904! Attributes:
2905! Language: Fortran 90.
2906!
2907!$$$
2908 implicit none
2909 integer jx
2910 real(krealfp) xmin,xmax,xinc,x,p
2911! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2912 xmin=0._krealfp
2913 xmax=fpkapx(110000._krealfp)
2914 xinc=(xmax-xmin)/(nxrkap-1)
2915 c1xrkap=1.-xmin/xinc
2916 c2xrkap=1./xinc
2917 do jx=1,nxrkap
2918 x=xmin+(jx-1)*xinc
2919 p=x
2920 tbrkap(jx)=frkapx(p)
2921 enddo
2922! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2923 end subroutine
2924!-------------------------------------------------------------------------------
2930!\param[out] frkap real, pressure in Pascals
2931 elemental function frkap(pkap)
2932!$$$ Subprogram Documentation Block
2933!
2934! Subprogram: frkap raise pressure to the 1/kappa power.
2935! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2936!
2937! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
2938! A linear interpolation is done between values in a lookup table
2939! computed in grkap. See documentation for frkapx for details.
2940! Input values outside table range are reset to table extrema.
2941! The interpolation accuracy is better than 7 decimal places.
2942! On the IBM, fpkap is about 4 times faster than exact calculation.
2943! This function should be expanded inline in the calling routine.
2944!
2945! Program History Log:
2946! 91-05-07 Iredell made into inlinable function
2947! 94-12-30 Iredell standardized kappa,
2948! increased range and accuracy
2949! 1999-03-01 Iredell f90 module
2950! 1999-03-24 Iredell table lookup
2951!
2952! Usage: p=frkap(pkap)
2953!
2954! Input argument list:
2955! pkap Real(krealfp) p over 1e5 Pa to the kappa power
2956!
2957! Output argument list:
2958! frkap Real(krealfp) pressure in Pascals
2959!
2960! Attributes:
2961! Language: Fortran 90.
2962!
2963!$$$
2964 implicit none
2965 real(krealfp) frkap
2966 real(krealfp),intent(in):: pkap
2967 integer jx
2968 real(krealfp) xj
2969! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2970 xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
2971 jx=min(xj,nxrkap-1._krealfp)
2972 frkap=tbrkap(jx)+(xj-jx)*(tbrkap(jx+1)-tbrkap(jx))
2973! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2974 end function
2975!-------------------------------------------------------------------------------
2981!\param[out] frkapq real, pressure in Pascals
2982 elemental function frkapq(pkap)
2983!$$$ Subprogram Documentation Block
2984!
2985! Subprogram: frkapq raise pressure to the 1/kappa power.
2986! Author: N Phillips w/NMC2X2 Date: 30 dec 82
2987!
2988! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
2989! A quadratic interpolation is done between values in a lookup table
2990! computed in grkap. see documentation for frkapx for details.
2991! Input values outside table range are reset to table extrema.
2992! The interpolation accuracy is better than 11 decimal places.
2993! On the IBM, fpkap is almost 4 times faster than exact calculation.
2994! This function should be expanded inline in the calling routine.
2995!
2996! Program History Log:
2997! 91-05-07 Iredell made into inlinable function
2998! 94-12-30 Iredell standardized kappa,
2999! increased range and accuracy
3000! 1999-03-01 Iredell f90 module
3001! 1999-03-24 Iredell table lookup
3002!
3003! Usage: p=frkapq(pkap)
3004!
3005! Input argument list:
3006! pkap Real(krealfp) p over 1e5 Pa to the kappa power
3007!
3008! Output argument list:
3009! frkapq Real(krealfp) pressure in Pascals
3010!
3011! Attributes:
3012! Language: Fortran 90.
3013!
3014!$$$
3015 implicit none
3016 real(krealfp) frkapq
3017 real(krealfp),intent(in):: pkap
3018 integer jx
3019 real(krealfp) xj,dxj,fj1,fj2,fj3
3020! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3021 xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
3022 jx=min(max(nint(xj),2),nxrkap-1)
3023 dxj=xj-jx
3024 fj1=tbrkap(jx-1)
3025 fj2=tbrkap(jx)
3026 fj3=tbrkap(jx+1)
3027 frkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
3028! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3029 end function
3030!-------------------------------------------------------------------------------
3034!\param[out] frkapx real, pressure in Pascals
3035 elemental function frkapx(pkap)
3036!$$$ Subprogram documentation block
3037!
3038! Subprogram: frkapx raise pressure to the 1/kappa power.
3039! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
3040!
3041! Abstract: raise pressure over 1e5 Pa to the 1/kappa power.
3042! Kappa is equal to rd/cp where rd and cp are physical constants.
3043! This function should be expanded inline in the calling routine.
3044!
3045! Program History Log:
3046! 94-12-30 Iredell made into inlinable function
3047! 1999-03-01 Iredell f90 module
3048!
3049! Usage: p=frkapx(pkap)
3050!
3051! Input argument list:
3052! pkap Real(krealfp) p over 1e5 Pa to the kappa power
3053!
3054! Output argument list:
3055! frkapx Real(krealfp) pressure in Pascals
3056!
3057! Attributes:
3058! Language: Fortran 90.
3059!
3060!$$$
3061 implicit none
3062 real(krealfp) frkapx
3063 real(krealfp),intent(in):: pkap
3064! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3065 frkapx=pkap**(1/con_rocp)*1.e5_krealfp
3066! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3067 end function
3068!-------------------------------------------------------------------------------
3076 subroutine gtlcl
3077!$$$ Subprogram Documentation Block
3078!
3079! Subprogram: gtlcl Compute equivalent potential temperature table
3080! Author: N Phillips w/NMC2X2 Date: 30 dec 82
3081!
3082! Abstract: Compute lifting condensation level temperature table
3083! as a function of temperature and dewpoint depression for function ftlcl.
3084! Lifting condensation level temperature is calculated in subprogram ftlclx
3085! The current implementation computes a table with a first dimension
3086! of 151 for temperatures ranging from 180.0 to 330.0 Kelvin
3087! and a second dimension of 61 for dewpoint depression ranging from
3088! 0 to 60 Kelvin.
3089!
3090! Program History Log:
3091! 1999-03-01 Iredell f90 module
3092!
3093! Usage: call gtlcl
3094!
3095! Subprograms called:
3096! (ftlclx) inlinable function to compute LCL temperature
3097!
3098! Attributes:
3099! Language: Fortran 90.
3100!
3101!$$$
3102 implicit none
3103 integer jx,jy
3104 real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,tdpd,t
3105! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3106 xmin=180._krealfp
3107 xmax=330._krealfp
3108 ymin=0._krealfp
3109 ymax=60._krealfp
3110 xinc=(xmax-xmin)/(nxtlcl-1)
3111 c1xtlcl=1.-xmin/xinc
3112 c2xtlcl=1./xinc
3113 yinc=(ymax-ymin)/(nytlcl-1)
3114 c1ytlcl=1.-ymin/yinc
3115 c2ytlcl=1./yinc
3116 do jy=1,nytlcl
3117 y=ymin+(jy-1)*yinc
3118 tdpd=y
3119 do jx=1,nxtlcl
3120 x=xmin+(jx-1)*xinc
3121 t=x
3122 tbtlcl(jx,jy)=ftlclx(t,tdpd)
3123 enddo
3124 enddo
3125! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3126 end subroutine
3127!-------------------------------------------------------------------------------
3135!\param[out] ftlcl real, temperature at the LCL in Kelvin
3136 elemental function ftlcl(t,tdpd)
3137!$$$ Subprogram Documentation Block
3138!
3139! Subprogram: ftlcl Compute LCL temperature
3140! Author: N Phillips w/NMC2X2 Date: 30 dec 82
3141!
3142! Abstract: Compute temperature at the lifting condensation level
3143! from temperature and dewpoint depression.
3144! A bilinear interpolation is done between values in a lookup table
3145! computed in gtlcl. See documentation for ftlclx for details.
3146! Input values outside table range are reset to table extrema.
3147! The interpolation accuracy is better than 0.0005 Kelvin.
3148! On the Cray, ftlcl is ? times faster than exact calculation.
3149! This function should be expanded inline in the calling routine.
3150!
3151! Program History Log:
3152! 1999-03-01 Iredell f90 module
3153!
3154! Usage: tlcl=ftlcl(t,tdpd)
3155!
3156! Input argument list:
3157! t Real(krealfp) LCL temperature in Kelvin
3158! tdpd Real(krealfp) dewpoint depression in Kelvin
3159!
3160! Output argument list:
3161! ftlcl Real(krealfp) temperature at the LCL in Kelvin
3162!
3163! Attributes:
3164! Language: Fortran 90.
3165!
3166!$$$
3167 implicit none
3168 real(krealfp) ftlcl
3169 real(krealfp),intent(in):: t,tdpd
3170 integer jx,jy
3171 real(krealfp) xj,yj,ftx1,ftx2
3172! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3173 xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
3174 yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
3175 jx=min(xj,nxtlcl-1._krealfp)
3176 jy=min(yj,nytlcl-1._krealfp)
3177 ftx1=tbtlcl(jx,jy)+(xj-jx)*(tbtlcl(jx+1,jy)-tbtlcl(jx,jy))
3178 ftx2=tbtlcl(jx,jy+1)+(xj-jx)*(tbtlcl(jx+1,jy+1)-tbtlcl(jx,jy+1))
3179 ftlcl=ftx1+(yj-jy)*(ftx2-ftx1)
3180! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3181 end function
3182!-------------------------------------------------------------------------------
3190!\param[out] ftlcl real, temperature at the LCL in Kelvin
3191 elemental function ftlclq(t,tdpd)
3192!$$$ Subprogram Documentation Block
3193!
3194! Subprogram: ftlclq Compute LCL temperature
3195! Author: N Phillips w/NMC2X2 Date: 30 dec 82
3196!
3197! Abstract: Compute temperature at the lifting condensation level
3198! from temperature and dewpoint depression.
3199! A biquadratic interpolation is done between values in a lookup table
3200! computed in gtlcl. see documentation for ftlclx for details.
3201! Input values outside table range are reset to table extrema.
3202! The interpolation accuracy is better than 0.000003 Kelvin.
3203! On the Cray, ftlclq is ? times faster than exact calculation.
3204! This function should be expanded inline in the calling routine.
3205!
3206! Program History Log:
3207! 1999-03-01 Iredell f90 module
3208!
3209! Usage: tlcl=ftlclq(t,tdpd)
3210!
3211! Input argument list:
3212! t Real(krealfp) LCL temperature in Kelvin
3213! tdpd Real(krealfp) dewpoint depression in Kelvin
3214!
3215! Output argument list:
3216! ftlcl Real(krealfp) temperature at the LCL in Kelvin
3217!
3218! Attributes:
3219! Language: Fortran 90.
3220!
3221!$$$
3222 implicit none
3223 real(krealfp) ftlclq
3224 real(krealfp),intent(in):: t,tdpd
3225 integer jx,jy
3226 real(krealfp) xj,yj,dxj,dyj
3227 real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
3228 real(krealfp) ftx1,ftx2,ftx3
3229! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3230 xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
3231 yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
3232 jx=min(max(nint(xj),2),nxtlcl-1)
3233 jy=min(max(nint(yj),2),nytlcl-1)
3234 dxj=xj-jx
3235 dyj=yj-jy
3236 ft11=tbtlcl(jx-1,jy-1)
3237 ft12=tbtlcl(jx-1,jy)
3238 ft13=tbtlcl(jx-1,jy+1)
3239 ft21=tbtlcl(jx,jy-1)
3240 ft22=tbtlcl(jx,jy)
3241 ft23=tbtlcl(jx,jy+1)
3242 ft31=tbtlcl(jx+1,jy-1)
3243 ft32=tbtlcl(jx+1,jy)
3244 ft33=tbtlcl(jx+1,jy+1)
3245 ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
3246 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
3247 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
3248 ftlclq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
3249! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3250 end function
3251!-------------------------------------------------------------------------------
3258!\param[out] ftlclo real, temperature at the LCL in Kelvin
3259 function ftlclo(t,tdpd)
3260!$$$ Subprogram documentation block
3261!
3262! Subprogram: ftlclo Compute LCL temperature.
3263! Author: Phillips org: w/NMC2X2 Date: 29 dec 82
3264!
3265! Abstract: Compute temperature at the lifting condensation level
3266! from temperature and dewpoint depression. the formula used is
3267! a polynomial taken from Phillips mstadb routine which empirically
3268! approximates the original exact implicit relationship.
3269! (This kind of approximation is customary (inman, 1969), but
3270! the original source for this particular one is not yet known. -MI)
3271! Its accuracy is about 0.03 Kelvin for a dewpoint depression of 30.
3272! This function should be expanded inline in the calling routine.
3273!
3274! Program History Log:
3275! 91-05-07 Iredell made into inlinable function
3276! 1999-03-01 Iredell f90 module
3277!
3278! Usage: tlcl=ftlclo(t,tdpd)
3279!
3280! Input argument list:
3281! t Real(krealfp) temperature in Kelvin
3282! tdpd Real(krealfp) dewpoint depression in Kelvin
3283!
3284! Output argument list:
3285! ftlclo Real(krealfp) temperature at the LCL in Kelvin
3286!
3287! Attributes:
3288! Language: Fortran 90.
3289!
3290!$$$
3291 implicit none
3292 real(krealfp) ftlclo
3293 real(krealfp),intent(in):: t,tdpd
3294 real(krealfp),parameter:: clcl1= 0.954442e+0,clcl2= 0.967772e-3,&
3295 clcl3=-0.710321e-3,clcl4=-0.270742e-5
3296! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3297 ftlclo=t-tdpd*(clcl1+clcl2*t+tdpd*(clcl3+clcl4*t))
3298! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3299 end function
3300!-------------------------------------------------------------------------------
3322!\param[out] ftlclx real, temperature at the LCL in Kelvin
3323 elemental function ftlclx(t,tdpd)
3324!$$$ Subprogram documentation block
3325!
3326! Subprogram: ftlclx Compute LCL temperature.
3327! Author: Iredell org: w/NMC2X2 Date: 25 March 1999
3328!
3329! Abstract: Compute temperature at the lifting condensation level
3330! from temperature and dewpoint depression. A parcel lifted
3331! adiabatically becomes saturated at the lifting condensation level.
3332! The water model assumes a perfect gas, constant specific heats
3333! for gas and liquid, and neglects the volume of the liquid.
3334! The model does account for the variation of the latent heat
3335! of condensation with temperature. The ice option is not included.
3336! The Clausius-Clapeyron equation is integrated from the triple point
3337! to get the formulas
3338! pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl))
3339! pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew))
3340! where pvlcl is the saturated parcel vapor pressure at the LCL,
3341! pvdew is the unsaturated parcel vapor pressure initially,
3342! trlcl is ttp/tlcl and trdew is ttp/tdew. The adiabatic lifting
3343! of the parcel is represented by the following formula
3344! pvdew=pvlcl*(t/tlcl)**(1/kappa)
3345! This formula is inverted by iterating Newtonian approximations
3346! until tlcl is found to within 1.e-6 Kelvin. Note that the minimum
3347! returned temperature is 180 Kelvin.
3348!
3349! Program History Log:
3350! 1999-03-25 Iredell
3351!
3352! Usage: tlcl=ftlclx(t,tdpd)
3353!
3354! Input argument list:
3355! t Real(krealfp) temperature in Kelvin
3356! tdpd Real(krealfp) dewpoint depression in Kelvin
3357!
3358! Output argument list:
3359! ftlclx Real(krealfp) temperature at the LCL in Kelvin
3360!
3361! Attributes:
3362! Language: Fortran 90.
3363!
3364!$$$
3365 implicit none
3366 real(krealfp) ftlclx
3367 real(krealfp),intent(in):: t,tdpd
3368 real(krealfp),parameter:: terrm=1.e-4,tlmin=180.,tlminx=tlmin-5.
3369 real(krealfp) tr,pvdew,tlcl,ta,pvlcl,el,dpvlcl,terr,terrp
3370 integer i
3371! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3372 tr=con_ttp/(t-tdpd)
3373 pvdew=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))
3374 tlcl=t-tdpd
3375 do i=1,100
3376 tr=con_ttp/tlcl
3377 ta=t/tlcl
3378 pvlcl=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))*ta**(1/con_rocp)
3379 el=con_hvap+con_dldt*(tlcl-con_ttp)
3380 dpvlcl=(el/(con_rv*t**2)+1/(con_rocp*tlcl))*pvlcl
3381 terr=(pvlcl-pvdew)/dpvlcl
3382 tlcl=tlcl-terr
3383 if(abs(terr).le.terrm.or.tlcl.lt.tlminx) exit
3384 enddo
3385 ftlclx=max(tlcl,tlmin)
3386! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3387 end function
3388!-------------------------------------------------------------------------------
3393 subroutine gfuncphys
3394!$$$ Subprogram Documentation Block
3395!
3396! Subprogram: gfuncphys Compute all physics function tables
3397! Author: N Phillips w/NMC2X2 Date: 30 dec 82
3398!
3399! Abstract: Compute all physics function tables. Lookup tables are
3400! set up for computing saturation vapor pressure, dewpoint temperature,
3401! equivalent potential temperature, moist adiabatic temperature and humidity,
3402! pressure to the kappa, and lifting condensation level temperature.
3403!
3404! Program History Log:
3405! 1999-03-01 Iredell f90 module
3406!
3407! Usage: call gfuncphys
3408!
3409! Subprograms called:
3410! gpvsl compute saturation vapor pressure over liquid table
3411! gpvsi compute saturation vapor pressure over ice table
3412! gpvs compute saturation vapor pressure table
3413! gtdpl compute dewpoint temperature over liquid table
3414! gtdpi compute dewpoint temperature over ice table
3415! gtdp compute dewpoint temperature table
3416! gthe compute equivalent potential temperature table
3417! gtma compute moist adiabat tables
3418! gpkap compute pressure to the kappa table
3419! grkap compute pressure to the 1/kappa table
3420! gtlcl compute LCL temperature table
3421!
3422! Attributes:
3423! Language: Fortran 90.
3424!
3425!$$$
3426 implicit none
3427! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3428 call gpvsl
3429 call gpvsi
3430 call gpvs
3431 call gtdpl
3432 call gtdpi
3433 call gtdp
3434 call gthe
3435 call gtma
3436 call gpkap
3437 call grkap
3438 call gtlcl
3439! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3440 end subroutine
3441!-------------------------------------------------------------------------------
3442end module
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36