FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
sat_vapor_pres.F90
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @defgroup sat_vapor_pres_mod sat_vapor_pres_mod
20!> @ingroup sat_vapor_pres
21!> @brief Routines for computing the saturation vapor pressure (es),
22!! the specific humidity (qs) and vapor mixing ratio (mrs)
23!> Given a specified relative humidity, calculates es, qs, and mrs, as well as their
24!! derivatives with respect to temperature, and also includes routines
25!! to initialize the look-up table.
26!! This module contains routines for determining the saturation vapor
27!! pressure (<TT>ES</TT>) from lookup tables constructed using equations given
28!! in the Smithsonian tables. The <TT>ES</TT> lookup tables are valid between
29!! -160C and +100C (approx 113K to 373K).
30!!
31!! The values of <TT>ES</TT> are computed over ice from -160C to -20C,
32!! over water from 0C to 100C, and a blended value (over water and ice)
33!! from -20C to 0C.
34!!
35!! Routines are also included to calculate the saturation specific
36!! humidity and saturation mixing ratio for vapor, and their deriv-
37!! atives with respect to temperature. By default, the values returned
38!! are those at saturation; optionally, values of q and mr at a spec-
39!! ified relative humidity may instead be returned. Two forms are
40!! available; the approximate form that has been traditionally used in
41!! GCMs, and an exact form provided by SJ Lin in which saturation is
42!! reached while maintaining constant pressure and temperature.
43!!
44!! This version was written for non-vector machines.
45!! See the <LINK SRC="#NOTES">notes</LINK> section for details on vectorization.
46!!
47!! arguments
48!! ---------
49!! temp intent in temperature in degrees kelvin
50!! es intent out saturation vapor pressure in Pascals
51!! des intent out derivative of saturation vapor pressure
52!! with respect to temperature
53!! (Pascals/degree)
54!! press intent in atmospheric pressure in Pascals
55!! qs intent out specific humidity at relative humidity hc
56!! (kg(vapor) / kg(moist air)
57!! mrs intent out mixing ratio at relative humidity hc
58!! (kg(vapor) / kg(dry air)
59!!
60!! optional arguments
61!! ------------------
62!! q intent in vapor specific humidity
63!! (kg(vapor) / kg(moist air)
64!! hc intent in relative humidity at which output
65!! fields are desired: default is 100 %
66!! dqsdT intent out derivative of saturation specific
67!! humidity with respect to temperature
68!! (kg(vapor) / kg(moist air) /degree)
69!! mr intent in vapor mixing ratio
70!! (kg(vapor) / kg(dry air)
71!! dmrsdT intent out derivative of saturation mixing ratio
72!! with respect to temperature
73!! (kg(vapor) / kg(dry air) /degree)
74!! esat intent out saturation vapor pressure
75!! (Pascals)
76!! err_msg intent out character string to hold error message
77!! es_over_liq
78!! intent in use es table wrt liquid only
79!!
80!! Example Usages:
81!!
82!! call lookup_es (temp, es, err_msg)
83!!
84!! call lookup_des (temp, des, err_msg)
85!!
86!! call lookup_es_des (temp, es, des, err_msg)
87!!
88!! call lookup_es2 (temp, es, err_msg)
89!!
90!! call lookup_des2 (temp, des, err_msg)
91!!
92!! call lookup_es2_des2 (temp, es, des, err_msg)
93!!
94!! call compute_qs (temp, press, qs, q, hc, dqsdT, esat,
95!! err_msg, es_over_liq)
96!!
97!! call compute_mrs (temp, press, mrs, mr, hc, dmrsdT, esat,
98!! err_msg, es_over_liq)
99
100module sat_vapor_pres_mod
101
102!-----------------------------------------------------------------------
103!
104!
105! arguments
106! ---------
107! temp intent in temperature in degrees kelvin
108! es intent out saturation vapor pressure in Pascals
109! des intent out derivative of saturation vapor pressure
110! with respect to temperature
111! (Pascals/degree)
112! press intent in atmospheric pressure in Pascals
113! qs intent out specific humidity at relative humidity hc
114! (kg(vapor) / kg(moist air)
115! mrs intent out mixing ratio at relative humidity hc
116! (kg(vapor) / kg(dry air)
117!
118! optional arguments
119! ------------------
120! q intent in vapor specific humidity
121! (kg(vapor) / kg(moist air)
122! hc intent in relative humidity at which output
123! fields are desired: default is 100 %
124! dqsdT intent out derivative of saturation specific
125! humidity with respect to temperature
126! (kg(vapor) / kg(moist air) /degree)
127! mr intent in vapor mixing ratio
128! (kg(vapor) / kg(dry air)
129! dmrsdT intent out derivative of saturation mixing ratio
130! with respect to temperature
131! (kg(vapor) / kg(dry air) /degree)
132! esat intent out saturation vapor pressure
133! (Pascals)
134! err_msg intent out character string to hold error message
135! es_over_liq
136! intent in use es table wrt liquid only
137!
138!-----------------------------------------------------------------------
139
140! <CONTACT EMAIL="Bruce.Wyman@noaa.gov">
141! Bruce Wyman
142! </CONTACT>
143
144! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>
145
146! <OVERVIEW>
147! Routines for determining the saturation vapor pressure
148! (<TT>ES</TT>), saturation vapor specific humidity and saturation
149! vapor mixing ratio, and their derivatives with respect to
150! temperature.
151! </OVERVIEW>
152
153! <DESCRIPTION>
154! This module contains routines for determining the saturation vapor
155! pressure (<TT>ES</TT>) from lookup tables constructed using equations given
156! in the Smithsonian tables. The <TT>ES</TT> lookup tables are valid between
157! -160C and +100C (approx 113K to 373K).
158
159! The values of <TT>ES</TT> are computed over ice from -160C to -20C,
160! over water from 0C to 100C, and a blended value (over water and ice)
161! from -20C to 0C.
162
163! Routines are also included to calculate the saturation specific
164! humidity and saturation mixing ratio for vapor, and their deriv-
165! atives with respect to temperature. By default, the values returned
166! are those at saturation; optionally, values of q and mr at a spec-
167! ified relative humidity may instead be returned. Two forms are
168! available; the approximate form that has been traditionally used in
169! GCMs, and an exact form provided by SJ Lin in which saturation is
170! reached while maintaining constant pressure and temperature.
171
172! This version was written for non-vector machines.
173! See the <LINK SRC="#NOTES">notes</LINK> section for details on vectorization.
174
175! </DESCRIPTION>
176
177! <PUBLIC>
178! Description summarizing public interface.
179! </PUBLIC>
180
181 use constants_mod, only: tfreeze, rdgas, rvgas, hlv, es0
182 use fms_mod, only: write_version_number, stdout, stdlog, mpp_pe, mpp_root_pe, &
183 mpp_error, fatal, warning, fms_error_handler, &
184 error_mesg, check_nml_error
185 use mpp_mod, only: input_nml_file
186 use sat_vapor_pres_k_mod, only: sat_vapor_pres_init_k, lookup_es_k, &
188 lookup_es2_k, &
190 lookup_es3_k, &
193 use platform_mod, only: r4_kind, r8_kind
194
195implicit none
196private
197
198 public :: lookup_es, lookup_des, sat_vapor_pres_init
199 public :: lookup_es2, lookup_des2, lookup_es2_des2
200 public :: lookup_es3, lookup_des3, lookup_es3_des3
201 public :: lookup_es_des, compute_qs, compute_mrs
202!public :: compute_es
203 public :: escomp, descomp ! for backward compatibility
204 ! use lookup_es, lookup_des instead
205
206!-----------------------------------------------------------------------
207
208! <INTERFACE NAME="lookup_es">
209
210! <OVERVIEW>
211! For the given temperatures, returns the saturation vapor pressures.
212! </OVERVIEW>
213! <DESCRIPTION>
214! For the given temperatures these routines return the
215! saturation vapor pressure (esat). The return values are derived from
216! lookup tables (see notes below).
217! </DESCRIPTION>
218! <TEMPLATE>
219! call lookup_es( temp, esat, err_msg )
220! </TEMPLATE>
221! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
222! Temperature in degrees Kelvin.
223! </IN>
224! <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
225! Saturation vapor pressure in pascals.
226! May be a scalar, 1d, 2d, or 3d array.
227! Must have the same order and size as temp.
228! </OUT>
229! <OUT NAME="err_msg" UNITS=" " TYPE="character">
230! Character string containing error message to be returned to
231! calling routine.
232! </OUT>
233! <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
234! Temperature(s) provided to the saturation vapor pressure lookup
235! are outside the valid range of the lookup table (-160 to 100 deg C).
236! This may be due to a numerical instability in the model.
237! Information should have been printed to standard output to help
238! determine where the instability may have occurred.
239! If the lookup table needs a larger temperature range,
240! then parameters in the module header must be modified.
241! </ERROR> *
242
243 !> @brief For the given temperatures, returns the saturation vapor pressures
244 !!
245 !> For the given temperatures these routines return the saturation vapor pressure(esat).
246 !! The return values are derived from lookup tables.
247 !! Example usage:
248 !! @code{.F90} call lookup_es( temp, esat, err_msg ) @endcode
249 !!
250 !! @param temp Temperature in degrees Kelvin.
251 !! @param esat Saturation vapor pressure in pascals.
252 !! May be a scalar, 1d, 2d, or 3d array
253 !! Must have the same order and size as temp.
254 !! @param err_msg Character string containing error message to be returned to
255 !! calling routine.
256 !! @throws FATAL table overflow, nbad=##
257 !! Temperature(s) provided to the saturation vapor pressure lookup
258 !! are outside the valid range of the lookup table (-160 to 100 deg C).
259 !! This may be due to a numerical instability in the model.
260 !! Information should have been printed to standard output to help
261 !! determine where the instability may have occurred.
262 !! If the lookup table needs a larger temperature range,
263 !! then parameters in the module header must be modified.
264 !> @ingroup sat_vapor_pres_mod
265 interface lookup_es
266 module procedure lookup_es_0d_r4, lookup_es_0d_r8
267 module procedure lookup_es_1d_r4, lookup_es_1d_r8
268 module procedure lookup_es_2d_r4, lookup_es_2d_r8
269 module procedure lookup_es_3d_r4, lookup_es_3d_r8
270 end interface lookup_es
271 !> Provided for backward compatibility (to be removed soon)
272 !> @ingroup sat_vapor_pres_mod
273 interface escomp
274 module procedure lookup_es_0d_r4, lookup_es_0d_r8
275 module procedure lookup_es_1d_r4, lookup_es_1d_r8
276 module procedure lookup_es_2d_r4, lookup_es_2d_r8
277 module procedure lookup_es_3d_r4, lookup_es_3d_r8
278 end interface escomp
279! </INTERFACE>
280!-----------------------------------------------------------------------
281! <INTERFACE NAME="lookup_des">
282
283! <OVERVIEW>
284! For the given temperatures, returns the derivative of saturation vapor pressure
285! with respect to temperature.
286! </OVERVIEW>
287! <DESCRIPTION>
288! For the given temperatures these routines return the derivative of esat w.r.t.
289! temperature (desat). The return values are derived from
290! lookup tables (see notes below).
291! </DESCRIPTION>
292! <TEMPLATE>
293! call lookup_des( temp, desat )
294! </TEMPLATE>
295! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
296! Temperature in degrees Kelvin.
297! </IN>
298! <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
299! Derivative of saturation vapor pressure w.r.t. temperature
300! in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
301! Must have the same order and size as temp.
302! </OUT>
303! <OUT NAME="err_msg" UNITS=" " TYPE="character">
304! Character string containing error message to be returned to
305! calling routine.
306! </OUT>
307! <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
308! Temperature(s) provided to the saturation vapor pressure lookup
309! are outside the valid range of the lookup table (-160 to 100 deg C).
310! This may be due to a numerical instability in the model.
311! Information should have been printed to standard output to help
312! determine where the instability may have occurred.
313! If the lookup table needs a larger temperature range,
314! then parameters in the module header must be modified.
315! </ERROR> *
316
317 !> For the given temperatures, returns the derivative of saturation vapor pressure
318 !! with respect to temperature.
319 !!
320 !! For the given temperatures these routines return the derivtive of esat w.r.t. temperature
321 !! (desat). The return values are derived from lookup tables.
322 !!
323 !! @param [in] temp Temperature in degrees kelvin
324 !! @param [out] desat Derivative of saturation vapor pressure w.r.t. temperature
325 !! in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
326 !! Must have the same order and size as temp.
327 !! @param [out] err_msg Character string containing error message to be returned to
328 !! calling routine.
329 !!
330 !! @error FATAL table overflow, nbad=##
331 !! Temperature(s) provided to the saturation vapor pressure lookup
332 !! are outside the valid range of the lookup table (-160 to 100 deg C).
333 !! This may be due to a numerical instability in the model.
334 !! Information should have been printed to standard output to help
335 !! determine where the instability may have occurred.
336 !! If the lookup table needs a larger temperature range,
337 !! then parameters in the module header must be modified.
338 !!
339 !! <br>Example usage:
340 !! @code{.F90} call lookup_des( temp, desat) @endcode
341 !> @ingroup sat_vapor_pres_mod
342 interface lookup_des
343 module procedure lookup_des_0d_r4, lookup_des_0d_r8
344 module procedure lookup_des_1d_r4, lookup_des_1d_r8
345 module procedure lookup_des_2d_r4, lookup_des_2d_r8
346 module procedure lookup_des_3d_r4, lookup_des_3d_r8
347 end interface lookup_des
348! </INTERFACE>
349 !> Provided for backward compatibility (to be removed soon)
350 !> @ingroup sat_vapor_pres_mod
351 interface descomp
352 module procedure lookup_des_0d_r4, lookup_des_0d_r8
353 module procedure lookup_des_1d_r4, lookup_des_1d_r8
354 module procedure lookup_des_2d_r4, lookup_des_2d_r8
355 module procedure lookup_des_3d_r4, lookup_des_3d_r8
356 end interface descomp
357
358!-----------------------------------------------------------------------
359
360! <INTERFACE NAME="lookup_es_des">
361
362! <OVERVIEW>
363! For the given temperatures, returns the saturation vapor pressure
364! and the derivative of saturation vapor pressure with respect to
365! temperature.
366! </OVERVIEW>
367! <DESCRIPTION>
368! For the given temperatures these routines return the
369! saturation vapor pressure (esat) and the derivative of esat w.r.t
370! temperature (desat). The return values are derived from
371! lookup tables (see notes below).
372! </DESCRIPTION>
373! <TEMPLATE>
374! call lookup_es_des( temp, esat, desat, err_msg )
375! </TEMPLATE>
376! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
377! Temperature in degrees Kelvin.
378! </IN>
379! <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
380! Saturation vapor pressure in pascals.
381! May be a scalar, 1d, 2d, or 3d array.
382! Must have the same order and size as temp.
383! </OUT>
384! <OUT NAME="desat" UNITS="pascal/ degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
385! Derivative of saturation vapor pressure w.r.t. temperature
386! in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
387! Must have the same order and size as temp.
388! </OUT>
389! <OUT NAME="err_msg" UNITS=" " TYPE="character">
390! Character string containing error message to be returned to
391! calling routine.
392! </OUT>
393! <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
394! Temperature(s) provided to the saturation vapor pressure lookup
395! are outside the valid range of the lookup table (-160 to 100 deg C).
396! This may be due to a numerical instability in the model.
397! Information should have been printed to standard output to help
398! determine where the instability may have occurred.
399! If the lookup table needs a larger temperature range,
400! then parameters in the module header must be modified.
401! </ERROR> *
402
403 !> @brief For the given temperatures, returns the saturation vapor pressure
404 !! and the derivative of saturation vapor pressure with respect to
405 !! temperature.
406 !!
407 !> For the given temperatures these routines return the
408 !! saturation vapor pressure (esat) and the derivative of esat w.r.t
409 !! temperature (desat). The return values are derived from
410 !! lookup tables (see notes below).
411 !!
412 !! <br>Example usage:
413 !! @code{.F90} call lookup_es_des( temp, esat, desat, err_msg ) @endcode
414 !!
415 !! @param temp Temperature in degrees Kelvin.
416 !! @param [out] esat Saturation vapor pressure in pascals. May be a scalar, 1d, 2d, or 3d array.
417 !! Must have the same order and size as temp.
418 !! @param [out] desat Derivative of saturation vapor pressure w.r.t. temperature
419 !! in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
420 !! Must have the same order and size as temp.
421 !! @param [out] err_msg Character string containing error message to be returned to
422 !! calling routine.
423 !! @error FATAL table overflow, nbad=##
424 !! Temperature(s) provided to the saturation vapor pressure lookup
425 !! are outside the valid range of the lookup table (-160 to 100 deg C).
426 !! This may be due to a numerical instability in the model.
427 !! Information should have been printed to standard output to help
428 !! determine where the instability may have occurred.
429 !! If the lookup table needs a larger temperature range,
430 !! then parameters in the module header must be modified.
431 !> @ingroup sat_vapor_pres_mod
432 interface lookup_es_des
433 module procedure lookup_es_des_0d_r4, lookup_es_des_0d_r8
434 module procedure lookup_es_des_1d_r4, lookup_es_des_1d_r8
435 module procedure lookup_es_des_2d_r4, lookup_es_des_2d_r8
436 module procedure lookup_es_des_3d_r4, lookup_es_des_3d_r8
437 end interface lookup_es_des
438
439 !> @ingroup sat_vapor_pres_mod
440 interface lookup_es2
441 module procedure lookup_es2_0d_r4, lookup_es2_0d_r8
442 module procedure lookup_es2_1d_r4, lookup_es2_1d_r8
443 module procedure lookup_es2_2d_r4, lookup_es2_2d_r8
444 module procedure lookup_es2_3d_r4, lookup_es2_3d_r8
445 end interface lookup_es2
446
447 !> @ingroup sat_vapor_pres_mod
448 interface lookup_des2
449 module procedure lookup_des2_0d_r4, lookup_des2_0d_r8
450 module procedure lookup_des2_1d_r4, lookup_des2_1d_r8
451 module procedure lookup_des2_2d_r4, lookup_des2_2d_r8
452 module procedure lookup_des2_3d_r4, lookup_des2_3d_r8
453 end interface lookup_des2
454
455 !> @ingroup sat_vapor_pres_mod
456 interface lookup_es2_des2
457 module procedure lookup_es2_des2_0d_r4, lookup_es2_des2_0d_r8
458 module procedure lookup_es2_des2_1d_r4, lookup_es2_des2_1d_r8
459 module procedure lookup_es2_des2_2d_r4, lookup_es2_des2_2d_r8
460 module procedure lookup_es2_des2_3d_r4, lookup_es2_des2_3d_r8
461 end interface lookup_es2_des2
462
463 !> @ingroup sat_vapor_pres_mod
464 interface lookup_es3
465 module procedure lookup_es3_0d_r4, lookup_es3_0d_r8
466 module procedure lookup_es3_1d_r4, lookup_es3_1d_r8
467 module procedure lookup_es3_2d_r4, lookup_es3_2d_r8
468 module procedure lookup_es3_3d_r4, lookup_es3_3d_r8
469 end interface lookup_es3
470
471 !> @ingroup sat_vapor_pres_mod
472 interface lookup_des3
473 module procedure lookup_des3_0d_r4, lookup_des3_0d_r8
474 module procedure lookup_des3_1d_r4, lookup_des3_1d_r8
475 module procedure lookup_des3_2d_r4, lookup_des3_2d_r8
476 module procedure lookup_des3_3d_r4, lookup_des3_3d_r8
477 end interface lookup_des3
478
479 !> @ingroup sat_vapor_pres_mod
480 interface lookup_es3_des3
481 module procedure lookup_es3_des3_0d_r4, lookup_es3_des3_0d_r8
482 module procedure lookup_es3_des3_1d_r4, lookup_es3_des3_1d_r8
483 module procedure lookup_es3_des3_2d_r4, lookup_es3_des3_2d_r8
484 module procedure lookup_es3_des3_3d_r4, lookup_es3_des3_3d_r8
485 end interface lookup_es3_des3
486
487!-----------------------------------------------------------------------
488
489! <INTERFACE NAME="compute_qs">
490
491! <OVERVIEW>
492! For the given temperatures, pressures and optionally vapor
493! specific humidity, returns the specific humidity at saturation
494! (optionally at relative humidity hc instead of at saturation) and
495! optionally the derivative of saturation specific humidity w.r.t.
496! temperature, and the saturation vapor pressure.
497! </OVERVIEW>
498! <DESCRIPTION>
499! For the input temperature and pressure these routines return the
500! specific humidity (qsat) at saturation (unless optional argument
501! hc is used to specify the relative humidity at which qsat should
502! apply) and, if desired, the derivative of qsat w.r.t temperature
503! (dqsdT) and / or the saturation vapor pressure (esat). If the
504! optional input argument specific humidity (q) is present, the
505! exact expression for qs is used; if q is not present the tradit-
506! ional form (valid at saturation) is used. if the optional qsat
507! derivative argument is present, the derivative of qsat w.r.t.
508! temperature will also be returned, defined consistent with the
509! expression used for qsat. The return values are derived from
510! lookup tables (see notes below).
511! </DESCRIPTION>
512! <TEMPLATE>
513! call compute_qs( temp, press, qsat, q, hc, dqsdT, esat, err_msg )
514! </TEMPLATE>
515! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
516! Temperature in degrees Kelvin.
517! </IN>
518! <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
519! Air pressure in Pascals.
520! </IN>
521! <OUT NAME="qsat" UNITS="kg(vapor) / kg(moist air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
522! Specific humidity in kg (vapor) / kg (moist air)
523! May be a scalar, 1d, 2d, or 3d array.
524! Must have the same order and size as temp.
525! </OUT>
526! <IN NAME="q" UNIT="kg(vapor) / kg (moist air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
527! Vapor specific humidity in kg (vapor) / kg (moist air).
528! If present, exact formulation for qsat and dqsdT will be used.
529! </IN>
530! <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)">
531! Relative humidity at which output variables are desired.
532! If not present, values will apply at saturation.
533! </IN>
534! <OUT NAME="dqsdT" UNITS="kg(vapor) / kg(moist air) / degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
535! Derivative of saturation specific humidity w.r.t. temperature
536! in kg(vapor) / kg(moist air) / degree. May be a
537! scalar, 1d, 2d, or 3d array.
538! Must have the same order and size as temp.
539! </OUT>
540! <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
541! Saturation vapor pressure. May be a scalar, 1d, 2d, or 3d array.
542! Must have the same order and size as temp.
543! </OUT>
544! <OUT NAME="err_msg" UNITS=" " TYPE="character">
545! Character string containing error message to be returned to
546! calling routine.
547! </OUT>
548! <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
549! Temperature(s) provided to the saturation vapor pressure lookup
550! are outside the valid range of the lookup table (-160 to 100 deg C).
551! This may be due to a numerical instability in the model.
552! Information should have been printed to standard output to help
553! determine where the instability may have occurred.
554! If the lookup table needs a larger temperature range,
555! then parameters in the module header must be modified.
556! </ERROR> *
557
558 !> @brief For the given temperatures, pressures and optionally vapor
559 !! specific humidity, returns the specific humidity at saturation
560 !! (optionally at relative humidity hc instead of at saturation) and
561 !! optionally the derivative of saturation specific humidity w.r.t.
562 !! temperature, and the saturation vapor pressure.
563 !!
564 !! For the input temperature and pressure these routines return the
565 !! specific humidity (qsat) at saturation (unless optional argument
566 !! hc is used to specify the relative humidity at which qsat should
567 !! apply) and, if desired, the derivative of qsat w.r.t temperature
568 !! (dqsdT) and / or the saturation vapor pressure (esat). If the
569 !! optional input argument specific humidity (q) is present, the
570 !! exact expression for qs is used; if q is not present the tradit-
571 !! ional form (valid at saturation) is used. if the optional qsat
572 !! derivative argument is present, the derivative of qsat w.r.t.
573 !! temperature will also be returned, defined consistent with the
574 !! expression used for qsat. The return values are derived from
575 !! lookup tables (see notes below).
576 !!
577 !! Example usage:
578 !! @code{.F90} call compute_qs( temp, press, qsat, q, hc, dqsdT, esat, err_msg ) @endcode
579 !!
580 !> @ingroup sat_vapor_pres_mod
581 interface compute_qs
582 module procedure compute_qs_0d_r4, compute_qs_0d_r8
583 module procedure compute_qs_1d_r4, compute_qs_1d_r8
584 module procedure compute_qs_2d_r4, compute_qs_2d_r8
585 module procedure compute_qs_3d_r4, compute_qs_3d_r8
586 end interface compute_qs
587
588!-----------------------------------------------------------------------
589
590! <INTERFACE NAME="compute_mrs">
591
592! <OVERVIEW>
593! For the given temperatures, pressures and optionally vapor
594! mixing ratio, returns the vapor mixing ratio at saturation
595! (optionally at relative humidity hc instead of at saturation) and
596! optionally the derivative of saturation vapor mixing ratio w.r.t.
597! temperature, and the saturation vapor pressure.
598! </OVERVIEW>
599! <DESCRIPTION>
600! For the input temperature and pressure these routines return the
601! vapor mixing ratio (mrsat) at saturation (unless optional argument
602! hc is used to specify the relative humidity at which mrsat should
603! apply) and, if desired, the derivative of mrsat w.r.t temperature
604! (dmrsdT) and / or the saturation vapor pressure (esat). If the
605! optional input argument specific humidity (mr) is present, the
606! exact expression for mrs is used; if qr is not present the tradit-
607! ional form (valid at saturation) is used. if the optional mrsat
608! derivative argument is present, the derivative of mrsat w.r.t.
609! temperature will also be returned, defined consistent with the
610! expression used for mrsat. The return values are derived from
611! lookup tables (see notes below).
612! </DESCRIPTION>
613! <TEMPLATE>
614! call compute_mrs( temp, press, mrsat, mr, hc, dmrsdT, esat,
615! err_msg )
616! </TEMPLATE>
617! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
618! Temperature in degrees Kelvin.
619! </IN>
620! <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
621! Air pressure in Pascals.
622! </IN>
623! <OUT NAME="mrsat" UNITS="kg(vapor) / kg (dry air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
624! Vapor mixing ratio in kg (vapor) / kg (dry air)
625! May be a scalar, 1d, 2d, or 3d array.
626! Must have the same order and size as temp.
627! </OUT>
628! <IN NAME="mr" UNIT="kg(vapor) / kg (dry air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
629! Vapor mixing ratio in kg (vapor) / kg (dry air).
630! If present, exact formulation for mrsat and dmrsdT will be used.
631! </IN>
632! <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)">
633! Relative humidity at which output variables are desired.
634! If not present, values will apply at saturation.
635! </IN>
636! <OUT NAME="dmrsdT" UNITS="kg(vapor) / kg(dry air) / degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
637! Derivative of saturation vapor mixing ratio w.r.t. temperature
638! in kg(vapor) / kg(dry air) / degree. May be a
639! scalar, 1d, 2d, or 3d array.
640! Must have the same order and size as temp.
641! </OUT>
642! <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
643! Saturation vapor pressure. May be a scalar, 1d, 2d, or 3d array.
644! Must have the same order and size as temp.
645! </OUT>
646! <OUT NAME="err_msg" UNITS=" " TYPE="character">
647! Character string containing error message to be returned to
648! calling routine.
649! </OUT>
650! <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
651! Temperature(s) provided to the saturation vapor pressure lookup
652! are outside the valid range of the lookup table (-160 to 100 deg C).
653! This may be due to a numerical instability in the model.
654! Information should have been printed to standard output to help
655! determine where the instability may have occurred.
656! If the lookup table needs a larger temperature range,
657! then parameters in the module header must be modified.
658! </ERROR> *
659
660 !> For the given temperatures, pressures and optionally vapor
661 !! mixing ratio, returns the vapor mixing ratio at saturation
662 !! (optionally at relative humidity hc instead of at saturation) and
663 !! optionally the derivative of saturation vapor mixing ratio w.r.t.
664 !! temperature, and the saturation vapor pressure.
665 !!
666 !! For the input temperature and pressure these routines return the
667 !! vapor mixing ratio (mrsat) at saturation (unless optional argument
668 !! hc is used to specify the relative humidity at which mrsat should
669 !! apply) and, if desired, the derivative of mrsat w.r.t temperature
670 !! (dmrsdT) and / or the saturation vapor pressure (esat). If the
671 !! optional input argument specific humidity (mr) is present, the
672 !! exact expression for mrs is used; if qr is not present the tradit-
673 !! ional form (valid at saturation) is used. if the optional mrsat
674 !! derivative argument is present, the derivative of mrsat w.r.t.
675 !! temperature will also be returned, defined consistent with the
676 !! expression used for mrsat. The return values are derived from
677 !! lookup tables (see notes below).
678 !!
679 !! <br>Example usage:
680 !! @code{.F90} call compute_mrs( temp, press, mrsat, mr, hc, dmrsdT, esat,
681 !! err_msg ) @endcode
682 !> @ingroup sat_vapor_pres_mod
683 interface compute_mrs
684 module procedure compute_mrs_0d_r4, compute_mrs_0d_r8
685 module procedure compute_mrs_1d_r4, compute_mrs_1d_r8
686 module procedure compute_mrs_2d_r4, compute_mrs_2d_r8
687 module procedure compute_mrs_3d_r4, compute_mrs_3d_r8
688 end interface compute_mrs
689
690!-----------------------------------------------------------------------
691! <INTERFACE NAME="compute_es">
692
693! <OVERVIEW>
694! For the given temperatures, computes the saturation vapor pressures.
695! </OVERVIEW>
696! <DESCRIPTION>
697! Computes saturation vapor pressure for the given temperature using
698! the equations given in the Smithsonian Meteorological Tables.
699! Between -20C and 0C a blended value over ice and water is returned.
700! </DESCRIPTION>
701! <TEMPLATE>
702! es = compute_es ( temp )
703! </TEMPLATE>
704! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
705! Temperature in degrees Kelvin.
706! </IN>
707! <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
708! Saturation vapor pressure in pascals.
709! May be a scalar, 1d, 2d, or 3d array.
710! Must have the same order and size as temp.
711! </OUT>
712
713!interface compute_es
714! module procedure compute_es_0d, compute_es_1d, compute_es_2d, compute_es_3d
715!end interface
716! </INTERFACE>
717!-----------------------------------------------------------------------
718 !> @ingroup sat_vapor_pres_mod
719 interface check_1d
720 module procedure check_1d_r4, check_1d_r8
721 end interface check_1d
722
723 interface check_2d
724 module procedure check_2d_r4, check_2d_r8
725 end interface check_2d
726
727 !> @ingroup sat_vapor_pres_mod
728 interface temp_check
729 module procedure temp_check_1d_r4, temp_check_1d_r8
730 module procedure temp_check_2d_r4, temp_check_2d_r8
731 module procedure temp_check_3d_r4, temp_check_3d_r8
732 end interface temp_check
733
734 !> @ingroup sat_vapor_pres_mod
735 interface show_all_bad
736 module procedure show_all_bad_0d_r4, show_all_bad_0d_r8
737 module procedure show_all_bad_1d_r4, show_all_bad_1d_r8
738 module procedure show_all_bad_2d_r4, show_all_bad_2d_r8
739 module procedure show_all_bad_3d_r4, show_all_bad_3d_r8
740 end interface show_all_bad
741
742!> @addtogroup sat_vapor_pres_mod
743!> @{
744!-----------------------------------------------------------------------
745! Include variable "version" to be written to log file.
746#include<file_version.h>
747
748 logical :: module_is_initialized = .false.
749
750!-----------------------------------------------------------------------
751! parameters for use in computing qs and mrs
752
753 real(r8_kind), parameter :: EPSILO = real(rdgas,r8_kind)/real(rvgas, r8_kind)
754 real(r8_kind), parameter :: ZVIR = real(rvgas,r8_kind)/real(rdgas,r8_kind) - 1.0_r8_kind
755
756!-----------------------------------------------------------------------
757! parameters for table size and resolution
758
759 integer, public :: tcmin = -160 ! minimum temperature (degC) in lookup table
760 integer, public :: tcmax = 100 ! maximum temperature (degC) in lookup table
761 integer :: esres = 10 ! table resolution (increments per degree)
762 integer :: nsize ! (tcmax-tcmin)*esres+1 ! lookup table size
763 integer :: nlim ! nsize-1
764
765 integer :: stdoutunit=0
766!-----------------------------------------------------------------------
767! variables needed by temp_check
768 real(r8_kind) :: tmin, dtinv, teps
769
770! The default values below preserve the behavior of omsk and earlier revisions.
771 logical :: show_bad_value_count_by_slice=.true.
772 logical :: show_all_bad_values=.false.
773 logical :: use_exact_qs = .false.
774 logical :: do_simple =.false.
775 logical :: construct_table_wrt_liq = .false.
776 logical :: construct_table_wrt_liq_and_ice = .false.
777
778 namelist / sat_vapor_pres_nml / show_bad_value_count_by_slice, show_all_bad_values, &
779 use_exact_qs, do_simple, &
780 construct_table_wrt_liq, &
781 construct_table_wrt_liq_and_ice
782
783contains
784
785 subroutine sat_vapor_pres_init(err_msg)
786
787! =================================================================
788! + +
789! + construction of the es table +
790! + +
791! + this table is constructed from es equations from the +
792! + smithsonian tables. the es input is computed from values +
793! + (in one-tenth of a degree increments) of es over ice +
794! + from -153c to 0c and values of es over water from 0c to 102c. +
795! + output table contains these data interleaved with their +
796! + derivatives with respect to temperature except between -20c +
797! + and 0c where blended (over water and over ice) es values and +
798! + derivatives are calculated. +
799! + note: all es computation is done in pascals +
800! =================================================================
801
802 character(len=*), intent(out), optional :: err_msg
803 character(len=128) :: err_msg_local
804 integer :: iunit, ierr, io
805
806! return silently if this routine has already been called
807 if (module_is_initialized) return
808
809!---- read namelist input ----
810 read (input_nml_file, sat_vapor_pres_nml, iostat=io)
811 ierr = check_nml_error(io,'sat_vapor_pres_nml')
812
813! write version number and namelist to log file
814 call write_version_number("SAT_VAPOR_PRES_MOD", version)
815 iunit = stdlog()
816 stdoutunit = stdout()
817 if (mpp_pe() == mpp_root_pe()) write (iunit, nml=sat_vapor_pres_nml)
818
819 if(do_simple) then
820 tcmin = -173
821 tcmax = 350
822 endif
823 nsize = (tcmax-tcmin)*esres+1
824 nlim = nsize-1
825 call sat_vapor_pres_init_k(nsize, real(tcmin,r8_kind), real(tcmax,r8_kind), &
826 real(TFREEZE,r8_kind), real(HLV,r8_kind),&
827 real(RVGAS,r8_kind), real(ES0,r8_kind), err_msg_local, use_exact_qs, do_simple,&
828 construct_table_wrt_liq, &
829 construct_table_wrt_liq_and_ice, &
830 teps, tmin, dtinv)
831 if ( err_msg_local == '' ) then
832 if(present(err_msg)) err_msg = ''
833 else
834 if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
835 endif
836
837 module_is_initialized = .true.
838
839end subroutine sat_vapor_pres_init
840
841#include "sat_vapor_pres_r4.fh"
842#include "sat_vapor_pres_r8.fh"
843
844!#######################################################################
845!#######################################################################
846!-------------------------------------------------------------------
847! Computation of the es values
848!
849! Saturation vapor pressure (es) values are computed from
850! equations in the Smithsonian meteorological tables page 350.
851! For temperatures < 0C, sat vapor pres is computed over ice.
852! For temperatures > -20C, sat vapor pres is computed over water.
853! Between -20C and 0C the returned value is blended (over water
854! and over ice). All sat vapor pres values are returned in pascals.
855!
856! Reference: Smithsonian meteorological tables, page 350.
857!-------------------------------------------------------------------
858
859! <FUNCTION NAME="compute_es_1d" INTERFACE="compute_es">
860! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
861! <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
862! </FUNCTION>
863!function compute_es_1d (tem) result (es)
864!real, intent(in) :: tem(:)
865!real :: es(size(tem,1))
866
867!es = compute_es_k(tem, TFREEZE)
868
869!end function compute_es_1d
870!--------------------------------------------------------
871
872! <FUNCTION NAME="compute_es_0d" INTERFACE="compute_es">
873! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
874! <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
875! </FUNCTION>
876!function compute_es_0d (tem) result (es)
877!real, intent(in) :: tem
878!real :: es
879!real, dimension(1) :: tem1, es1
880
881! tem1(1) = tem
882! es1 = compute_es_1d (tem1)
883! es = es1(1)
884
885!end function compute_es_0d
886
887!--------------------------
888
889! <FUNCTION NAME="compute_es_2d" INTERFACE="compute_es">
890! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
891! <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
892! </FUNCTION>
893!function compute_es_2d (tem) result (es)
894!real, intent(in) :: tem(:,:)
895!real, dimension(size(tem,1),size(tem,2)) :: es
896!integer :: j
897
898! do j = 1, size(tem,2)
899! es(:,j) = compute_es_1d (tem(:,j))
900! enddo
901
902!end function compute_es_2d
903
904!--------------------------
905! <FUNCTION NAME="compute_es_3d" INTERFACE="compute_es">
906! <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
907! <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
908! </FUNCTION>
909!function compute_es_3d (tem) result (es)
910!real, intent(in) :: tem(:,:,:)
911!real, dimension(size(tem,1),size(tem,2),size(tem,3)) :: es
912!integer :: j, k
913
914! do k = 1, size(tem,3)
915! do j = 1, size(tem,2)
916! es(:,j,k) = compute_es_1d (tem(:,j,k))
917! enddo
918! enddo
919
920!end function compute_es_3d
921
922!#######################################################################
923end module sat_vapor_pres_mod
924!#######################################################################
925
926! <INFO>
927
928! <REFERENCE>
929! Smithsonian Meteorological Tables Page 350.
930! </REFERENCE>
931
932! <BUG>
933! No error checking is done to make sure that the size of the
934! input and output fields match.
935! </BUG>
936
937! <NOTE>
938! 1. <B>Vectorization</B><BR/>
939! To create a vector version the lookup routines need to be modified.
940! The local variables: tmp, del, ind, should be changed to arrays
941! with the same size and order as input array temp.
942!
943! 2. <B>Construction of the <TT>ES</TT> tables</B><BR/>
944! The tables are constructed using the saturation vapor pressure (<TT>ES</TT>)
945! equations in the Smithsonian tables. The tables are valid between
946! -160C to +100C with increments at 1/10 degree. Between -160C and -20C
947! values of <TT>ES</TT> over ice are used, between 0C and 100C values of<TT> ES</TT>
948! over water are used, between -20C and 0C blended values of <TT>ES</TT>
949! (over water and over ice) are used.
950!
951! There are three tables constructed: <TT>ES</TT>, first derivative
952! (<TT>ES'</TT>), and
953! second derivative (<TT>ES</TT>''). The ES table is constructed directly from
954! the equations in the Smithsonian tables. The <TT>ES</TT>' table is constructed
955! by bracketing temperature values at +/- 0.01 degrees. The <TT>ES</TT>'' table
956! is estimated by using centered differencing of the <TT>ES</TT>' table.
957!
958! 3. <B>Determination of <TT>es</TT> and <TT>es'</TT> from lookup tables</B><BR/>
959! Values of the saturation vapor pressure (<TT>es</TT>) and the
960! derivative (<TT>es'</TT>) are determined at temperature (T) from the lookup
961! tables (<TT>ES</TT>, <TT>ES'</TT>, <TT>ES''</TT>)
962! using the following formula.
963!<PRE>
964! es (T) = ES(t) + ES'(t) * dt + 0.5 * ES''(t) * dt**2
965! es'(T) = ES'(t) + ES''(t) * dt
966!
967! where t = lookup table temperature closest to T
968! dt = T - t
969!</PRE>
970!
971! 4. Internal (private) parameters<BR/>
972! These parameters can be modified to increase/decrease the size/range
973! of the lookup tables.
974!<PRE>
975!! tcmin The minimum temperature (in deg C) in the lookup tables.
976!! [integer, default: tcmin = -160]
977!!
978!! tcmax The maximum temperature (in deg C) in the lookup tables.
979!! [integer, default: tcmin = +100]
980!!</PRE>
981!! </NOTE>
982!
983!! <TESTPROGRAM NAME="test_sat_vapor_pres">
984!<PRE>
985!use sat_vapor_pres_mod
986!implicit none
987!
988!integer, parameter :: ipts=500, jpts=100, kpts=50, nloop=1
989!real, dimension(ipts,jpts,kpts) :: t,es,esn,des,desn
990!integer :: n
991!
992!! generate temperatures between 120K and 340K
993! call random_number (t)
994! t = 130. + t * 200.
995!
996!! initialize the tables (optional)
997! call sat_vapor_pres_init
998!
999!! compute actual es and "almost" actual des
1000! es = compute_es (t)
1001! des = compute_des (t)
1002!
1003!do n = 1, nloop
1004!! es and des
1005! call lookup_es (t, esn)
1006! call lookup_des (t,desn)
1007!enddo
1008!
1009!! terminate, print deviation from actual
1010! print *, 'size=',ipts,jpts,kpts,nloop
1011! print *, 'err es = ', sum((esn-es)**2)
1012! print *, 'err des = ', sum((desn-des)**2)
1013!
1014!contains
1015!
1016!!----------------------------------
1017!! routine to estimate derivative
1018!
1019! function compute_des (tem) result (des)
1020! real, intent(in) :: tem(:,:,:)
1021! real, dimension(size(tem,1),size(tem,2),size(tem,3)) :: des,esp,esm
1022! real, parameter :: tdel = .01
1023! esp = compute_es (tem+tdel)
1024! esm = compute_es (tem-tdel)
1025! des = (esp-esm)/(2*tdel)
1026! end function compute_des
1027!!----------------------------------
1028!
1029!end program test_sat_vapor_pres
1030!</PRE>
1031! </TESTPROGRAM>
1032! </INFO>
1033!> @}
1034! close documentation grouping