sat_vapor_pres/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
 100 module 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, 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, &
 187                                    lookup_des_k, lookup_es_des_k, &
 188                                    lookup_es2_k,  &
 189                                    lookup_des2_k, lookup_es2_des2_k, &
 190                                    lookup_es3_k,  &
 191                                    lookup_des3_k, lookup_es3_des3_k, &
 192                                    compute_qs_k, compute_mrs_k
 193
 194  use platform_mod, only: r4_kind, r8_kind
 195
 196 implicit none
 197 private
 198
 199  public :: lookup_es, lookup_des, sat_vapor_pres_init
 200  public :: lookup_es2, lookup_des2, lookup_es2_des2
 201  public :: lookup_es3, lookup_des3, lookup_es3_des3
 202  public :: lookup_es_des, compute_qs, compute_mrs
 203 !public :: compute_es
 204  public :: escomp, descomp ! for backward compatibility
 205                            ! use lookup_es, lookup_des instead
 206
 207 !-----------------------------------------------------------------------
 208
 209 ! <INTERFACE NAME="lookup_es">
 210
 211 !   <OVERVIEW>
 212 !     For the given temperatures, returns the saturation vapor pressures.
 213 !   </OVERVIEW>
 214 !   <DESCRIPTION>
 215 !     For the given temperatures these routines return the
 216 !     saturation vapor pressure (esat). The return values are derived from
 217 !     lookup tables (see notes below).
 218 !   </DESCRIPTION>
 219 !   <TEMPLATE>
 220 !     call lookup_es( temp, esat, err_msg )
 221 !   </TEMPLATE>
 222 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 223 !     Temperature in degrees Kelvin.
 224 !   </IN>
 225 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 226 !     Saturation vapor pressure in pascals.
 227 !             May be a scalar, 1d, 2d, or 3d array.
 228 !             Must have the same order and size as temp.
 229 !   </OUT>
 230 !   <OUT NAME="err_msg" UNITS="      " TYPE="character">
 231 !     Character string containing error message to be returned to
 232 !     calling routine.
 233 !   </OUT>
 234 !   <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
 235 !     Temperature(s) provided to the saturation vapor pressure lookup
 236 !          are outside the valid range of the lookup table (-160 to 100 deg C).
 237 !          This may be due to a numerical instability in the model.
 238 !          Information should have been printed to standard output to help
 239 !          determine where the instability may have occurred.
 240 !          If the lookup table needs a larger temperature range,
 241 !          then parameters in the module header must be modified.
 242 !   </ERROR> *
 243
 244  !> @brief For the given temperatures, returns the saturation vapor pressures
 245  !!
 246  !> For the given temperatures these routines return the saturation vapor pressure(esat).
 247  !! The return values are derived from lookup tables.
 248  !! Example usage:
 249  !! @code{.F90} call lookup_es( temp, esat, err_msg ) @endcode
 250  !!
 251  !! @param temp Temperature in degrees Kelvin.
 252  !! @param esat Saturation vapor pressure in pascals.
 253  !!             May be a scalar, 1d, 2d, or 3d array
 254  !!             Must have the same order and size as temp.
 255  !! @param err_msg Character string containing error message to be returned to
 256  !!                calling routine.
 257  !! @throws FATAL table overflow, nbad=##
 258  !!     Temperature(s) provided to the saturation vapor pressure lookup
 259  !!          are outside the valid range of the lookup table (-160 to 100 deg C).
 260  !!          This may be due to a numerical instability in the model.
 261  !!          Information should have been printed to standard output to help
 262  !!          determine where the instability may have occurred.
 263  !!          If the lookup table needs a larger temperature range,
 264  !!          then parameters in the module header must be modified.
 265  !> @ingroup sat_vapor_pres_mod
 266  interface lookup_es
 267    module procedure lookup_es_0d, lookup_es_1d, lookup_es_2d, lookup_es_3d
 268  end interface
 269  !> Provided for backward compatibility (to be removed soon)
 270  !> @ingroup sat_vapor_pres_mod
 271  interface escomp
 272    module procedure lookup_es_0d, lookup_es_1d, lookup_es_2d, lookup_es_3d
 273  end interface
 274 ! </INTERFACE>
 275 !-----------------------------------------------------------------------
 276 ! <INTERFACE NAME="lookup_des">
 277
 278 !   <OVERVIEW>
 279 !     For the given temperatures, returns the derivative of saturation vapor pressure
 280 !     with respect to temperature.
 281 !   </OVERVIEW>
 282 !   <DESCRIPTION>
 283 !     For the given temperatures these routines return the derivative of esat w.r.t.
 284 !     temperature (desat). The return values are derived from
 285 !     lookup tables (see notes below).
 286 !   </DESCRIPTION>
 287 !   <TEMPLATE>
 288 !     call lookup_des( temp, desat )
 289 !   </TEMPLATE>
 290 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 291 !     Temperature in degrees Kelvin.
 292 !   </IN>
 293 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 294 !     Derivative of saturation vapor pressure w.r.t. temperature
 295 !                 in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
 296 !                 Must have the same order and size as temp.
 297 !   </OUT>
 298 !   <OUT NAME="err_msg" UNITS="      " TYPE="character">
 299 !     Character string containing error message to be returned to
 300 !     calling routine.
 301 !   </OUT>
 302 !   <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
 303 !     Temperature(s) provided to the saturation vapor pressure lookup
 304 !          are outside the valid range of the lookup table (-160 to 100 deg C).
 305 !          This may be due to a numerical instability in the model.
 306 !          Information should have been printed to standard output to help
 307 !          determine where the instability may have occurred.
 308 !          If the lookup table needs a larger temperature range,
 309 !          then parameters in the module header must be modified.
 310 !   </ERROR> *
 311
 312  !> For the given temperatures, returns the derivative of saturation vapor pressure
 313  !! with respect to temperature.
 314  !!
 315  !! For the given temperatures these routines return the derivtive of esat w.r.t. temperature
 316  !! (desat). The return values are derived from lookup tables.
 317  !!
 318  !! @param [in] temp Temperature in degrees kelvin
 319  !! @param [out] desat Derivative of saturation vapor pressure w.r.t. temperature
 320  !!                 in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
 321  !!                 Must have the same order and size as temp.
 322  !! @param [out] err_msg Character string containing error message to be returned to
 323  !!     calling routine.
 324  !!
 325  !! @error FATAL table overflow, nbad=##
 326  !!        Temperature(s) provided to the saturation vapor pressure lookup
 327  !!        are outside the valid range of the lookup table (-160 to 100 deg C).
 328  !!        This may be due to a numerical instability in the model.
 329  !!        Information should have been printed to standard output to help
 330  !!        determine where the instability may have occurred.
 331  !!        If the lookup table needs a larger temperature range,
 332  !!        then parameters in the module header must be modified.
 333  !!
 334  !! <br>Example usage:
 335  !! @code{.F90} call lookup_des( temp, desat) @endcode
 336  !> @ingroup sat_vapor_pres_mod
 337  interface lookup_des
 338    module procedure lookup_des_0d, lookup_des_1d, lookup_des_2d, lookup_des_3d
 339  end interface
 340 ! </INTERFACE>
 341  !> Provided for backward compatibility (to be removed soon)
 342  !> @ingroup sat_vapor_pres_mod
 343  interface descomp
 344    module procedure lookup_des_0d, lookup_des_1d, lookup_des_2d, lookup_des_3d
 345  end interface
 346
 347 !-----------------------------------------------------------------------
 348
 349 ! <INTERFACE NAME="lookup_es_des">
 350
 351 !   <OVERVIEW>
 352 !     For the given temperatures, returns the saturation vapor pressure
 353 !     and the derivative of saturation vapor pressure with respect to
 354 !     temperature.
 355 !   </OVERVIEW>
 356 !   <DESCRIPTION>
 357 !     For the given temperatures these routines return the
 358 !     saturation vapor pressure (esat) and the derivative of esat w.r.t
 359 !     temperature (desat). The return values are derived from
 360 !     lookup tables (see notes below).
 361 !   </DESCRIPTION>
 362 !   <TEMPLATE>
 363 !     call lookup_es_des( temp, esat, desat, err_msg )
 364 !   </TEMPLATE>
 365 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 366 !     Temperature in degrees Kelvin.
 367 !   </IN>
 368 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 369 !     Saturation vapor pressure in pascals.
 370 !             May be a scalar, 1d, 2d, or 3d array.
 371 !             Must have the same order and size as temp.
 372 !   </OUT>
 373 !   <OUT NAME="desat" UNITS="pascal/ degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 374 !     Derivative of saturation vapor pressure w.r.t. temperature
 375 !                 in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
 376 !                 Must have the same order and size as temp.
 377 !   </OUT>
 378 !   <OUT NAME="err_msg" UNITS="      " TYPE="character">
 379 !     Character string containing error message to be returned to
 380 !     calling routine.
 381 !   </OUT>
 382 !   <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
 383 !     Temperature(s) provided to the saturation vapor pressure lookup
 384 !          are outside the valid range of the lookup table (-160 to 100 deg C).
 385 !          This may be due to a numerical instability in the model.
 386 !          Information should have been printed to standard output to help
 387 !          determine where the instability may have occurred.
 388 !          If the lookup table needs a larger temperature range,
 389 !          then parameters in the module header must be modified.
 390 !   </ERROR> *
 391
 392  !> @brief For the given temperatures, returns the saturation vapor pressure
 393  !! and the derivative of saturation vapor pressure with respect to
 394  !! temperature.
 395  !!
 396  !> For the given temperatures these routines return the
 397  !! saturation vapor pressure (esat) and the derivative of esat w.r.t
 398  !! temperature (desat). The return values are derived from
 399  !! lookup tables (see notes below).
 400  !!
 401  !! <br>Example usage:
 402  !! @code{.F90} call lookup_es_des( temp, esat, desat, err_msg ) @endcode
 403  !!
 404  !! @param temp Temperature in degrees Kelvin.
 405  !! @param [out] esat Saturation vapor pressure in pascals. May be a scalar, 1d, 2d, or 3d array.
 406  !!              Must have the same order and size as temp.
 407  !! @param [out] desat Derivative of saturation vapor pressure w.r.t. temperature
 408  !!                    in pascals/degree. May be a scalar, 1d, 2d, or 3d array.
 409  !!                    Must have the same order and size as temp.
 410  !! @param [out] err_msg Character string containing error message to be returned to
 411  !!                      calling routine.
 412  !! @error FATAL table overflow, nbad=##
 413  !! Temperature(s) provided to the saturation vapor pressure lookup
 414  !! are outside the valid range of the lookup table (-160 to 100 deg C).
 415  !! This may be due to a numerical instability in the model.
 416  !! Information should have been printed to standard output to help
 417  !! determine where the instability may have occurred.
 418  !! If the lookup table needs a larger temperature range,
 419  !! then parameters in the module header must be modified.
 420  !> @ingroup sat_vapor_pres_mod
 421  interface lookup_es_des
 422    module procedure lookup_es_des_0d, lookup_es_des_1d, lookup_es_des_2d, lookup_es_des_3d
 423  end interface
 424
 425  !> @ingroup sat_vapor_pres_mod
 426  interface lookup_es2
 427    module procedure lookup_es2_0d, lookup_es2_1d, lookup_es2_2d, lookup_es2_3d
 428  end interface
 429
 430  !> @ingroup sat_vapor_pres_mod
 431  interface lookup_des2
 432    module procedure lookup_des2_0d, lookup_des2_1d, lookup_des2_2d, lookup_des2_3d
 433  end interface
 434
 435  !> @ingroup sat_vapor_pres_mod
 436  interface lookup_es2_des2
 437    module procedure lookup_es2_des2_0d, lookup_es2_des2_1d, lookup_es2_des2_2d, lookup_es2_des2_3d
 438  end interface
 439
 440  !> @ingroup sat_vapor_pres_mod
 441  interface lookup_es3
 442    module procedure lookup_es3_0d, lookup_es3_1d, lookup_es3_2d, lookup_es3_3d
 443  end interface
 444
 445  !> @ingroup sat_vapor_pres_mod
 446  interface lookup_des3
 447    module procedure lookup_des3_0d, lookup_des3_1d, lookup_des3_2d, lookup_des3_3d
 448  end interface
 449
 450  !> @ingroup sat_vapor_pres_mod
 451  interface lookup_es3_des3
 452    module procedure lookup_es3_des3_0d, lookup_es3_des3_1d, lookup_es3_des3_2d, lookup_es3_des3_3d
 453  end interface
 454
 455 !-----------------------------------------------------------------------
 456
 457 ! <INTERFACE NAME="compute_qs">
 458
 459 !   <OVERVIEW>
 460 !     For the given temperatures, pressures and optionally vapor
 461 !     specific humidity, returns the specific humidity at saturation
 462 !     (optionally at relative humidity hc instead of at saturation) and
 463 !     optionally the derivative of saturation specific humidity w.r.t.
 464 !     temperature, and the saturation vapor pressure.
 465 !   </OVERVIEW>
 466 !   <DESCRIPTION>
 467 !     For the input temperature and pressure these routines return the
 468 !     specific humidity (qsat) at saturation (unless optional argument
 469 !     hc is used to specify the relative humidity at which qsat should
 470 !     apply) and, if desired, the derivative of qsat w.r.t temperature
 471 !     (dqsdT) and / or the saturation vapor pressure (esat). If the
 472 !     optional input argument specific humidity (q) is present, the
 473 !     exact expression for qs is used; if q is not present the tradit-
 474 !     ional form (valid at saturation) is used. if the optional qsat
 475 !     derivative argument is present, the derivative of qsat w.r.t.
 476 !     temperature will also be returned, defined consistent with the
 477 !     expression used for qsat. The return values are derived from
 478 !     lookup tables (see notes below).
 479 !   </DESCRIPTION>
 480 !   <TEMPLATE>
 481 !     call compute_qs( temp, press, qsat, q, hc, dqsdT, esat, err_msg )
 482 !   </TEMPLATE>
 483 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 484 !     Temperature in degrees Kelvin.
 485 !   </IN>
 486 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 487 !     Air pressure in Pascals.
 488 !   </IN>
 489 !   <OUT NAME="qsat" UNITS="kg(vapor) / kg(moist air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 490 !     Specific humidity in kg (vapor) / kg (moist air)
 491 !             May be a scalar, 1d, 2d, or 3d array.
 492 !             Must have the same order and size as temp.
 493 !   </OUT>
 494 !   <IN NAME="q" UNIT="kg(vapor) / kg (moist air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 495 !     Vapor specific humidity in kg (vapor) / kg (moist air).
 496 !     If present, exact formulation for qsat and dqsdT will be used.
 497 !   </IN>
 498 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)">
 499 !     Relative humidity at which output variables are desired.
 500 !     If not present, values will apply at saturation.
 501 !   </IN>
 502 !   <OUT NAME="dqsdT" UNITS="kg(vapor) / kg(moist air) / degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 503 !     Derivative of saturation specific humidity w.r.t. temperature
 504 !                 in kg(vapor) / kg(moist air) / degree. May be a
 505 !                 scalar, 1d, 2d, or 3d array.
 506 !                 Must have the same order and size as temp.
 507 !   </OUT>
 508 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 509 !     Saturation vapor pressure. May be a scalar, 1d, 2d, or 3d array.
 510 !                 Must have the same order and size as temp.
 511 !   </OUT>
 512 !   <OUT NAME="err_msg" UNITS="      " TYPE="character">
 513 !     Character string containing error message to be returned to
 514 !     calling routine.
 515 !   </OUT>
 516 !   <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
 517 !     Temperature(s) provided to the saturation vapor pressure lookup
 518 !          are outside the valid range of the lookup table (-160 to 100 deg C).
 519 !          This may be due to a numerical instability in the model.
 520 !          Information should have been printed to standard output to help
 521 !          determine where the instability may have occurred.
 522 !          If the lookup table needs a larger temperature range,
 523 !          then parameters in the module header must be modified.
 524 !   </ERROR> *
 525
 526  !> @brief For the given temperatures, pressures and optionally vapor
 527  !! specific humidity, returns the specific humidity at saturation
 528  !! (optionally at relative humidity hc instead of at saturation) and
 529  !! optionally the derivative of saturation specific humidity w.r.t.
 530  !! temperature, and the saturation vapor pressure.
 531  !!
 532  !! For the input temperature and pressure these routines return the
 533  !! specific humidity (qsat) at saturation (unless optional argument
 534  !! hc is used to specify the relative humidity at which qsat should
 535  !! apply) and, if desired, the derivative of qsat w.r.t temperature
 536  !! (dqsdT) and / or the saturation vapor pressure (esat). If the
 537  !! optional input argument specific humidity (q) is present, the
 538  !! exact expression for qs is used; if q is not present the tradit-
 539  !! ional form (valid at saturation) is used. if the optional qsat
 540  !! derivative argument is present, the derivative of qsat w.r.t.
 541  !! temperature will also be returned, defined consistent with the
 542  !! expression used for qsat. The return values are derived from
 543  !! lookup tables (see notes below).
 544  !!
 545  !! Example usage:
 546  !! @code{.F90} call compute_qs( temp, press, qsat, q, hc, dqsdT, esat, err_msg ) @endcode
 547  !!
 548  !> @ingroup sat_vapor_pres_mod
 549  interface compute_qs
 550    module procedure compute_qs_0d, compute_qs_1d, compute_qs_2d, compute_qs_3d
 551  end interface
 552
 553 !-----------------------------------------------------------------------
 554
 555 ! <INTERFACE NAME="compute_mrs">
 556
 557 !   <OVERVIEW>
 558 !     For the given temperatures, pressures and optionally vapor
 559 !     mixing ratio, returns the  vapor mixing ratio at saturation
 560 !     (optionally at relative humidity hc instead of at saturation) and
 561 !     optionally the derivative of saturation vapor mixing ratio w.r.t.
 562 !     temperature, and the saturation vapor pressure.
 563 !   </OVERVIEW>
 564 !   <DESCRIPTION>
 565 !     For the input temperature and pressure these routines return the
 566 !     vapor mixing ratio (mrsat) at saturation (unless optional argument
 567 !     hc is used to specify the relative humidity at which mrsat should
 568 !     apply) and, if desired, the derivative of mrsat w.r.t temperature
 569 !     (dmrsdT) and / or the saturation vapor pressure (esat). If the
 570 !     optional input argument specific humidity (mr) is present, the
 571 !     exact expression for mrs is used; if qr is not present the tradit-
 572 !     ional form (valid at saturation) is used. if the optional mrsat
 573 !     derivative argument is present, the derivative of mrsat w.r.t.
 574 !     temperature will also be returned, defined consistent with the
 575 !     expression used for mrsat. The return values are derived from
 576 !     lookup tables (see notes below).
 577 !   </DESCRIPTION>
 578 !   <TEMPLATE>
 579 !     call compute_mrs( temp, press, mrsat, mr, hc, dmrsdT, esat,
 580 !                       err_msg )
 581 !   </TEMPLATE>
 582 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 583 !     Temperature in degrees Kelvin.
 584 !   </IN>
 585 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 586 !     Air pressure in Pascals.
 587 !   </IN>
 588 !   <OUT NAME="mrsat" UNITS="kg(vapor) / kg (dry air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 589 !     Vapor mixing ratio in kg (vapor) / kg (dry air)
 590 !             May be a scalar, 1d, 2d, or 3d array.
 591 !             Must have the same order and size as temp.
 592 !   </OUT>
 593 !   <IN NAME="mr" UNIT="kg(vapor) / kg (dry air)" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 594 !     Vapor mixing ratio in kg (vapor) / kg (dry air).
 595 !     If present, exact formulation for mrsat and dmrsdT will be used.
 596 !   </IN>
 597 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)">
 598 !     Relative humidity at which output variables are desired.
 599 !     If not present, values will apply at saturation.
 600 !   </IN>
 601 !   <OUT NAME="dmrsdT" UNITS="kg(vapor) / kg(dry air) / degree" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 602 !     Derivative of saturation vapor mixing ratio w.r.t. temperature
 603 !                 in kg(vapor) / kg(dry air) / degree. May be a
 604 !                 scalar, 1d, 2d, or 3d array.
 605 !                 Must have the same order and size as temp.
 606 !   </OUT>
 607 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 608 !     Saturation vapor pressure. May be a scalar, 1d, 2d, or 3d array.
 609 !                 Must have the same order and size as temp.
 610 !   </OUT>
 611 !   <OUT NAME="err_msg" UNITS="      " TYPE="character">
 612 !     Character string containing error message to be returned to
 613 !     calling routine.
 614 !   </OUT>
 615 !   <ERROR MSG="table overflow, nbad=##" STATUS="FATAL">
 616 !     Temperature(s) provided to the saturation vapor pressure lookup
 617 !          are outside the valid range of the lookup table (-160 to 100 deg C).
 618 !          This may be due to a numerical instability in the model.
 619 !          Information should have been printed to standard output to help
 620 !          determine where the instability may have occurred.
 621 !          If the lookup table needs a larger temperature range,
 622 !          then parameters in the module header must be modified.
 623 !   </ERROR> *
 624
 625  !> For the given temperatures, pressures and optionally vapor
 626  !! mixing ratio, returns the  vapor mixing ratio at saturation
 627  !! (optionally at relative humidity hc instead of at saturation) and
 628  !! optionally the derivative of saturation vapor mixing ratio w.r.t.
 629  !! temperature, and the saturation vapor pressure.
 630  !!
 631  !! For the input temperature and pressure these routines return the
 632  !! vapor mixing ratio (mrsat) at saturation (unless optional argument
 633  !! hc is used to specify the relative humidity at which mrsat should
 634  !! apply) and, if desired, the derivative of mrsat w.r.t temperature
 635  !! (dmrsdT) and / or the saturation vapor pressure (esat). If the
 636  !! optional input argument specific humidity (mr) is present, the
 637  !! exact expression for mrs is used; if qr is not present the tradit-
 638  !! ional form (valid at saturation) is used. if the optional mrsat
 639  !! derivative argument is present, the derivative of mrsat w.r.t.
 640  !! temperature will also be returned, defined consistent with the
 641  !! expression used for mrsat. The return values are derived from
 642  !! lookup tables (see notes below).
 643  !!
 644  !! <br>Example usage:
 645  !! @code{.F90} call compute_mrs( temp, press, mrsat, mr, hc, dmrsdT, esat,
 646  !!                       err_msg ) @endcode
 647  !> @ingroup sat_vapor_pres_mod
 648  interface compute_mrs
 649    module procedure compute_mrs_0d, compute_mrs_1d, compute_mrs_2d, compute_mrs_3d
 650  end interface
 651
 652 !-----------------------------------------------------------------------
 653 ! <INTERFACE NAME="compute_es">
 654
 655 !   <OVERVIEW>
 656 !     For the given temperatures, computes the saturation vapor pressures.
 657 !   </OVERVIEW>
 658 !   <DESCRIPTION>
 659 !     Computes saturation vapor pressure for the given temperature using
 660 !     the equations given in the Smithsonian Meteorological Tables.
 661 !     Between -20C and 0C a blended value over ice and water is returned.
 662 !   </DESCRIPTION>
 663 !   <TEMPLATE>
 664 !     es = compute_es ( temp )
 665 !   </TEMPLATE>
 666 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 667 !     Temperature in degrees Kelvin.
 668 !   </IN>
 669 !   <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(scalar),(:),(:,:),(:,:,:)">
 670 !     Saturation vapor pressure in pascals.
 671 !             May be a scalar, 1d, 2d, or 3d array.
 672 !             Must have the same order and size as temp.
 673 !   </OUT>
 674
 675 !interface compute_es
 676 !  module procedure compute_es_0d, compute_es_1d, compute_es_2d, compute_es_3d
 677 !end interface
 678 ! </INTERFACE>
 679 !-----------------------------------------------------------------------
 680  !> @ingroup sat_vapor_pres_mod
 681  interface temp_check
 682    module procedure temp_check_1d, temp_check_2d, temp_check_3d
 683  end interface
 684
 685  !> @ingroup sat_vapor_pres_mod
 686  interface show_all_bad
 687    module procedure show_all_bad_0d, show_all_bad_1d, show_all_bad_2d, show_all_bad_3d
 688  end interface
 689
 690 !> @addtogroup sat_vapor_pres_mod
 691 !> @{
 692 !-----------------------------------------------------------------------
 693 ! Include variable "version" to be written to log file.
 694 #include<file_version.h>
 695
 696  logical :: module_is_initialized = .false.
 697
 698 !-----------------------------------------------------------------------
 699 !  parameters for use in computing qs and mrs
 700
 701  real, parameter    :: EPSILO = RDGAS/RVGAS
 702  real, parameter    :: ZVIR = RVGAS/RDGAS - 1.0
 703
 704 !-----------------------------------------------------------------------
 705 !  parameters for table size and resolution
 706
 707  integer, public :: tcmin = -160  ! minimum temperature (degC) in lookup table
 708  integer, public :: tcmax =  100  ! maximum temperature (degC) in lookup table
 709  integer :: esres =  10   ! table resolution (increments per degree)
 710  integer :: nsize  ! (tcmax-tcmin)*esres+1    !  lookup table size
 711  integer :: nlim   ! nsize-1
 712
 713  integer :: stdoutunit=0
 714 !-----------------------------------------------------------------------
 715 !  variables needed by temp_check
 716  real :: tmin, dtinv, teps
 717
 718 ! The default values below preserve the behavior of omsk and earlier revisions.
 719  logical :: show_bad_value_count_by_slice=.true.
 720  logical :: show_all_bad_values=.false.
 721  logical :: use_exact_qs = .false.
 722  logical :: do_simple             =.false.
 723  logical :: construct_table_wrt_liq = .false.
 724  logical :: construct_table_wrt_liq_and_ice = .false.
 725
 726  namelist / sat_vapor_pres_nml / show_bad_value_count_by_slice, show_all_bad_values, &
 727                                  use_exact_qs, do_simple, &
 728                                  construct_table_wrt_liq, &
 729                                  construct_table_wrt_liq_and_ice
 730
 731 contains
 732
 733 !#######################################################################
 734 ! <SUBROUTINE NAME="lookup_es_0d" INTERFACE="lookup_es">
 735 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
 736 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
 737 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 738 ! </SUBROUTINE>
 739  subroutine lookup_es_0d ( temp, esat, err_msg )
 740
 741  class(*), intent(in)  :: temp
 742  class(*), intent(out) :: esat
 743  character(len=*), intent(out), optional :: err_msg
 744
 745  integer :: nbad
 746  character(len=128) :: err_msg_local
 747
 748    if (.not.module_is_initialized) then
 749       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 750    endif
 751
 752    call lookup_es_k(temp, esat, nbad)
 753
 754    if ( nbad == 0 ) then
 755      if(present(err_msg)) err_msg = ''
 756    else
 757      if(show_all_bad_values) call show_all_bad ( temp )
 758      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 759      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
 760    endif
 761
 762  end subroutine lookup_es_0d
 763
 764 !#######################################################################
 765
 766 ! <SUBROUTINE NAME="lookup_es_1d" INTERFACE="lookup_es">
 767 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
 768 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
 769 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 770 ! </SUBROUTINE>
 771  subroutine lookup_es_1d ( temp, esat, err_msg )
 772
 773  class(*), intent(in)  :: temp(:)
 774  class(*), intent(out) :: esat(:)
 775  character(len=*), intent(out), optional :: err_msg
 776
 777  character(len=54) :: err_msg_local
 778  integer :: nbad
 779 !-----------------------------------------------
 780
 781    if (.not.module_is_initialized) then
 782       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 783    endif
 784
 785    call lookup_es_k(temp, esat, nbad)
 786
 787    if ( nbad == 0 ) then
 788      if(present(err_msg)) err_msg = ''
 789    else
 790      if(show_bad_value_count_by_slice) call temp_check ( temp )
 791      if(show_all_bad_values) call show_all_bad ( temp )
 792      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 793      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
 794    endif
 795
 796 !-----------------------------------------------
 797
 798  end subroutine lookup_es_1d
 799
 800 !#######################################################################
 801
 802 ! <SUBROUTINE NAME="lookup_es_2d" INTERFACE="lookup_es">
 803 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
 804 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
 805 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 806 ! </SUBROUTINE>
 807  subroutine lookup_es_2d ( temp, esat, err_msg )
 808
 809  class(*), intent(in)  :: temp(:,:)
 810  class(*), intent(out) :: esat(:,:)
 811  character(len=*), intent(out), optional :: err_msg
 812
 813  character(len=54) :: err_msg_local
 814  integer :: nbad
 815 !-----------------------------------------------
 816
 817    if (.not.module_is_initialized) then
 818       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 819    endif
 820
 821    call lookup_es_k(temp, esat, nbad)
 822
 823    if ( nbad == 0 ) then
 824      if(present(err_msg)) err_msg = ''
 825    else
 826      if(show_bad_value_count_by_slice) call temp_check ( temp )
 827      if(show_all_bad_values) call show_all_bad ( temp )
 828      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 829      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
 830    endif
 831
 832 !-----------------------------------------------
 833
 834  end subroutine lookup_es_2d
 835
 836 !#######################################################################
 837
 838 ! <SUBROUTINE NAME="lookup_es_3d" INTERFACE="lookup_es">
 839 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
 840 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
 841 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 842 ! </SUBROUTINE>
 843  subroutine lookup_es_3d ( temp, esat, err_msg )
 844
 845  class(*), intent(in)  :: temp(:,:,:)
 846  class(*), intent(out) :: esat(:,:,:)
 847  character(len=*), intent(out), optional :: err_msg
 848
 849  integer :: nbad
 850  character(len=128) :: err_msg_tmp
 851
 852    if (.not.module_is_initialized) then
 853       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 854    endif
 855
 856    call lookup_es_k(temp, esat, nbad)
 857
 858    if ( nbad == 0 ) then
 859      if(present(err_msg)) err_msg = ''
 860    else
 861      if(show_bad_value_count_by_slice) call temp_check ( temp )
 862      if(show_all_bad_values) call show_all_bad ( temp )
 863      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 864      if(fms_error_handler('lookup_es',err_msg_tmp,err_msg)) return
 865    endif
 866
 867  end subroutine lookup_es_3d
 868
 869
 870 !#######################################################################
 871 ! <SUBROUTINE NAME="lookup_es2_0d" INTERFACE="lookup_es2">
 872 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
 873 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
 874 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 875 ! </SUBROUTINE>
 876  subroutine lookup_es2_0d ( temp, esat, err_msg )
 877
 878  real, intent(in)  :: temp
 879  real, intent(out) :: esat
 880  character(len=*), intent(out), optional :: err_msg
 881
 882  integer :: nbad
 883  character(len=128) :: err_msg_local
 884
 885    if (.not.module_is_initialized) then
 886       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 887    endif
 888
 889    call lookup_es2_k(temp, esat, nbad)
 890
 891    if ( nbad == 0 ) then
 892      if(present(err_msg)) err_msg = ''
 893    else
 894      if(show_all_bad_values) call show_all_bad ( temp )
 895      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 896      if(fms_error_handler('lookup_es2',err_msg_local,err_msg)) return
 897    endif
 898
 899  end subroutine lookup_es2_0d
 900
 901 !#######################################################################
 902
 903 ! <SUBROUTINE NAME="lookup_es2_1d" INTERFACE="lookup_es2">
 904 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
 905 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
 906 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 907 ! </SUBROUTINE>
 908  subroutine lookup_es2_1d ( temp, esat, err_msg )
 909
 910  real, intent(in)  :: temp(:)
 911  real, intent(out) :: esat(:)
 912  character(len=*), intent(out), optional :: err_msg
 913
 914  character(len=54) :: err_msg_local
 915  integer :: nbad
 916 !-----------------------------------------------
 917
 918    if (.not.module_is_initialized) then
 919       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 920    endif
 921
 922    call lookup_es2_k(temp, esat, nbad)
 923
 924    if ( nbad == 0 ) then
 925      if(present(err_msg)) err_msg = ''
 926    else
 927      if(show_bad_value_count_by_slice) call temp_check ( temp )
 928      if(show_all_bad_values) call show_all_bad ( temp )
 929      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 930      if(fms_error_handler('lookup_es2',err_msg_local,err_msg)) return
 931    endif
 932
 933 !-----------------------------------------------
 934
 935  end subroutine lookup_es2_1d
 936
 937 !#######################################################################
 938
 939 ! <SUBROUTINE NAME="lookup_es2_2d" INTERFACE="lookup_es2">
 940 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
 941 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
 942 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 943 ! </SUBROUTINE>
 944  subroutine lookup_es2_2d ( temp, esat, err_msg )
 945
 946  real, intent(in)  :: temp(:,:)
 947  real, intent(out) :: esat(:,:)
 948  character(len=*), intent(out), optional :: err_msg
 949
 950  character(len=54) :: err_msg_local
 951  integer :: nbad
 952 !-----------------------------------------------
 953
 954    if (.not.module_is_initialized) then
 955       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 956    endif
 957
 958    call lookup_es2_k(temp, esat, nbad)
 959
 960    if ( nbad == 0 ) then
 961      if(present(err_msg)) err_msg = ''
 962    else
 963      if(show_bad_value_count_by_slice) call temp_check ( temp )
 964      if(show_all_bad_values) call show_all_bad ( temp )
 965      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
 966      if(fms_error_handler('lookup_es2',err_msg_local,err_msg)) return
 967    endif
 968
 969 !-----------------------------------------------
 970
 971  end subroutine lookup_es2_2d
 972
 973 !#######################################################################
 974
 975 ! <SUBROUTINE NAME="lookup_es2_3d" INTERFACE="lookup_es2">
 976 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
 977 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
 978 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
 979 ! </SUBROUTINE>
 980  subroutine lookup_es2_3d ( temp, esat, err_msg )
 981
 982  real, intent(in)  :: temp(:,:,:)
 983  real, intent(out) :: esat(:,:,:)
 984  character(len=*), intent(out), optional :: err_msg
 985
 986  integer :: nbad
 987  character(len=128) :: err_msg_tmp
 988
 989    if (.not.module_is_initialized) then
 990       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
 991    endif
 992
 993    call lookup_es2_k(temp, esat, nbad)
 994
 995    if ( nbad == 0 ) then
 996      if(present(err_msg)) err_msg = ''
 997    else
 998      if(show_bad_value_count_by_slice) call temp_check ( temp )
 999      if(show_all_bad_values) call show_all_bad ( temp )
1000      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1001      if(fms_error_handler('lookup_es2',err_msg_tmp,err_msg)) return
1002    endif
1003
1004  end subroutine lookup_es2_3d
1005
1006
1007 !#######################################################################
1008 ! <SUBROUTINE NAME="lookup_es3_0d" INTERFACE="lookup_es3">
1009 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1010 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1011 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1012 ! </SUBROUTINE>
1013  subroutine lookup_es3_0d ( temp, esat, err_msg )
1014
1015  real, intent(in)  :: temp
1016  real, intent(out) :: esat
1017  character(len=*), intent(out), optional :: err_msg
1018
1019  integer :: nbad
1020  character(len=128) :: err_msg_local
1021
1022    if (.not.module_is_initialized) then
1023       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1024    endif
1025
1026    call lookup_es3_k(temp, esat, nbad)
1027
1028    if ( nbad == 0 ) then
1029      if(present(err_msg)) err_msg = ''
1030    else
1031      if(show_all_bad_values) call show_all_bad ( temp )
1032      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1033      if(fms_error_handler('lookup_es3',err_msg_local,err_msg)) return
1034    endif
1035
1036  end subroutine lookup_es3_0d
1037
1038 !#######################################################################
1039
1040 ! <SUBROUTINE NAME="lookup_es3_1d" INTERFACE="lookup_es3">
1041 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1042 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1043 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1044 ! </SUBROUTINE>
1045  subroutine lookup_es3_1d ( temp, esat, err_msg )
1046
1047  real, intent(in)  :: temp(:)
1048  real, intent(out) :: esat(:)
1049  character(len=*), intent(out), optional :: err_msg
1050
1051  character(len=54) :: err_msg_local
1052  integer :: nbad
1053 !-----------------------------------------------
1054
1055    if (.not.module_is_initialized) then
1056       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1057    endif
1058
1059    call lookup_es3_k(temp, esat, nbad)
1060
1061    if ( nbad == 0 ) then
1062      if(present(err_msg)) err_msg = ''
1063    else
1064      if(show_bad_value_count_by_slice) call temp_check ( temp )
1065      if(show_all_bad_values) call show_all_bad ( temp )
1066      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1067      if(fms_error_handler('lookup_es3',err_msg_local,err_msg)) return
1068    endif
1069
1070 !-----------------------------------------------
1071
1072  end subroutine lookup_es3_1d
1073
1074 !#######################################################################
1075
1076 ! <SUBROUTINE NAME="lookup_es3_2d" INTERFACE="lookup_es3">
1077 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1078 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1079 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1080 ! </SUBROUTINE>
1081  subroutine lookup_es3_2d ( temp, esat, err_msg )
1082
1083  real, intent(in)  :: temp(:,:)
1084  real, intent(out) :: esat(:,:)
1085  character(len=*), intent(out), optional :: err_msg
1086
1087  character(len=54) :: err_msg_local
1088  integer :: nbad
1089 !-----------------------------------------------
1090
1091    if (.not.module_is_initialized) then
1092       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1093    endif
1094
1095    call lookup_es3_k(temp, esat, nbad)
1096
1097    if ( nbad == 0 ) then
1098      if(present(err_msg)) err_msg = ''
1099    else
1100      if(show_bad_value_count_by_slice) call temp_check ( temp )
1101      if(show_all_bad_values) call show_all_bad ( temp )
1102      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1103      if(fms_error_handler('lookup_es3',err_msg_local,err_msg)) return
1104    endif
1105
1106 !-----------------------------------------------
1107
1108  end subroutine lookup_es3_2d
1109
1110 !#######################################################################
1111
1112 ! <SUBROUTINE NAME="lookup_es3_3d" INTERFACE="lookup_es3">
1113 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1114 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1115 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1116 ! </SUBROUTINE>
1117  subroutine lookup_es3_3d ( temp, esat, err_msg )
1118
1119  real, intent(in)  :: temp(:,:,:)
1120  real, intent(out) :: esat(:,:,:)
1121  character(len=*), intent(out), optional :: err_msg
1122
1123  integer :: nbad
1124  character(len=128) :: err_msg_tmp
1125
1126    if (.not.module_is_initialized) then
1127       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1128    endif
1129
1130    call lookup_es3_k(temp, esat, nbad)
1131
1132    if ( nbad == 0 ) then
1133      if(present(err_msg)) err_msg = ''
1134    else
1135      if(show_bad_value_count_by_slice) call temp_check ( temp )
1136      if(show_all_bad_values) call show_all_bad ( temp )
1137      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1138      if(fms_error_handler('lookup_es3',err_msg_tmp,err_msg)) return
1139    endif
1140
1141  end subroutine lookup_es3_3d
1142
1143
1144 !#######################################################################
1145 !  routines for computing derivative of es
1146 !#######################################################################
1147
1148 ! <SUBROUTINE NAME="lookup_des_0d" INTERFACE="lookup_des">
1149 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1150 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1151 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1152 ! </SUBROUTINE>
1153  subroutine lookup_des_0d ( temp, desat, err_msg )
1154
1155  real, intent(in)  :: temp
1156  real, intent(out) :: desat
1157  character(len=*), intent(out), optional :: err_msg
1158
1159  integer :: nbad
1160  character(len=128) :: err_msg_local
1161
1162    if (.not.module_is_initialized) then
1163       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1164    endif
1165
1166    call lookup_des_k( temp, desat, nbad)
1167
1168    if ( nbad == 0 ) then
1169      if(present(err_msg)) err_msg = ''
1170    else
1171      if(show_all_bad_values) call show_all_bad ( temp )
1172      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1173      if(fms_error_handler('lookup_des',err_msg_local,err_msg)) return
1174    endif
1175
1176  end subroutine lookup_des_0d
1177
1178 !#######################################################################
1179
1180 ! <SUBROUTINE NAME="lookup_des_1d" INTERFACE="lookup_des">
1181 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1182 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1183 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1184 ! </SUBROUTINE>
1185  subroutine lookup_des_1d ( temp, desat, err_msg )
1186
1187  real, intent(in)  :: temp (:)
1188  real, intent(out) :: desat(:)
1189  character(len=*), intent(out), optional :: err_msg
1190
1191  character(len=54) :: err_msg_local
1192  integer :: nbad
1193 !-----------------------------------------------
1194
1195    if (.not.module_is_initialized) then
1196       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1197    endif
1198
1199    if(present(err_msg)) err_msg=''
1200
1201    call lookup_des_k(temp, desat, nbad)
1202
1203    if ( nbad == 0 ) then
1204      if(present(err_msg)) err_msg = ''
1205    else
1206      if(show_bad_value_count_by_slice) call temp_check ( temp )
1207      if(show_all_bad_values) call show_all_bad ( temp )
1208      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1209      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1210    endif
1211 !-----------------------------------------------
1212
1213  end subroutine lookup_des_1d
1214
1215 !#######################################################################
1216
1217 ! <SUBROUTINE NAME="lookup_des_2d" INTERFACE="lookup_des">
1218 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1219 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1220 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1221 ! </SUBROUTINE>
1222  subroutine lookup_des_2d ( temp, desat, err_msg )
1223
1224  real, intent(in)  :: temp (:,:)
1225  real, intent(out) :: desat(:,:)
1226  character(len=*), intent(out), optional :: err_msg
1227
1228  character(len=54) :: err_msg_local
1229  integer :: nbad
1230 !-----------------------------------------------
1231
1232    if (.not.module_is_initialized) then
1233       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1234    endif
1235
1236    call lookup_des_k(temp, desat, nbad)
1237
1238    if ( nbad == 0 ) then
1239      if(present(err_msg)) err_msg = ''
1240    else
1241      if(show_bad_value_count_by_slice) call temp_check ( temp )
1242      if(show_all_bad_values) call show_all_bad ( temp )
1243      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1244      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1245    endif
1246 !-----------------------------------------------
1247
1248  end subroutine lookup_des_2d
1249
1250 !#######################################################################
1251 ! <SUBROUTINE NAME="lookup_des_3d" INTERFACE="lookup_des">
1252 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1253 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1254 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1255 ! </SUBROUTINE>
1256  subroutine lookup_des_3d ( temp, desat, err_msg )
1257
1258  real, intent(in)  :: temp (:,:,:)
1259  real, intent(out) :: desat(:,:,:)
1260  character(len=*), intent(out), optional :: err_msg
1261
1262  integer :: nbad
1263  character(len=128) :: err_msg_tmp
1264
1265    if (.not.module_is_initialized) then
1266       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1267    endif
1268
1269    call lookup_des_k( temp, desat, nbad )
1270
1271    if ( nbad == 0 ) then
1272      if(present(err_msg)) err_msg=''
1273    else
1274      if(show_bad_value_count_by_slice) call temp_check ( temp )
1275      if(show_all_bad_values) call show_all_bad ( temp )
1276      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1277      if(fms_error_handler('lookup_des',err_msg_tmp,err_msg)) return
1278    endif
1279
1280  end subroutine lookup_des_3d
1281
1282
1283 ! <SUBROUTINE NAME="lookup_des2_0d" INTERFACE="lookup_des2">
1284 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1285 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1286 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1287 ! </SUBROUTINE>
1288  subroutine lookup_des2_0d ( temp, desat, err_msg )
1289
1290  real, intent(in)  :: temp
1291  real, intent(out) :: desat
1292  character(len=*), intent(out), optional :: err_msg
1293
1294  integer :: nbad
1295  character(len=128) :: err_msg_local
1296
1297    if (.not.module_is_initialized) then
1298       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1299    endif
1300
1301    call lookup_des2_k( temp, desat, nbad)
1302
1303    if ( nbad == 0 ) then
1304      if(present(err_msg)) err_msg = ''
1305    else
1306      if(show_all_bad_values) call show_all_bad ( temp )
1307      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1308      if(fms_error_handler('lookup_des2',err_msg_local,err_msg)) return
1309    endif
1310
1311  end subroutine lookup_des2_0d
1312
1313 !#######################################################################
1314
1315 ! <SUBROUTINE NAME="lookup_des2_1d" INTERFACE="lookup_des2">
1316 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1317 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1318 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1319 ! </SUBROUTINE>
1320  subroutine lookup_des2_1d ( temp, desat, err_msg )
1321
1322  real, intent(in)  :: temp (:)
1323  real, intent(out) :: desat(:)
1324  character(len=*), intent(out), optional :: err_msg
1325
1326  character(len=54) :: err_msg_local
1327  integer :: nbad
1328 !-----------------------------------------------
1329
1330    if (.not.module_is_initialized) then
1331       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1332    endif
1333
1334    if(present(err_msg)) err_msg=''
1335
1336    call lookup_des2_k(temp, desat, nbad)
1337
1338    if ( nbad == 0 ) then
1339      if(present(err_msg)) err_msg = ''
1340    else
1341      if(show_bad_value_count_by_slice) call temp_check ( temp )
1342      if(show_all_bad_values) call show_all_bad ( temp )
1343      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1344      if(fms_error_handler('lookup_des2',err_msg_local,err_msg)) return
1345    endif
1346 !-----------------------------------------------
1347
1348  end subroutine lookup_des2_1d
1349
1350 !#######################################################################
1351
1352 ! <SUBROUTINE NAME="lookup_des2_2d" INTERFACE="lookup_des2">
1353 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1354 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1355 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1356 ! </SUBROUTINE>
1357  subroutine lookup_des2_2d ( temp, desat, err_msg )
1358
1359  real, intent(in)  :: temp (:,:)
1360  real, intent(out) :: desat(:,:)
1361  character(len=*), intent(out), optional :: err_msg
1362
1363  character(len=54) :: err_msg_local
1364  integer :: nbad
1365 !-----------------------------------------------
1366
1367    if (.not.module_is_initialized) then
1368       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1369    endif
1370
1371    call lookup_des2_k(temp, desat, nbad)
1372
1373    if ( nbad == 0 ) then
1374      if(present(err_msg)) err_msg = ''
1375    else
1376      if(show_bad_value_count_by_slice) call temp_check ( temp )
1377      if(show_all_bad_values) call show_all_bad ( temp )
1378      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1379      if(fms_error_handler('lookup_des2',err_msg_local,err_msg)) return
1380    endif
1381 !-----------------------------------------------
1382
1383  end subroutine lookup_des2_2d
1384
1385 !#######################################################################
1386 ! <SUBROUTINE NAME="lookup_des2_3d" INTERFACE="lookup_des2">
1387 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1388 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1389 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1390 ! </SUBROUTINE>
1391  subroutine lookup_des2_3d ( temp, desat, err_msg )
1392
1393  real, intent(in)  :: temp (:,:,:)
1394  real, intent(out) :: desat(:,:,:)
1395  character(len=*), intent(out), optional :: err_msg
1396
1397  integer :: nbad
1398  character(len=128) :: err_msg_tmp
1399
1400    if (.not.module_is_initialized) then
1401       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1402    endif
1403
1404    call lookup_des2_k( temp, desat, nbad )
1405
1406    if ( nbad == 0 ) then
1407      if(present(err_msg)) err_msg=''
1408    else
1409      if(show_bad_value_count_by_slice) call temp_check ( temp )
1410      if(show_all_bad_values) call show_all_bad ( temp )
1411      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1412      if(fms_error_handler('lookup_des2',err_msg_tmp,err_msg)) return
1413    endif
1414
1415  end subroutine lookup_des2_3d
1416
1417
1418 ! <SUBROUTINE NAME="lookup_des3_0d" INTERFACE="lookup_des3">
1419 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1420 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1421 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1422 ! </SUBROUTINE>
1423  subroutine lookup_des3_0d ( temp, desat, err_msg )
1424
1425  real, intent(in)  :: temp
1426  real, intent(out) :: desat
1427  character(len=*), intent(out), optional :: err_msg
1428
1429  integer :: nbad
1430  character(len=128) :: err_msg_local
1431
1432    if (.not.module_is_initialized) then
1433       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1434    endif
1435
1436    call lookup_des3_k( temp, desat, nbad)
1437
1438    if ( nbad == 0 ) then
1439      if(present(err_msg)) err_msg = ''
1440    else
1441      if(show_all_bad_values) call show_all_bad ( temp )
1442      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1443      if(fms_error_handler('lookup_des3',err_msg_local,err_msg)) return
1444    endif
1445
1446  end subroutine lookup_des3_0d
1447
1448 !#######################################################################
1449
1450 ! <SUBROUTINE NAME="lookup_des3_1d" INTERFACE="lookup_des3">
1451 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1452 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1453 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1454 ! </SUBROUTINE>
1455  subroutine lookup_des3_1d ( temp, desat, err_msg )
1456
1457  real, intent(in)  :: temp (:)
1458  real, intent(out) :: desat(:)
1459  character(len=*), intent(out), optional :: err_msg
1460
1461  character(len=54) :: err_msg_local
1462  integer :: nbad
1463 !-----------------------------------------------
1464
1465    if (.not.module_is_initialized) then
1466       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1467    endif
1468
1469    if(present(err_msg)) err_msg=''
1470
1471    call lookup_des3_k(temp, desat, nbad)
1472
1473    if ( nbad == 0 ) then
1474      if(present(err_msg)) err_msg = ''
1475    else
1476      if(show_bad_value_count_by_slice) call temp_check ( temp )
1477      if(show_all_bad_values) call show_all_bad ( temp )
1478      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1479      if(fms_error_handler('lookup_des3',err_msg_local,err_msg)) return
1480    endif
1481 !-----------------------------------------------
1482
1483  end subroutine lookup_des3_1d
1484
1485 !#######################################################################
1486
1487 ! <SUBROUTINE NAME="lookup_des3_2d" INTERFACE="lookup_des3">
1488 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1489 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1490 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1491 ! </SUBROUTINE>
1492  subroutine lookup_des3_2d ( temp, desat, err_msg )
1493
1494  real, intent(in)  :: temp (:,:)
1495  real, intent(out) :: desat(:,:)
1496  character(len=*), intent(out), optional :: err_msg
1497
1498  character(len=54) :: err_msg_local
1499  integer :: nbad
1500 !-----------------------------------------------
1501
1502    if (.not.module_is_initialized) then
1503       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1504    endif
1505
1506    call lookup_des3_k(temp, desat, nbad)
1507
1508    if ( nbad == 0 ) then
1509      if(present(err_msg)) err_msg = ''
1510    else
1511      if(show_bad_value_count_by_slice) call temp_check ( temp )
1512      if(show_all_bad_values) call show_all_bad ( temp )
1513      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1514      if(fms_error_handler('lookup_des3',err_msg_local,err_msg)) return
1515    endif
1516 !-----------------------------------------------
1517
1518  end subroutine lookup_des3_2d
1519
1520 !#######################################################################
1521 ! <SUBROUTINE NAME="lookup_des3_3d" INTERFACE="lookup_des3">
1522 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1523 !   <OUT NAME="desat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1524 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1525 ! </SUBROUTINE>
1526  subroutine lookup_des3_3d ( temp, desat, err_msg )
1527
1528  real, intent(in)  :: temp (:,:,:)
1529  real, intent(out) :: desat(:,:,:)
1530  character(len=*), intent(out), optional :: err_msg
1531
1532  integer :: nbad
1533  character(len=128) :: err_msg_tmp
1534
1535    if (.not.module_is_initialized) then
1536       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1537    endif
1538
1539    call lookup_des3_k( temp, desat, nbad )
1540
1541    if ( nbad == 0 ) then
1542      if(present(err_msg)) err_msg=''
1543    else
1544      if(show_bad_value_count_by_slice) call temp_check ( temp )
1545      if(show_all_bad_values) call show_all_bad ( temp )
1546      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1547      if(fms_error_handler('lookup_des3',err_msg_tmp,err_msg)) return
1548    endif
1549
1550  end subroutine lookup_des3_3d
1551
1552 !========================================================================================================
1553
1554 !#######################################################################
1555
1556 ! <SUBROUTINE NAME="lookup_es_des_0d" INTERFACE="lookup_es_des">
1557 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1558 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1559 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(scalar)"></OUT>
1560 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1561 ! </SUBROUTINE>
1562  subroutine lookup_es_des_0d ( temp, esat, desat, err_msg )
1563
1564  real, intent(in)  :: temp
1565  real, intent(out) :: esat, desat
1566  character(len=*), intent(out), optional :: err_msg
1567
1568  integer :: nbad
1569  character(len=128) :: err_msg_local
1570
1571    if (.not.module_is_initialized) then
1572       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1573    endif
1574
1575    call lookup_es_des_k(temp, esat, desat, nbad)
1576
1577    if ( nbad == 0 ) then
1578      if(present(err_msg)) err_msg = ''
1579    else
1580      if(show_all_bad_values) call show_all_bad ( temp )
1581      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1582      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1583    endif
1584
1585  end subroutine lookup_es_des_0d
1586
1587 !#######################################################################
1588
1589 ! <SUBROUTINE NAME="lookup_es_des_1d" INTERFACE="lookup_es_des">
1590 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1591 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1592 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:)"></OUT>
1593 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1594 ! </SUBROUTINE>
1595  subroutine lookup_es_des_1d ( temp, esat, desat, err_msg )
1596
1597  real, dimension(:), intent(in)  :: temp
1598  real, dimension(:), intent(out) :: esat, desat
1599  character(len=*), intent(out), optional :: err_msg
1600
1601  integer :: nbad
1602  character(len=128) :: err_msg_local
1603
1604    if (.not.module_is_initialized) then
1605       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1606    endif
1607
1608    call lookup_es_des_k(temp, esat, desat, nbad)
1609
1610    if ( nbad == 0 ) then
1611      if(present(err_msg)) err_msg = ''
1612    else
1613      if(show_bad_value_count_by_slice) call temp_check ( temp )
1614      if(show_all_bad_values) call show_all_bad ( temp )
1615      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1616      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1617    endif
1618
1619  end subroutine lookup_es_des_1d
1620
1621 !#######################################################################
1622
1623 ! <SUBROUTINE NAME="lookup_es_des_2d" INTERFACE="lookup_es_des">
1624 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1625 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1626 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:)"></OUT>
1627 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1628 ! </SUBROUTINE>
1629  subroutine lookup_es_des_2d ( temp, esat, desat, err_msg )
1630
1631  real, dimension(:,:), intent(in)  :: temp
1632  real, dimension(:,:), intent(out) :: esat, desat
1633  character(len=*), intent(out), optional :: err_msg
1634
1635  integer :: nbad
1636  character(len=128) :: err_msg_local
1637
1638    if (.not.module_is_initialized) then
1639       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1640    endif
1641
1642    call lookup_es_des_k(temp, esat, desat, nbad)
1643
1644    if ( nbad == 0 ) then
1645      if(present(err_msg)) err_msg = ''
1646    else
1647      if(show_bad_value_count_by_slice) call temp_check ( temp )
1648      if(show_all_bad_values) call show_all_bad ( temp )
1649      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1650      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1651    endif
1652
1653  end subroutine lookup_es_des_2d
1654
1655 !#######################################################################
1656
1657 ! <SUBROUTINE NAME="lookup_es_des_3d" INTERFACE="lookup_es_des">
1658 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1659 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1660 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:,:)"></OUT>
1661 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1662 ! </SUBROUTINE>
1663  subroutine lookup_es_des_3d ( temp, esat, desat, err_msg )
1664
1665  real, dimension(:,:,:), intent(in)  :: temp
1666  real, dimension(:,:,:), intent(out) :: esat, desat
1667  character(len=*), intent(out), optional :: err_msg
1668
1669  integer :: nbad
1670  character(len=128) :: err_msg_local
1671
1672    if (.not.module_is_initialized) then
1673       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1674    endif
1675
1676    call lookup_es_des_k(temp, esat, desat, nbad)
1677
1678    if ( nbad == 0 ) then
1679      if(present(err_msg)) err_msg = ''
1680    else
1681      if(show_bad_value_count_by_slice) call temp_check ( temp )
1682      if(show_all_bad_values) call show_all_bad ( temp )
1683      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1684      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
1685    endif
1686
1687  end subroutine lookup_es_des_3d
1688
1689 !#######################################################################
1690 !#######################################################################
1691
1692 ! <SUBROUTINE NAME="lookup_es2_des2_0d" INTERFACE="lookup_es2_des2">
1693 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1694 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1695 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(scalar)"></OUT>
1696 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1697 ! </SUBROUTINE>
1698  subroutine lookup_es2_des2_0d ( temp, esat, desat, err_msg )
1699
1700  real, intent(in)  :: temp
1701  real, intent(out) :: esat, desat
1702  character(len=*), intent(out), optional :: err_msg
1703
1704  integer :: nbad
1705  character(len=128) :: err_msg_local
1706
1707    if (.not.module_is_initialized) then
1708       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1709    endif
1710
1711    call lookup_es2_des2_k(temp, esat, desat, nbad)
1712
1713    if ( nbad == 0 ) then
1714      if(present(err_msg)) err_msg = ''
1715    else
1716      if(show_all_bad_values) call show_all_bad ( temp )
1717      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1718      if(fms_error_handler('lookup_es2_des2',err_msg_local,err_msg)) return
1719    endif
1720
1721  end subroutine lookup_es2_des2_0d
1722
1723 !#######################################################################
1724
1725 ! <SUBROUTINE NAME="lookup_es2_des2_1d" INTERFACE="lookup_es2_des2">
1726 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1727 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1728 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:)"></OUT>
1729 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1730 ! </SUBROUTINE>
1731  subroutine lookup_es2_des2_1d ( temp, esat, desat, err_msg )
1732
1733  real, dimension(:), intent(in)  :: temp
1734  real, dimension(:), intent(out) :: esat, desat
1735  character(len=*), intent(out), optional :: err_msg
1736
1737  integer :: nbad
1738  character(len=128) :: err_msg_local
1739
1740    if (.not.module_is_initialized) then
1741       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1742    endif
1743
1744    call lookup_es2_des2_k(temp, esat, desat, nbad)
1745
1746    if ( nbad == 0 ) then
1747      if(present(err_msg)) err_msg = ''
1748    else
1749      if(show_bad_value_count_by_slice) call temp_check ( temp )
1750      if(show_all_bad_values) call show_all_bad ( temp )
1751      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1752      if(fms_error_handler('lookup_es2_des2',err_msg_local,err_msg)) return
1753    endif
1754
1755  end subroutine lookup_es2_des2_1d
1756
1757 !#######################################################################
1758
1759 ! <SUBROUTINE NAME="lookup_es2_des2_2d" INTERFACE="lookup_es2_des2">
1760 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1761 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1762 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:)"></OUT>
1763 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1764 ! </SUBROUTINE>
1765  subroutine lookup_es2_des2_2d ( temp, esat, desat, err_msg )
1766
1767  real, dimension(:,:), intent(in)  :: temp
1768  real, dimension(:,:), intent(out) :: esat, desat
1769  character(len=*), intent(out), optional :: err_msg
1770
1771  integer :: nbad
1772  character(len=128) :: err_msg_local
1773
1774    if (.not.module_is_initialized) then
1775       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1776    endif
1777
1778    call lookup_es2_des2_k(temp, esat, desat, nbad)
1779
1780    if ( nbad == 0 ) then
1781      if(present(err_msg)) err_msg = ''
1782    else
1783      if(show_bad_value_count_by_slice) call temp_check ( temp )
1784      if(show_all_bad_values) call show_all_bad ( temp )
1785      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1786      if(fms_error_handler('lookup_es2_des2',err_msg_local,err_msg)) return
1787    endif
1788
1789  end subroutine lookup_es2_des2_2d
1790
1791 !#######################################################################
1792
1793 ! <SUBROUTINE NAME="lookup_es2_des2_3d" INTERFACE="lookup_es2_des2">
1794 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1795 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1796 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:,:)"></OUT>
1797 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1798 ! </SUBROUTINE>
1799  subroutine lookup_es2_des2_3d ( temp, esat, desat, err_msg )
1800
1801  real, dimension(:,:,:), intent(in)  :: temp
1802  real, dimension(:,:,:), intent(out) :: esat, desat
1803  character(len=*), intent(out), optional :: err_msg
1804
1805  integer :: nbad
1806  character(len=128) :: err_msg_local
1807
1808    if (.not.module_is_initialized) then
1809       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1810    endif
1811
1812    call lookup_es2_des2_k(temp, esat, desat, nbad)
1813
1814    if ( nbad == 0 ) then
1815      if(present(err_msg)) err_msg = ''
1816    else
1817      if(show_bad_value_count_by_slice) call temp_check ( temp )
1818      if(show_all_bad_values) call show_all_bad ( temp )
1819      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1820      if(fms_error_handler('lookup_es2_des2',err_msg_local,err_msg)) return
1821    endif
1822
1823  end subroutine lookup_es2_des2_3d
1824
1825
1826 !#######################################################################
1827 !#######################################################################
1828
1829 ! <SUBROUTINE NAME="lookup_es3_des3_0d" INTERFACE="lookup_es3_des3">
1830 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
1831 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
1832 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(scalar)"></OUT>
1833 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1834 ! </SUBROUTINE>
1835  subroutine lookup_es3_des3_0d ( temp, esat, desat, err_msg )
1836
1837  real, intent(in)  :: temp
1838  real, intent(out) :: esat, desat
1839  character(len=*), intent(out), optional :: err_msg
1840
1841  integer :: nbad
1842  character(len=128) :: err_msg_local
1843
1844    if (.not.module_is_initialized) then
1845       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1846    endif
1847
1848    call lookup_es3_des3_k(temp, esat, desat, nbad)
1849
1850    if ( nbad == 0 ) then
1851      if(present(err_msg)) err_msg = ''
1852    else
1853      if(show_all_bad_values) call show_all_bad ( temp )
1854      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1855      if(fms_error_handler('lookup_es3_des3',err_msg_local,err_msg)) return
1856    endif
1857
1858  end subroutine lookup_es3_des3_0d
1859
1860 !#######################################################################
1861
1862 ! <SUBROUTINE NAME="lookup_es3_des3_1d" INTERFACE="lookup_es3_des3">
1863 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
1864 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
1865 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:)"></OUT>
1866 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1867 ! </SUBROUTINE>
1868  subroutine lookup_es3_des3_1d ( temp, esat, desat, err_msg )
1869
1870  real, dimension(:), intent(in)  :: temp
1871  real, dimension(:), intent(out) :: esat, desat
1872  character(len=*), intent(out), optional :: err_msg
1873
1874  integer :: nbad
1875  character(len=128) :: err_msg_local
1876
1877    if (.not.module_is_initialized) then
1878       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1879    endif
1880
1881    call lookup_es3_des3_k(temp, esat, desat, nbad)
1882
1883    if ( nbad == 0 ) then
1884      if(present(err_msg)) err_msg = ''
1885    else
1886      if(show_bad_value_count_by_slice) call temp_check ( temp )
1887      if(show_all_bad_values) call show_all_bad ( temp )
1888      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1889      if(fms_error_handler('lookup_es3_des3',err_msg_local,err_msg)) return
1890    endif
1891
1892  end subroutine lookup_es3_des3_1d
1893
1894 !#######################################################################
1895
1896 ! <SUBROUTINE NAME="lookup_es3_des3_2d" INTERFACE="lookup_es3_des3">
1897 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
1898 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
1899 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:)"></OUT>
1900 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1901 ! </SUBROUTINE>
1902  subroutine lookup_es3_des3_2d ( temp, esat, desat, err_msg )
1903
1904  real, dimension(:,:), intent(in)  :: temp
1905  real, dimension(:,:), intent(out) :: esat, desat
1906  character(len=*), intent(out), optional :: err_msg
1907
1908  integer :: nbad
1909  character(len=128) :: err_msg_local
1910
1911    if (.not.module_is_initialized) then
1912       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1913    endif
1914
1915    call lookup_es3_des3_k(temp, esat, desat, nbad)
1916
1917    if ( nbad == 0 ) then
1918      if(present(err_msg)) err_msg = ''
1919    else
1920      if(show_bad_value_count_by_slice) call temp_check ( temp )
1921      if(show_all_bad_values) call show_all_bad ( temp )
1922      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1923      if(fms_error_handler('lookup_es3_des3',err_msg_local,err_msg)) return
1924    endif
1925
1926  end subroutine lookup_es3_des3_2d
1927
1928 !#######################################################################
1929
1930 ! <SUBROUTINE NAME="lookup_es3_des3_3d" INTERFACE="lookup_es3_des3">
1931 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
1932 !   <OUT NAME="esat" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
1933 !   <OUT NAME="desat" UNITS="pascal / degree" TYPE="real" DIM="(:,:,:)"></OUT>
1934 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1935 ! </SUBROUTINE>
1936  subroutine lookup_es3_des3_3d ( temp, esat, desat, err_msg )
1937
1938  real, dimension(:,:,:), intent(in)  :: temp
1939  real, dimension(:,:,:), intent(out) :: esat, desat
1940  character(len=*), intent(out), optional :: err_msg
1941
1942  integer :: nbad
1943  character(len=128) :: err_msg_local
1944
1945    if (.not.module_is_initialized) then
1946       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1947    endif
1948
1949    call lookup_es3_des3_k(temp, esat, desat, nbad)
1950
1951    if ( nbad == 0 ) then
1952      if(present(err_msg)) err_msg = ''
1953    else
1954      if(show_bad_value_count_by_slice) call temp_check ( temp )
1955      if(show_all_bad_values) call show_all_bad ( temp )
1956      write(err_msg_local,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
1957      if(fms_error_handler('lookup_es3_des3',err_msg_local,err_msg)) return
1958    endif
1959
1960  end subroutine lookup_es3_des3_3d
1961
1962 !#######################################################################
1963
1964 ! <SUBROUTINE NAME="compute_qs_0d" INTERFACE="compute_qs">
1965 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(SCALAR)"></IN>
1966 !   <IN NAME="press UNIT="Pascals" TYPE="real" DIM="(SCALAR)"></IN>
1967 !   <OUT NAME="qsat" UNITS="kg(vapor)/kg(moist air)" TYPE="real" DIM="(SCALAR)"></OUT>
1968 !   <IN NAME="q" UNIT="kg(vapor)/kg(moistair)" TYPE="real" DIM="(SCALAR)"></IN>
1969 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
1970 !   <OUT NAME="dqsdT" UNIT="kg(vapor)/kg(moistair)/ degree Kelvin" TYPE="real" DIM="(SCALAR)"></OUT>
1971 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar)"> </OUT>
1972 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
1973 ! </SUBROUTINE>
1974  subroutine compute_qs_0d ( temp, press, qsat, q, hc, dqsdT, esat, &
1975                             err_msg, es_over_liq, es_over_liq_and_ice )
1976
1977  class(*), intent(in)                    :: temp, press
1978  class(*), intent(out)                   :: qsat
1979  class(*), intent(in),          optional :: q, hc
1980  class(*), intent(out),         optional :: dqsdT, esat
1981  character(len=*), intent(out), optional :: err_msg
1982  logical,intent(in),            optional :: es_over_liq
1983  logical,intent(in),            optional :: es_over_liq_and_ice
1984
1985  integer :: nbad
1986  character(len=128) :: err_msg_tmp
1987
1988    if (.not.module_is_initialized) then
1989       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
1990    endif
1991
1992    if (present(es_over_liq)) then
1993      if (.not. (construct_table_wrt_liq)) then
1994        call error_mesg ('compute_qs', &
1995           'requesting es wrt liq, but that table not constructed', &
1996                                                                 FATAL)
1997      endif
1998    endif
1999    if (present(es_over_liq_and_ice)) then
2000      if (.not. (construct_table_wrt_liq_and_ice)) then
2001        call error_mesg ('compute_qs', &
2002       'requesting es wrt liq and ice, but that table not constructed', &
2003                                                                 FATAL)
2004      endif
2005    endif
2006
2007    call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, hc, &
2008                        dqsdT, esat, es_over_liq, es_over_liq_and_ice)
2009
2010    if ( nbad == 0 ) then
2011      if(present(err_msg)) err_msg = ''
2012    else
2013      if(show_all_bad_values) call show_all_bad ( temp )
2014      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2015      if(fms_error_handler('compute_qs',err_msg_tmp,err_msg)) return
2016    endif
2017
2018  end subroutine compute_qs_0d
2019
2020 !#######################################################################
2021
2022 ! <SUBROUTINE NAME="compute_qs_1d" INTERFACE="compute_qs">
2023 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
2024 !   <IN NAME="press UNIT="Pascals" TYPE="real" DIM="(:)"></IN>
2025 !   <OUT NAME="qsat" UNITS="kg(vapor)/kg(moist air)" TYPE="real" DIM="(:)"></OUT>
2026 !   <IN NAME="q" UNIT="kg(vapor)/kg(moistair)" TYPE="real" DIM="(:)"></IN>
2027 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2028 !   <OUT NAME="dqsdT" UNIT="kg(vapor)/kg(moistair)/ degree Kelvin" TYPE="real" DIM="(:)"></OUT>
2029 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:)"> </OUT>
2030 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2031 ! </SUBROUTINE>
2032  subroutine compute_qs_1d ( temp, press, qsat, q, hc, dqsdT, esat, &
2033                             err_msg, es_over_liq, es_over_liq_and_ice )
2034
2035  class(*), intent(in)                    :: temp(:), press(:)
2036  class(*), intent(out)                   :: qsat(:)
2037  class(*), intent(in),          optional :: q(:)
2038  class(*), intent(in),          optional :: hc
2039  class(*), intent(out),         optional :: dqsdT(:), esat(:)
2040  character(len=*), intent(out), optional :: err_msg
2041  logical,intent(in),            optional :: es_over_liq
2042  logical,intent(in),            optional :: es_over_liq_and_ice
2043
2044  integer :: nbad
2045  character(len=128) :: err_msg_tmp
2046
2047    if (.not.module_is_initialized) then
2048       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2049    endif
2050
2051    if (present(es_over_liq)) then
2052      if (.not. (construct_table_wrt_liq)) then
2053        call error_mesg ('compute_qs', &
2054           'requesting es wrt liq, but that table not constructed', &
2055                                                                 FATAL)
2056      endif
2057    endif
2058    if (present(es_over_liq_and_ice)) then
2059      if (.not. (construct_table_wrt_liq_and_ice)) then
2060        call error_mesg ('compute_qs', &
2061       'requesting es wrt liq and ice, but that table not constructed', &
2062                                                                 FATAL)
2063      endif
2064    endif
2065
2066 !  call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, dqsdT)
2067    call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, hc, &
2068                        dqsdT, esat, es_over_liq, es_over_liq_and_ice)
2069
2070    if ( nbad == 0 ) then
2071      if(present(err_msg)) err_msg = ''
2072    else
2073      if(show_bad_value_count_by_slice) call temp_check ( temp )
2074      if(show_all_bad_values) call show_all_bad ( temp )
2075      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2076      if(fms_error_handler('compute_qs',err_msg_tmp,err_msg)) return
2077    endif
2078
2079  end subroutine compute_qs_1d
2080
2081
2082 !#######################################################################
2083
2084 ! <SUBROUTINE NAME="compute_qs_2d" INTERFACE="compute_qs">
2085 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
2086 !   <IN NAME="press UNIT="Pascals" TYPE="real" DIM="(:,:)"></IN>
2087 !   <OUT NAME="qsat" UNITS="kg(vapor)/kg(moist air)" TYPE="real" DIM="(;,:)"></OUT>
2088 !   <IN NAME="q" UNIT="kg(vapor)/kg(moistair)" TYPE="real" DIM="(:,:)"></IN>
2089 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2090 !   <OUT NAME="dqsdT" UNIT="kg(vapor)/kg(moistair)/ degree Kelvin" TYPE="real" DIM="(:,:)"></OUT>
2091 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:,:)"> </OUT>
2092 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2093 ! </SUBROUTINE>
2094  subroutine compute_qs_2d ( temp, press, qsat, q, hc, dqsdT, esat, &
2095                             err_msg, es_over_liq, es_over_liq_and_ice )
2096
2097  class(*), intent(in)                    :: temp(:,:), press(:,:)
2098  class(*), intent(out)                   :: qsat(:,:)
2099  class(*), intent(in),          optional :: q(:,:)
2100  class(*), intent(in),          optional :: hc
2101  class(*), intent(out),         optional :: dqsdT(:,:), esat(:,:)
2102  character(len=*), intent(out), optional :: err_msg
2103  logical,intent(in),            optional :: es_over_liq
2104  logical,intent(in),            optional :: es_over_liq_and_ice
2105
2106  integer :: nbad
2107  character(len=128) :: err_msg_tmp
2108
2109    if (.not.module_is_initialized) then
2110       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2111    endif
2112
2113    if (present(es_over_liq)) then
2114      if (.not. (construct_table_wrt_liq)) then
2115        call error_mesg ('compute_qs', &
2116           'requesting es wrt liq, but that table not constructed', &
2117                                                                 FATAL)
2118      endif
2119    endif
2120    if (present(es_over_liq_and_ice)) then
2121      if (.not. (construct_table_wrt_liq_and_ice)) then
2122        call error_mesg ('compute_qs', &
2123       'requesting es wrt liq and ice, but that table not constructed', &
2124                                                                 FATAL)
2125      endif
2126    endif
2127
2128 !  call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, dqsdT)
2129    call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, hc, &
2130                        dqsdT, esat, es_over_liq, es_over_liq_and_ice)
2131
2132    if ( nbad == 0 ) then
2133      if(present(err_msg)) err_msg = ''
2134    else
2135      if(show_bad_value_count_by_slice) call temp_check ( temp )
2136      if(show_all_bad_values) call show_all_bad ( temp )
2137      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2138      if(fms_error_handler('compute_qs',err_msg_tmp,err_msg)) return
2139    endif
2140
2141  end subroutine compute_qs_2d
2142
2143 !#######################################################################
2144
2145 ! <SUBROUTINE NAME="compute_qs_3d" INTERFACE="compute_qs">
2146 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
2147 !   <IN NAME="press UNIT="Pascals" TYPE="real" DIM="(:,:,:)"></IN>
2148 !   <OUT NAME="qsat" UNITS="kg(vapor)/kg(moist air)" TYPE="real" DIM="(;,:,:)"></OUT>
2149 !   <IN NAME="q" UNIT="kg(vapor)/kg(moistair)" TYPE="real" DIM="(:,:,:)"></IN>
2150 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2151 !   <OUT NAME="dqsdT" UNIT="kg(vapor)/kg(moistair)/ degree Kelvin" TYPE="real" DIM="(:,:,:)"></OUT>
2152 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:,:,:)"> </OUT>
2153 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2154 ! </SUBROUTINE>
2155  subroutine compute_qs_3d ( temp, press, qsat, q, hc, dqsdT, esat, &
2156                             err_msg, es_over_liq, es_over_liq_and_ice )
2157
2158  class(*), intent(in)                    :: temp(:,:,:), press(:,:,:)
2159  class(*), intent(out)                   :: qsat(:,:,:)
2160  class(*), intent(in),          optional :: q(:,:,:)
2161  class(*), intent(in),          optional :: hc
2162  class(*), intent(out),         optional :: dqsdT(:,:,:), esat(:,:,:)
2163  character(len=*), intent(out), optional :: err_msg
2164  logical,intent(in),            optional :: es_over_liq
2165  logical,intent(in),            optional :: es_over_liq_and_ice
2166
2167  integer :: nbad
2168  character(len=128) :: err_msg_tmp
2169
2170    if (.not.module_is_initialized) then
2171       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2172    endif
2173
2174    if (present(es_over_liq)) then
2175      if (.not. (construct_table_wrt_liq)) then
2176        call error_mesg ('compute_qs', &
2177           'requesting es wrt liq, but that table not constructed', &
2178                                                                 FATAL)
2179      endif
2180    endif
2181    if (present(es_over_liq_and_ice)) then
2182      if (.not. (construct_table_wrt_liq_and_ice)) then
2183        call error_mesg ('compute_qs', &
2184       'requesting es wrt liq and ice, but that table not constructed', &
2185                                                                 FATAL)
2186      endif
2187    endif
2188
2189 !  call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, dqsdT)
2190    call compute_qs_k (temp, press,  EPSILO, ZVIR, qsat, nbad, q, hc, &
2191                        dqsdT, esat, es_over_liq, es_over_liq_and_ice)
2192
2193
2194    if ( nbad == 0 ) then
2195      if(present(err_msg)) err_msg = ''
2196    else
2197      if(show_bad_value_count_by_slice) call temp_check ( temp )
2198      if(show_all_bad_values) call show_all_bad ( temp )
2199      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2200      if(fms_error_handler('compute_qs',err_msg_tmp,err_msg)) return
2201    endif
2202
2203  end subroutine compute_qs_3d
2204
2205 !#######################################################################
2206 !#######################################################################
2207
2208 ! <SUBROUTINE NAME="compute_mrs_0d" INTERFACE="compute_mrs">
2209 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(SCALAR)"></IN>
2210 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(SCALAR)"></IN>
2211 !   <OUT NAME="mrsat" UNITS="kg(vapor)/kg(dry air)" TYPE="real" DIM="(SCALAR</OUT>
2212 !   <IN NAME="mr" UNIT="kg(vapor)/kg(dry air)" TYPE="real" DIM="(SCALAR)"></IN>
2213 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2214 !   <OUT NAME="dmrsdT" UNIT="kg(vapor)/kg(dry air)/ degree Kelvin" TYPE="real" DIM="(SCALAR)"></OUT>
2215 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(scalar)"> </OUT>
2216 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2217 ! </SUBROUTINE>
2218  subroutine compute_mrs_0d ( temp, press, mrsat, mr, hc, dmrsdT, esat, &
2219                             err_msg, es_over_liq, es_over_liq_and_ice )
2220
2221  real, intent(in)                        :: temp, press
2222  real, intent(out)                       :: mrsat
2223  real, intent(in),              optional :: mr, hc
2224  real, intent(out),             optional :: dmrsdT, esat
2225  character(len=*), intent(out), optional :: err_msg
2226  logical,intent(in),            optional :: es_over_liq
2227  logical,intent(in),            optional :: es_over_liq_and_ice
2228
2229  integer :: nbad
2230  character(len=128) :: err_msg_tmp
2231
2232    if (.not.module_is_initialized) then
2233       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2234    endif
2235
2236    if (present(es_over_liq)) then
2237      if (.not. (construct_table_wrt_liq)) then
2238        call error_mesg ('compute_mrs', &
2239           'requesting es wrt liq, but that table not constructed', &
2240                                                                 FATAL)
2241      endif
2242    endif
2243    if (present(es_over_liq_and_ice)) then
2244      if (.not. (construct_table_wrt_liq_and_ice)) then
2245        call error_mesg ('compute_qs', &
2246       'requesting es wrt liq and ice, but that table not constructed', &
2247                                                                 FATAL)
2248      endif
2249    endif
2250
2251    call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat, nbad, mr,  &
2252                      hc, dmrsdT, esat, es_over_liq, es_over_liq_and_ice)
2253
2254    if ( nbad == 0 ) then
2255      if(present(err_msg)) err_msg = ''
2256    else
2257      if(show_all_bad_values) call show_all_bad ( temp )
2258      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2259      if(fms_error_handler('compute_mrs',err_msg_tmp,err_msg)) return
2260    endif
2261
2262  end subroutine compute_mrs_0d
2263
2264 !#######################################################################
2265 !#######################################################################
2266
2267 ! <SUBROUTINE NAME="compute_mrs_1d" INTERFACE="compute_mrs">
2268 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
2269 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(:)"></IN>
2270 !   <OUT NAME="mrsat" UNITS="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:)"></OUT>
2271 !   <IN NAME="mr" UNIT="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:)"></IN>
2272 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2273 !   <OUT NAME="dmrsdT" UNIT="kg(vapor)/kg(dry air)/ degree Kelvin" TYPE="real" DIM="(:)"></OUT>
2274 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:)"> </OUT>
2275 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2276 ! </SUBROUTINE>
2277  subroutine compute_mrs_1d ( temp, press, mrsat, mr, hc, dmrsdT, esat,&
2278                             err_msg, es_over_liq, es_over_liq_and_ice )
2279
2280  real, intent(in)                        :: temp(:), press(:)
2281  real, intent(out)                       :: mrsat(:)
2282  real, intent(in),              optional :: mr(:)
2283  real, intent(in),              optional :: hc
2284  real, intent(out),             optional :: dmrsdT(:), esat(:)
2285  character(len=*), intent(out), optional :: err_msg
2286  logical,intent(in),            optional :: es_over_liq
2287  logical,intent(in),            optional :: es_over_liq_and_ice
2288
2289  integer :: nbad
2290  character(len=128) :: err_msg_tmp
2291
2292    if (.not.module_is_initialized) then
2293       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2294    endif
2295
2296    if (present(es_over_liq)) then
2297      if (.not. (construct_table_wrt_liq)) then
2298        call error_mesg ('compute_mrs', &
2299           'requesting es wrt liq, but that table not constructed', &
2300                                                                 FATAL)
2301      endif
2302    endif
2303    if (present(es_over_liq_and_ice)) then
2304      if (.not. (construct_table_wrt_liq_and_ice)) then
2305        call error_mesg ('compute_qs', &
2306       'requesting es wrt liq and ice, but that table not constructed', &
2307                                                                 FATAL)
2308      endif
2309    endif
2310
2311 !  call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat,  &
2312 !                                                     nbad, mr, dmrsdT)
2313    call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat, nbad, mr,  &
2314                      hc, dmrsdT, esat, es_over_liq, es_over_liq_and_ice)
2315
2316    if ( nbad == 0 ) then
2317      if(present(err_msg)) err_msg = ''
2318    else
2319      if(show_bad_value_count_by_slice) call temp_check ( temp )
2320      if(show_all_bad_values) call show_all_bad ( temp )
2321      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2322      if(fms_error_handler('compute_mrs',err_msg_tmp,err_msg)) return
2323    endif
2324
2325  end subroutine compute_mrs_1d
2326
2327 !#######################################################################
2328
2329 ! <SUBROUTINE NAME="compute_mrs_2d" INTERFACE="compute_mrs">
2330 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
2331 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(:,:)"></IN>
2332 !   <OUT NAME="mrsat" UNITS="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:,:)"></OUT>
2333 !   <IN NAME="mr" UNIT="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:,:)"></IN>
2334 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2335 !   <OUT NAME="dmrsdT" UNIT="kg(vapor)/kg(dry air)/ degree Kelvin" TYPE="real" DIM="(:,:)"></OUT>
2336 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:,:)"> </OUT>
2337 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2338 ! </SUBROUTINE>
2339  subroutine compute_mrs_2d ( temp, press, mrsat, mr, hc, dmrsdT, esat,&
2340                             err_msg, es_over_liq, es_over_liq_and_ice )
2341
2342  real, intent(in)                        :: temp(:,:), press(:,:)
2343  real, intent(out)                       :: mrsat(:,:)
2344  real, intent(in),              optional :: mr(:,:)
2345  real, intent(in),              optional :: hc
2346  real, intent(out),             optional :: dmrsdT(:,:), esat(:,:)
2347  character(len=*), intent(out), optional :: err_msg
2348  logical,intent(in),            optional :: es_over_liq
2349  logical,intent(in),            optional :: es_over_liq_and_ice
2350
2351  integer :: nbad
2352  character(len=128) :: err_msg_tmp
2353
2354    if (.not.module_is_initialized) then
2355       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2356    endif
2357
2358    if (present(es_over_liq)) then
2359      if (.not. (construct_table_wrt_liq)) then
2360        call error_mesg ('compute_mrs', &
2361           'requesting es wrt liq, but that table not constructed', &
2362                                                                 FATAL)
2363      endif
2364    endif
2365    if (present(es_over_liq_and_ice)) then
2366      if (.not. (construct_table_wrt_liq_and_ice)) then
2367        call error_mesg ('compute_qs', &
2368       'requesting es wrt liq and ice, but that table not constructed', &
2369                                                                 FATAL)
2370      endif
2371    endif
2372
2373 !  call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat,  &
2374 !                                                     nbad, mr, dmrsdT)
2375    call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat, nbad, mr,  &
2376                      hc, dmrsdT, esat, es_over_liq, es_over_liq_and_ice)
2377
2378    if ( nbad == 0 ) then
2379      if(present(err_msg)) err_msg = ''
2380    else
2381      if(show_bad_value_count_by_slice) call temp_check ( temp )
2382      if(show_all_bad_values) call show_all_bad ( temp )
2383      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2384      if(fms_error_handler('compute_mrs',err_msg_tmp,err_msg)) return
2385    endif
2386
2387  end subroutine compute_mrs_2d
2388
2389 !#######################################################################
2390
2391 ! <SUBROUTINE NAME="compute_mrs_3d" INTERFACE="compute_mrs">
2392 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
2393 !   <IN NAME="press" UNIT="Pascals" TYPE="real" DIM="(:,:,:)"></IN>
2394 !   <OUT NAME="mrsat" UNITS="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:,:,:)"></OUT>
2395 !   <IN NAME="mr" UNIT="kg(vapor)/kg(dry air)" TYPE="real" DIM="(:,:,:)"></IN>
2396 !   <IN NAME="hc" UNIT="fraction" TYPE="real" DIM="(scalar)"></IN>
2397 !   <OUT NAME="dmrsdT" UNIT="kg(vapor)/kg(dry air)/ degree Kelvin" TYPE="real" DIM="(:,:,:)"></OUT>
2398 !   <OUT NAME="esat" UNITS="Pascals" TYPE="real" DIM="(:,:,:)"> </OUT>
2399 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2400 ! </SUBROUTINE>
2401  subroutine compute_mrs_3d ( temp, press, mrsat, mr, hc, dmrsdT, esat,&
2402                             err_msg, es_over_liq, es_over_liq_and_ice )
2403
2404  real, intent(in)                        :: temp(:,:,:), press(:,:,:)
2405  real, intent(out)                       :: mrsat(:,:,:)
2406  real, intent(in),              optional :: mr(:,:,:)
2407  real, intent(in),              optional :: hc
2408  real, intent(out),             optional :: dmrsdT(:,:,:), esat(:,:,:)
2409  character(len=*), intent(out), optional :: err_msg
2410  logical,intent(in),            optional :: es_over_liq
2411  logical,intent(in),            optional :: es_over_liq_and_ice
2412
2413  integer :: nbad
2414  character(len=128) :: err_msg_tmp
2415
2416    if (.not.module_is_initialized) then
2417       if(fms_error_handler('lookup_es','sat_vapor_pres_init is not called' ,err_msg)) return
2418    endif
2419
2420    if (present(es_over_liq)) then
2421      if (.not. (construct_table_wrt_liq)) then
2422        call error_mesg ('compute_mrs', &
2423           'requesting es wrt liq, but that table not constructed', &
2424                                                                 FATAL)
2425      endif
2426    endif
2427    if (present(es_over_liq_and_ice)) then
2428      if (.not. (construct_table_wrt_liq_and_ice)) then
2429        call error_mesg ('compute_qs', &
2430       'requesting es wrt liq and ice, but that table not constructed', &
2431                                                                 FATAL)
2432      endif
2433    endif
2434
2435 !  call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat,   &
2436 !                                                    nbad, mr, dmrsdT)
2437    call compute_mrs_k (temp, press, EPSILO, ZVIR, mrsat, nbad, mr,  &
2438                      hc, dmrsdT, esat, es_over_liq, es_over_liq_and_ice)
2439
2440    if ( nbad == 0 ) then
2441      if(present(err_msg)) err_msg = ''
2442    else
2443      if(show_bad_value_count_by_slice) call temp_check ( temp )
2444      if(show_all_bad_values) call show_all_bad ( temp )
2445      write(err_msg_tmp,'(a47,i7)') 'saturation vapor pressure table overflow, nbad=', nbad
2446      if(fms_error_handler('compute_mrs',err_msg_tmp,err_msg)) return
2447    endif
2448
2449  end subroutine compute_mrs_3d
2450
2451
2452 !#######################################################################
2453
2454 !#######################################################################
2455
2456 ! <SUBROUTINE NAME="sat_vapor_pres_init">
2457
2458 !   <OVERVIEW>
2459 !     Initializes the lookup tables for saturation vapor pressure.
2460 !   </OVERVIEW>
2461 !   <DESCRIPTION>
2462 !     Initializes the lookup tables for saturation vapor pressure.
2463 !     This routine will be called automatically the first time
2464 !     <B>lookup_es</B> or <B>lookup_des</B> is called,
2465 !     the user does not need to call this routine.
2466 !     There are no arguments.
2467 !   </DESCRIPTION>
2468 !   <TEMPLATE>
2469 !     call sat_vapor_pres_init
2470 !   </TEMPLATE>
2471 !   <OUT NAME="err_msg" TYPE="character">  </OUT>
2472
2473 ! </SUBROUTINE>
2474  subroutine sat_vapor_pres_init(err_msg)
2475
2476 !  =================================================================
2477 !  +                                                               +
2478 !  +             construction of the es table                      +
2479 !  +                                                               +
2480 !  + this table is constructed from es equations from the          +
2481 !  + smithsonian tables.  the es input is computed from values     +
2482 !  + (in one-tenth of a degree increments) of es over ice          +
2483 !  + from -153c to 0c and values of es over water from 0c to 102c. +
2484 !  + output table contains these data interleaved with their       +
2485 !  + derivatives with respect to temperature except between -20c   +
2486 !  + and 0c where blended (over water and over ice) es values and  +
2487 !  + derivatives are calculated.                                   +
2488 !  +   note: all es computation is done in pascals                 +
2489 !  =================================================================
2490
2491   character(len=*), intent(out), optional :: err_msg
2492   character(len=128) :: err_msg_local
2493   integer :: unit, ierr, io
2494
2495 ! return silently if this routine has already been called
2496   if (module_is_initialized) return
2497
2498 !---- read namelist input ----
2499   read (input_nml_file, sat_vapor_pres_nml, iostat=io)
2500   ierr = check_nml_error(io,'sat_vapor_pres_nml')
2501
2502 ! write version number and namelist to log file
2503   call write_version_number("SAT_VAPOR_PRES_MOD", version)
2504   unit = stdlog()
2505   stdoutunit = stdout()
2506   if (mpp_pe() == mpp_root_pe()) write (unit, nml=sat_vapor_pres_nml)
2507
2508   if(do_simple) then
2509     tcmin = -173
2510     tcmax =  350
2511   endif
2512   nsize = (tcmax-tcmin)*esres+1
2513   nlim  = nsize-1
2514   call sat_vapor_pres_init_k(nsize, real(tcmin), real(tcmax), TFREEZE, HLV, &
2515                              RVGAS, ES0, err_msg_local, use_exact_qs, do_simple, &
2516                              construct_table_wrt_liq, &
2517                              construct_table_wrt_liq_and_ice, &
2518                              teps, tmin, dtinv)
2519   if ( err_msg_local == '' ) then
2520      if(present(err_msg)) err_msg = ''
2521   else
2522      if(fms_error_handler('lookup_es',err_msg_local,err_msg)) return
2523   endif
2524
2525   module_is_initialized = .true.
2526
2527 end subroutine sat_vapor_pres_init
2528
2529 !#######################################################################
2530 !#######################################################################
2531 !-------------------------------------------------------------------
2532 !                Computation of the es values
2533 !
2534 !   Saturation vapor pressure (es) values are computed from
2535 !   equations in the Smithsonian meteorological tables page 350.
2536 !   For temperatures < 0C, sat vapor pres is computed over ice.
2537 !   For temperatures > -20C, sat vapor pres is computed over water.
2538 !   Between -20C and 0C the returned value is blended (over water
2539 !   and over ice).  All sat vapor pres values are returned in pascals.
2540 !
2541 !   Reference:  Smithsonian meteorological tables, page 350.
2542 !-------------------------------------------------------------------
2543
2544 ! <FUNCTION NAME="compute_es_1d" INTERFACE="compute_es">
2545 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:)"></IN>
2546 !   <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:)"></OUT>
2547 ! </FUNCTION>
2548 !function compute_es_1d (tem) result (es)
2549 !real, intent(in) :: tem(:)
2550 !real :: es(size(tem,1))
2551
2552 !es = compute_es_k(tem, TFREEZE)
2553
2554 !end function compute_es_1d
2555 !--------------------------------------------------------
2556
2557 ! <FUNCTION NAME="compute_es_0d" INTERFACE="compute_es">
2558 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(scalar)"></IN>
2559 !   <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(scalar)"></OUT>
2560 ! </FUNCTION>
2561 !function compute_es_0d (tem) result (es)
2562 !real, intent(in) :: tem
2563 !real :: es
2564 !real, dimension(1) :: tem1, es1
2565
2566 !  tem1(1) = tem
2567 !  es1 = compute_es_1d (tem1)
2568 !  es = es1(1)
2569
2570 !end function compute_es_0d
2571
2572 !--------------------------
2573
2574 ! <FUNCTION NAME="compute_es_2d" INTERFACE="compute_es">
2575 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:)"></IN>
2576 !   <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:,:)"></OUT>
2577 ! </FUNCTION>
2578 !function compute_es_2d (tem) result (es)
2579 !real, intent(in) :: tem(:,:)
2580 !real, dimension(size(tem,1),size(tem,2)) :: es
2581 !integer :: j
2582
2583 !   do j = 1, size(tem,2)
2584 !     es(:,j) = compute_es_1d (tem(:,j))
2585 !   enddo
2586
2587 !end function compute_es_2d
2588
2589 !--------------------------
2590 ! <FUNCTION NAME="compute_es_3d" INTERFACE="compute_es">
2591 !   <IN NAME="temp" UNIT="degrees Kelvin" TYPE="real" DIM="(:,:,:)"></IN>
2592 !   <OUT NAME="es" UNITS="pascal" TYPE="real" DIM="(:,:,:)"></OUT>
2593 ! </FUNCTION>
2594 !function compute_es_3d (tem) result (es)
2595 !real, intent(in) :: tem(:,:,:)
2596 !real, dimension(size(tem,1),size(tem,2),size(tem,3)) :: es
2597 !integer :: j, k
2598
2599 !   do k = 1, size(tem,3)
2600 !   do j = 1, size(tem,2)
2601 !     es(:,j,k) = compute_es_1d (tem(:,j,k))
2602 !   enddo
2603 !   enddo
2604
2605 !end function compute_es_3d
2606
2607 !#######################################################################
2608
2609  function check_1d ( temp ) result ( nbad )
2610  class(*), intent(in)  :: temp(:)
2611  integer :: nbad, ind, i
2612
2613    nbad = 0
2614
2615    select type (temp)
2616    type is (real(kind=r4_kind))
2617      do i = 1, size(temp,1)
2618        ind = int(dtinv*(temp(i)-tmin+teps))
2619        if (ind < 0 .or. ind > nlim) nbad = nbad+1
2620      enddo
2621    type is (real(kind=r8_kind))
2622      do i = 1, size(temp,1)
2623        ind = int(dtinv*(temp(i)-tmin+teps))
2624        if (ind < 0 .or. ind > nlim) nbad = nbad+1
2625      enddo
2626    class default
2627      call error_mesg ('sat_vapor_pres_mod::check_1d',&
2628           & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2629    end select
2630
2631  end function check_1d
2632
2633 !------------------------------------------------
2634
2635  function check_2d ( temp ) result ( nbad )
2636  class(*), intent(in)  :: temp(:,:)
2637  integer :: nbad
2638  integer :: j
2639
2640    nbad = 0
2641
2642    select type (temp)
2643    type is (real(kind=r4_kind))
2644      do j = 1, size(temp,2)
2645        nbad = nbad + check_1d ( temp(:,j) )
2646      enddo
2647    type is (real(kind=r8_kind))
2648      do j = 1, size(temp,2)
2649        nbad = nbad + check_1d ( temp(:,j) )
2650      enddo
2651    class default
2652      call error_mesg ('sat_vapor_pres_mod::check_2d',&
2653           &  'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2654    end select
2655
2656  end function check_2d
2657
2658 !#######################################################################
2659
2660  subroutine temp_check_1d ( temp )
2661  class(*), intent(in) :: temp(:)
2662  integer :: i, unit
2663
2664    unit = stdoutunit
2665
2666    select type (temp)
2667    type is (real(kind=r4_kind))
2668      write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i:i)),i=1,size(temp,1))
2669    type is (real(kind=r8_kind))
2670      write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i:i)),i=1,size(temp,1))
2671    class default
2672      call error_mesg ('sat_vapor_pres_mod::temp_check_1d',&
2673           & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2674    end select
2675
2676  end subroutine temp_check_1d
2677
2678 !--------------------------------------------------------------
2679
2680  subroutine temp_check_2d ( temp )
2681  class(*), intent(in) :: temp(:,:)
2682  integer :: i, j, unit
2683
2684    unit = stdoutunit
2685
2686    select type (temp)
2687    type is (real(kind=r4_kind))
2688      write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i,:)),i=1,size(temp,1))
2689      write(unit,*) 'Bad temperatures (dimension 2): ', (check_1d(temp(:,j)),j=1,size(temp,2))
2690    type is (real(kind=r8_kind))
2691      write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i,:)),i=1,size(temp,1))
2692      write(unit,*) 'Bad temperatures (dimension 2): ', (check_1d(temp(:,j)),j=1,size(temp,2))
2693    class default
2694      call error_mesg ('sat_vapor_pres_mod::temp_check_2d',&
2695           & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2696    end select
2697
2698  end subroutine temp_check_2d
2699
2700 !--------------------------------------------------------------
2701
2702  subroutine temp_check_3d ( temp )
2703  class(*), intent(in)  :: temp(:,:,:)
2704  integer :: i, j, k, unit
2705
2706    unit = stdoutunit
2707
2708    select type (temp)
2709    type is (real(kind=r4_kind))
2710      write(unit,*) 'Bad temperatures (dimension 1): ', (check_2d(temp(i,:,:)),i=1,size(temp,1))
2711      write(unit,*) 'Bad temperatures (dimension 2): ', (check_2d(temp(:,j,:)),j=1,size(temp,2))
2712      write(unit,*) 'Bad temperatures (dimension 3): ', (check_2d(temp(:,:,k)),k=1,size(temp,3))
2713    type is (real(kind=r8_kind))
2714      write(unit,*) 'Bad temperatures (dimension 1): ', (check_2d(temp(i,:,:)),i=1,size(temp,1))
2715      write(unit,*) 'Bad temperatures (dimension 2): ', (check_2d(temp(:,j,:)),j=1,size(temp,2))
2716      write(unit,*) 'Bad temperatures (dimension 3): ', (check_2d(temp(:,:,k)),k=1,size(temp,3))
2717    class default
2718      call error_mesg ('sat_vapor_pres_mod::temp_check_3d',&
2719           & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2720    end select
2721
2722  end subroutine temp_check_3d
2723
2724 !#######################################################################
2725
2726 subroutine show_all_bad_0d ( temp )
2727  class(*), intent(in) :: temp
2728  integer :: ind, unit
2729
2730  unit = stdoutunit
2731
2732  select type (temp)
2733  type is (real(kind=r4_kind))
2734    ind = int(dtinv*(temp-tmin+teps))
2735    if (ind < 0 .or. ind > nlim) then
2736      write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe()
2737    endif
2738  type is (real(kind=r8_kind))
2739    ind = int(dtinv*(temp-tmin+teps))
2740    if (ind < 0 .or. ind > nlim) then
2741      write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe()
2742    endif
2743  class default
2744    call error_mesg ('sat_vapor_pres_mod::show_all_bad_0d',&
2745         & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2746  end select
2747
2748  end subroutine show_all_bad_0d
2749
2750 !--------------------------------------------------------------
2751
2752  subroutine show_all_bad_1d ( temp )
2753  class(*), intent(in) :: temp(:)
2754  integer :: i, ind, unit
2755
2756  unit = stdoutunit
2757
2758  select type (temp)
2759  type is (real(kind=r4_kind))
2760    do i=1,size(temp)
2761      ind = int(dtinv*(temp(i)-tmin+teps))
2762      if (ind < 0 .or. ind > nlim) then
2763        write(unit,'(a,e10.3,a,i4,a,i6)') 'Bad temperature=',temp(i),'  at i=',i,' pe=',mpp_pe()
2764      endif
2765    enddo
2766  type is (real(kind=r8_kind))
2767    do i=1,size(temp)
2768      ind = int(dtinv*(temp(i)-tmin+teps))
2769      if (ind < 0 .or. ind > nlim) then
2770        write(unit,'(a,e10.3,a,i4,a,i6)') 'Bad temperature=',temp(i),'  at i=',i,' pe=',mpp_pe()
2771      endif
2772    enddo
2773  class default
2774    call error_mesg ('sat_vapor_pres_mod::show_all_bad_1d',&
2775         & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2776  end select
2777
2778  end subroutine show_all_bad_1d
2779
2780 !--------------------------------------------------------------
2781
2782  subroutine show_all_bad_2d ( temp )
2783  class(*), intent(in) :: temp(:,:)
2784  integer :: i, j, ind, unit
2785
2786  unit = stdoutunit
2787
2788  select type (temp)
2789  type is (real(kind=r4_kind))
2790    do j=1,size(temp,2)
2791    do i=1,size(temp,1)
2792      ind = int(dtinv*(temp(i,j)-tmin+teps))
2793      if (ind < 0 .or. ind > nlim) then
2794        write(unit,'(a,e10.3,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j),'  at i=',i,' j=',j,' pe=',mpp_pe()
2795      endif
2796    enddo
2797    enddo
2798  type is (real(kind=r8_kind))
2799    do j=1,size(temp,2)
2800    do i=1,size(temp,1)
2801      ind = int(dtinv*(temp(i,j)-tmin+teps))
2802      if (ind < 0 .or. ind > nlim) then
2803        write(unit,'(a,e10.3,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j),'  at i=',i,' j=',j,' pe=',mpp_pe()
2804      endif
2805    enddo
2806    enddo
2807  class default
2808    call error_mesg ('sat_vapor_pres_mod::show_all_bad_2d',&
2809         & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2810  end select
2811
2812  end subroutine show_all_bad_2d
2813
2814 !--------------------------------------------------------------
2815
2816  subroutine show_all_bad_3d ( temp )
2817  class(*), intent(in)  :: temp(:,:,:)
2818  integer :: i, j, k, ind, unit
2819
2820  unit = stdoutunit
2821
2822  select type (temp)
2823  type is (real(kind=r4_kind))
2824    do k=1,size(temp,3)
2825    do j=1,size(temp,2)
2826    do i=1,size(temp,1)
2827      ind = int(dtinv*(temp(i,j,k)-tmin+teps))
2828      if (ind < 0 .or. ind > nlim) then
2829        write(unit,'(a,e10.3,a,i4,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j,k),&
2830             &'  at i=',i,' j=',j,' k=',k,' pe=',mpp_pe()
2831      endif
2832    enddo
2833    enddo
2834    enddo
2835  type is (real(kind=r8_kind))
2836    do k=1,size(temp,3)
2837    do j=1,size(temp,2)
2838    do i=1,size(temp,1)
2839      ind = int(dtinv*(temp(i,j,k)-tmin+teps))
2840      if (ind < 0 .or. ind > nlim) then
2841        write(unit,'(a,e10.3,a,i4,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j,k),&
2842             &'  at i=',i,' j=',j,' k=',k,' pe=',mpp_pe()
2843      endif
2844    enddo
2845    enddo
2846    enddo
2847  class default
2848    call error_mesg ('sat_vapor_pres_mod::show_all_bad_3d',&
2849         & 'The temp is not one of the supported types of real(kind=4) or real(kind=8)', FATAL)
2850  end select
2851
2852  end subroutine show_all_bad_3d
2853
2854 !#######################################################################
2855 end module sat_vapor_pres_mod
2856 !#######################################################################
2857
2858 ! <INFO>
2859
2860 !   <REFERENCE>
2861 !     Smithsonian Meteorological Tables Page 350.
2862 !   </REFERENCE>
2863
2864 !   <BUG>
2865 !     No error checking is done to make sure that the size of the
2866 !     input and output fields match.
2867 !   </BUG>
2868
2869 !   <NOTE>
2870 !     1. <B>Vectorization</B><BR/>
2871 !        To create a vector version the lookup routines need to be modified.
2872 !    The local variables: tmp, del, ind, should be changed to arrays
2873 !    with the same size and order as input array temp.
2874 !
2875 !     2. <B>Construction of the <TT>ES</TT> tables</B><BR/>
2876 !         The tables are constructed using the saturation vapor pressure (<TT>ES</TT>)
2877 !    equations in the Smithsonian tables. The tables are valid between
2878 !    -160C to +100C with increments at 1/10 degree. Between -160C and -20C
2879 !    values of <TT>ES</TT> over ice are used, between 0C and 100C values of<TT> ES</TT>
2880 !    over water are used, between -20C and 0C blended values of <TT>ES</TT>
2881 !    (over water and over ice) are used.
2882 !
2883 !    There are three tables constructed: <TT>ES</TT>, first derivative
2884 !       (<TT>ES'</TT>), and
2885 !    second derivative (<TT>ES</TT>'').  The ES table is constructed directly from
2886 !    the equations in the Smithsonian tables. The <TT>ES</TT>' table is constructed
2887 !    by bracketing temperature values at +/- 0.01 degrees. The <TT>ES</TT>'' table
2888 !    is estimated by using centered differencing of the <TT>ES</TT>' table.
2889 !
2890 !     3. <B>Determination of <TT>es</TT> and <TT>es'</TT> from lookup tables</B><BR/>
2891 !         Values of the saturation vapor pressure (<TT>es</TT>) and the
2892 !    derivative (<TT>es'</TT>) are determined at temperature (T) from the lookup
2893 !    tables (<TT>ES</TT>, <TT>ES'</TT>, <TT>ES''</TT>)
2894 !    using the following formula.
2895 !<PRE>
2896 !    es (T) = ES(t) + ES'(t) * dt + 0.5 * ES''(t) * dt**2
2897 !    es'(T) = ES'(t) + ES''(t) * dt
2898 !
2899 !    where     t = lookup table temperature closest to T
2900 !             dt = T - t
2901 !</PRE>
2902 !
2903 !     4. Internal (private) parameters<BR/>
2904 !       These parameters can be modified to increase/decrease the size/range
2905 !    of the lookup tables.
2906 !<PRE>
2907 !!    tcmin   The minimum temperature (in deg C) in the lookup tables.
2908 !!              [integer, default: tcmin = -160]
2909 !!
2910 !!    tcmax   The maximum temperature (in deg C) in the lookup tables.
2911 !!              [integer, default: tcmin = +100]
2912 !!</PRE>
2913 !!   </NOTE>
2914 !
2915 !!   <TESTPROGRAM NAME="test_sat_vapor_pres">
2916 !<PRE>
2917 !use sat_vapor_pres_mod
2918 !implicit none
2919 !
2920 !integer, parameter :: ipts=500, jpts=100, kpts=50, nloop=1
2921 !real, dimension(ipts,jpts,kpts) :: t,es,esn,des,desn
2922 !integer :: n
2923 !
2924 !! generate temperatures between 120K and 340K
2925 !  call random_number (t)
2926 !  t = 130. + t * 200.
2927 !
2928 !! initialize the tables (optional)
2929 !  call sat_vapor_pres_init
2930 !
2931 !! compute actual es and "almost" actual des
2932 !   es = compute_es  (t)
2933 !  des = compute_des (t)
2934 !
2935 !do n = 1, nloop
2936 !! es and des
2937 !  call lookup_es  (t, esn)
2938 !  call lookup_des (t,desn)
2939 !enddo
2940 !
2941 !! terminate, print deviation from actual
2942 !  print *, 'size=',ipts,jpts,kpts,nloop
2943 !  print *, 'err es  = ', sum((esn-es)**2)
2944 !  print *, 'err des = ', sum((desn-des)**2)
2945 !
2946 !contains
2947 !
2948 !!----------------------------------
2949 !! routine to estimate derivative
2950 !
2951 ! function compute_des (tem) result (des)
2952 ! real, intent(in) :: tem(:,:,:)
2953 ! real, dimension(size(tem,1),size(tem,2),size(tem,3)) :: des,esp,esm
2954 ! real, parameter :: tdel = .01
2955 !    esp = compute_es (tem+tdel)
2956 !    esm = compute_es (tem-tdel)
2957 !    des = (esp-esm)/(2*tdel)
2958 ! end function compute_des
2959 !!----------------------------------
2960 !
2961 !end program test_sat_vapor_pres
2962 !</PRE>
2963 !   </TESTPROGRAM>
2964 ! </INFO>
2965 !> @}
2966 ! close documentation grouping