LCOV - code coverage report
Current view: | top level - physics/cam - wv_saturation.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 177 | 311 | 56.9 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 18 | 29 | 62.1 %

              Line data    Source code
           1             : module wv_saturation
           2             :
           3             : !--------------------------------------------------------------------!
           4             : ! Module Overview:                                                   !
           5             : !                                                                    !
           6             : ! This module provides an interface to wv_sat_methods, providing     !
           7             : ! saturation vapor pressure and related calculations to CAM.         !
           8             : !                                                                    !
           9             : ! The original wv_saturation codes were introduced by J. J. Hack,    !
          10             : ! February 1990. The code has been extensively rewritten since then, !
          11             : ! including a total refactoring in Summer 2012.                      !
          12             : !                                                                    !
          13             : !--------------------------------------------------------------------!
          14             : ! Methods:                                                           !
          15             : !                                                                    !
          16             : ! Pure water/ice saturation vapor pressures are calculated on the    !
          17             : ! fly, with the specific method determined by a runtime option.      !
          18             : ! Mixed phase SVP is interpolated from the internal table, estbl,    !
          19             : ! which is created during initialization.                            !
          20             : !                                                                    !
          21             : ! The default method for calculating SVP is determined by a namelist !
          22             : ! option, and used whenever svp_water/ice or qsat are called.        !
          23             : !                                                                    !
          24             : !--------------------------------------------------------------------!
          25             :
          26             : use shr_kind_mod, only: r8 => shr_kind_r8
          27             : use physconst,    only: epsilo, &
          28             :                         latvap, &
          29             :                         latice, &
          30             :                         rh2o,   &
          31             :                         cpair,  &
          32             :                         tmelt,  &
          33             :                         h2otrip
          34             :
          35             : use wv_sat_methods, only: &
          36             :      svp_to_qsat => wv_sat_svp_to_qsat, &
          37             :      svp_to_qsat_vect => wv_sat_svp_to_qsat_vect
          38             :
          39             : implicit none
          40             : private
          41             : save
          42             :
          43             : ! Public interfaces
          44             : ! Namelist, initialization, finalization
          45             : public wv_sat_readnl
          46             : public wv_sat_init
          47             : public wv_sat_final
          48             :
          49             : ! Saturation vapor pressure calculations
          50             : public svp_water, svp_water_vect
          51             : public svp_ice, svp_ice_vect
          52             :
          53             : ! Mixed phase (water + ice) saturation vapor pressure table lookup
          54             : public estblf
          55             :
          56             : public svp_to_qsat
          57             :
          58             : ! Subroutines that return both SVP and humidity
          59             : ! Optional arguments do temperature derivatives
          60             : interface qsat
          61             :   module procedure qsat_line
          62             :   module procedure qsat_vect
          63             :   module procedure qsat_2D
          64             : end interface
          65             : public qsat           ! Mixed phase
          66             : interface qsat_water
          67             :   module procedure qsat_water_line
          68             :   module procedure qsat_water_vect
          69             :   module procedure qsat_water_2D
          70             : end interface
          71             : public qsat_water     ! SVP over water only
          72             : interface qsat_ice
          73             :   module procedure qsat_ice_line
          74             :   module procedure qsat_ice_vect
          75             :   module procedure qsat_ice_2D
          76             : end interface
          77             : public qsat_ice       ! SVP over ice only
          78             :
          79             : ! Wet bulb temperature solver
          80             : public :: findsp_vc, findsp
          81             :
          82             : ! Data
          83             :
          84             : ! This value is slightly high, but it seems to be the value for the
          85             : ! steam point of water originally (and most frequently) used in the
          86             : ! Goff & Gratch scheme.
          87             : real(r8), parameter :: tboil = 373.16_r8
          88             :
          89             : ! Table of saturation vapor pressure values (estbl) from tmin to
          90             : ! tmax+1 Kelvin, in one degree increments.  ttrice defines the
          91             : ! transition region, estbl contains a combination of ice & water
          92             : ! values.
          93             : ! Make these public parameters in case another module wants to see the
          94             : ! extent of the table.
          95             :   real(r8), public, parameter :: tmin = 127.16_r8
          96             :   real(r8), public, parameter :: tmax = 375.16_r8
          97             :
          98             :   real(r8), parameter :: ttrice = 20.00_r8  ! transition range from es over H2O to es over ice
          99             :
         100             :   integer :: plenest                             ! length of estbl
         101             :   real(r8), allocatable :: estbl(:)              ! table values of saturation vapor pressure
         102             :
         103             :   real(r8) :: omeps      ! 1.0_r8 - epsilo
         104             :
         105             :   real(r8) :: c3         ! parameter used by findsp
         106             :
         107             :   ! Set coefficients for polynomial approximation of difference
         108             :   ! between saturation vapor press over water and saturation pressure
         109             :   ! over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial is
         110             :   ! valid in the range -40 < t < 0 (degrees C).
         111             :   real(r8) :: pcf(5) = (/ &
         112             :        5.04469588506e-01_r8, &
         113             :        -5.47288442819e+00_r8, &
         114             :        -3.67471858735e-01_r8, &
         115             :        -8.95963532403e-03_r8, &
         116             :        -7.78053686625e-05_r8 /)
         117             :
         118             : !   --- Degree 6 approximation ---
         119             : !  real(r8) :: pcf(6) = (/ &
         120             : !       7.63285250063e-02, &
         121             : !       5.86048427932e+00, &
         122             : !       4.38660831780e-01, &
         123             : !       1.37898276415e-02, &
         124             : !       2.14444472424e-04, &
         125             : !       1.36639103771e-06 /)
         126             :
         127             :   integer, parameter :: VLENS = 128 ! vector length for a GPU kernel
         128             :
         129             :   !$acc declare create (plenest,estbl,omeps,c3,pcf)
         130             :
         131             : contains
         132             :
         133             : !---------------------------------------------------------------------
         134             : ! ADMINISTRATIVE FUNCTIONS
         135             : !---------------------------------------------------------------------
         136             :
         137        1536 : subroutine wv_sat_readnl(nlfile)
         138             :   !------------------------------------------------------------------!
         139             :   ! Purpose:                                                         !
         140             :   !   Get runtime options for wv_saturation.                         !
         141             :   !------------------------------------------------------------------!
         142             :
         143             :   use wv_sat_methods, only: wv_sat_get_scheme_idx, &
         144             :                             wv_sat_valid_idx, &
         145             :                             wv_sat_set_default
         146             :
         147             :   use spmd_utils,      only: masterproc
         148             :   use namelist_utils,  only: find_group_name
         149             :   use units,           only: getunit, freeunit
         150             :   use mpishorthand
         151             :   use cam_abortutils,  only: endrun
         152             :
         153             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         154             :
         155             :   ! Local variables
         156             :   integer :: unitn, ierr
         157             :
         158             :   character(len=32) :: wv_sat_scheme = "GoffGratch"
         159             :
         160             :   character(len=*), parameter :: subname = 'wv_sat_readnl'
         161             :
         162             :   namelist /wv_sat_nl/ wv_sat_scheme
         163             :   !-----------------------------------------------------------------------------
         164             :
         165        1536 :   if (masterproc) then
         166           2 :      unitn = getunit()
         167           2 :      open( unitn, file=trim(nlfile), status='old' )
         168           2 :      call find_group_name(unitn, 'wv_sat_nl', status=ierr)
         169           2 :      if (ierr == 0) then
         170           0 :         read(unitn, wv_sat_nl, iostat=ierr)
         171           0 :         if (ierr /= 0) then
         172           0 :            call endrun(subname // ':: ERROR reading namelist')
         173           0 :            return
         174             :         end if
         175             :      end if
         176           2 :      close(unitn)
         177           2 :      call freeunit(unitn)
         178             :
         179             :   end if
         180             :
         181             : #ifdef SPMD
         182        1536 :   call mpibcast(wv_sat_scheme, len(wv_sat_scheme) , mpichar, 0, mpicom)
         183             : #endif
         184             :
         185        1536 :   if (.not. wv_sat_set_default(wv_sat_scheme)) then
         186           0 :      call endrun('wv_sat_readnl :: Invalid wv_sat_scheme.')
         187           0 :      return
         188             :   end if
         189             :
         190             : end subroutine wv_sat_readnl
         191             :
         192        1536 : subroutine wv_sat_init
         193             :   !------------------------------------------------------------------!
         194             :   ! Purpose:                                                         !
         195             :   !   Initialize module (e.g. setting parameters, initializing the   !
         196             :   !   SVP lookup table).                                             !
         197             :   !------------------------------------------------------------------!
         198             :
         199             :   use wv_sat_methods, only: wv_sat_methods_init, &
         200             :                             wv_sat_get_scheme_idx, &
         201             :                             wv_sat_valid_idx
         202             :   use spmd_utils,     only: masterproc
         203             :   use cam_logfile,    only: iulog
         204             :   use cam_abortutils, only: endrun
         205             :   use shr_assert_mod, only: shr_assert_in_domain
         206             :   use error_messages, only: handle_errmsg
         207             :
         208             :   integer :: status
         209             :
         210             :   ! For wv_sat_methods error reporting.
         211             :   character(len=256) :: errstring
         212             :
         213             :   ! For generating internal SVP table.
         214             :   real(r8) :: t         ! Temperature
         215             :   integer  :: i         ! Increment counter
         216             :
         217             :   ! Precalculated because so frequently used.
         218        1536 :   omeps  = 1.0_r8 - epsilo
         219             :
         220             :   ! Transition range method is only valid for transition temperatures at:
         221             :   ! -40 deg C < T < 0 deg C
         222             :   call shr_assert_in_domain(ttrice, ge=0._r8, le=40._r8, varname="ttrice",&
         223        1536 :        msg="wv_sat_init: Invalid transition temperature range.")
         224             :
         225             : ! This parameter uses a hardcoded 287.04_r8?
         226        1536 :   c3 = 287.04_r8*(7.5_r8*log(10._r8))/cpair
         227             :
         228             : ! Init "methods" module containing actual SVP formulae.
         229             :
         230             :   call wv_sat_methods_init(r8, tmelt, h2otrip, tboil, ttrice, &
         231        1536 :        epsilo, errstring)
         232             :
         233        1536 :   call handle_errmsg(errstring, subname="wv_sat_methods_init")
         234             :
         235             :   ! Add two to make the table slightly too big, just in case.
         236        1536 :   plenest = ceiling(tmax-tmin) + 2
         237             :
         238             :   ! Allocate SVP table.
         239        1536 :   allocate(estbl(plenest), stat=status)
         240        1536 :   if (status /= 0) then
         241           0 :      call endrun('wv_sat_init :: ERROR allocating saturation vapor pressure table')
         242           0 :      return
         243             :   end if
         244             :
         245      387072 :   do i = 1, plenest
         246      387072 :      estbl(i) = svp_trans(tmin + real(i-1,r8))
         247             :   end do
         248             :
         249             :   !$acc update device (plenest,estbl,omeps,c3,pcf)
         250             :
         251        1536 :   if (masterproc) then
         252           2 :      write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
         253             :   end if
         254             :
         255             : end subroutine wv_sat_init
         256             :
         257        1536 : subroutine wv_sat_final
         258             :   !------------------------------------------------------------------!
         259             :   ! Purpose:                                                         !
         260             :   !   Deallocate global variables in module.                         !
         261             :   !------------------------------------------------------------------!
         262             :   use cam_abortutils, only: endrun
         263             :
         264             :   integer :: status
         265             :
         266        1536 :   if (allocated(estbl)) then
         267             :
         268        1536 :      deallocate(estbl, stat=status)
         269             :
         270             :      if (status /= 0) then
         271             :         call endrun('wv_sat_final :: ERROR deallocating table')
         272             :         return
         273             :      end if
         274             :
         275             :   end if
         276             :
         277             : end subroutine wv_sat_final
         278             :
         279             : !---------------------------------------------------------------------
         280             : ! DEFAULT SVP FUNCTIONS
         281             : !---------------------------------------------------------------------
         282             :
         283             : ! Compute saturation vapor pressure over water
         284           0 : function svp_water(t) result(es)
         285             :
         286             :   use wv_sat_methods, only: &
         287             :        wv_sat_svp_water
         288             :
         289             :   real(r8), intent(in)  :: t  ! Temperature (K)
         290             :   real(r8)              :: es ! SVP (Pa)
         291             :
         292           0 :   es = wv_sat_svp_water(t)
         293             :
         294           0 : end function svp_water
         295             :
         296             : ! Compute saturation vapor pressure over water
         297           0 : subroutine svp_water_vect(t, es, vlen)
         298             :
         299             :   use wv_sat_methods, only: &
         300             :        wv_sat_svp_water_vect
         301             :
         302             :   integer,  intent(in)  :: vlen
         303             :   real(r8), intent(in)  :: t(vlen)  ! Temperature (K)
         304             :   real(r8), intent(out) :: es(vlen) ! SVP (Pa)
         305             :
         306             :   !$acc data copyin (t) copyout (es)
         307             :
         308           0 :   call wv_sat_svp_water_vect(t, es, vlen)
         309             :
         310             :   !$acc end data
         311           0 : end subroutine svp_water_vect
         312             :
         313             : ! Compute saturation vapor pressure over ice
         314           0 : function svp_ice(t) result(es)
         315             :
         316             :   use wv_sat_methods, only: &
         317             :        wv_sat_svp_ice
         318             :
         319             :   real(r8), intent(in)  :: t  ! Temperature (K)
         320             :   real(r8)              :: es ! SVP (Pa)
         321             :
         322           0 :   es = wv_sat_svp_ice(t)
         323             :
         324           0 : end function svp_ice
         325             :
         326             : ! Compute saturation vapor pressure over ice
         327           0 : subroutine svp_ice_vect(t, es, vlen)
         328             :
         329             :   use wv_sat_methods, only: &
         330             :        wv_sat_svp_ice_vect
         331             :
         332             :   integer,  intent(in)  :: vlen
         333             :   real(r8), intent(in)  :: t(vlen)  ! Temperature (K)
         334             :   real(r8), intent(out) :: es(vlen) ! SVP (Pa)
         335             :
         336             :   !$acc data copyin(t) copyout(es)
         337             :
         338           0 :   call wv_sat_svp_ice_vect(t, es, vlen)
         339             :
         340             :   !$acc end data
         341           0 : end subroutine svp_ice_vect
         342             :
         343             : ! Compute saturation vapor pressure with an ice-water transition
         344      385536 : function svp_trans(t) result(es)
         345             :
         346             :   use wv_sat_methods, only: &
         347             :        wv_sat_svp_trans
         348             :
         349             :   real(r8), intent(in)  :: t   ! Temperature (K)
         350             :   real(r8)              :: es  ! SVP (Pa)
         351             :
         352      385536 :   es = wv_sat_svp_trans(t)
         353             :
         354      385536 : end function svp_trans
         355             :
         356             : ! Compute saturation vapor pressure with an ice-water transition
         357             : subroutine svp_trans_vect(t, es, vlen)
         358             :
         359             :   use wv_sat_methods, only: &
         360             :        wv_sat_svp_trans_vect
         361             :
         362             :   integer,  intent(in)  :: vlen
         363             :   real(r8), intent(in)  :: t(vlen)   ! Temperature (K)
         364             :   real(r8), intent(out) :: es(vlen)  ! SVP (Pa)
         365             :
         366             :   !$acc data copyin(t) copyout(es)
         367             :
         368             :   call wv_sat_svp_trans_vect(t, es, vlen)
         369             :
         370             :   !$acc end data
         371             : end subroutine svp_trans_vect
         372             :
         373             : !---------------------------------------------------------------------
         374             : ! UTILITIES
         375             : !---------------------------------------------------------------------
         376             :
         377             : ! Does linear interpolation from nearest values found
         378             : ! in the table (estbl).
         379  4913471668 : elemental function estblf(t) result(es)
         380             :
         381             :   real(r8), intent(in) :: t ! Temperature
         382             :   real(r8) :: es            ! SVP (Pa)
         383             :
         384             :   integer  :: i         ! Index for t in the table
         385             :   real(r8) :: t_tmp     ! intermediate temperature for es look-up
         386             :
         387             :   real(r8) :: weight ! Weight for interpolation
         388             :
         389  4913471668 :   t_tmp = max(min(t,tmax)-tmin, 0._r8)   ! Number of table entries above tmin
         390  4913471668 :   i = int(t_tmp) + 1                     ! Corresponding index.
         391  4913471668 :   weight = t_tmp - aint(t_tmp, r8)       ! Fractional part of t_tmp (for interpolation).
         392  4913471668 :   es = (1._r8 - weight)*estbl(i) + weight*estbl(i+1)
         393             :
         394  4913471668 : end function estblf
         395             :
         396             : ! Does linear interpolation from nearest values found
         397             : ! in the table (estbl).
         398   486229471 : subroutine estblf_vect(t, es, vlen)
         399             :
         400             :   integer,                   intent(in)  :: vlen
         401             :   real(r8), dimension(vlen), intent(in)  :: t     ! Temperature
         402             :   real(r8), dimension(vlen), intent(out) :: es    ! SVP (Pa)
         403             :
         404             :   integer  :: i         ! Index for t in the table
         405             :   integer  :: j
         406             :   real(r8) :: t_tmp     ! intermediate temperature for es look-up
         407             :
         408             :   real(r8) :: weight ! Weight for interpolation
         409             :
         410             :   !$acc data present (t,es)
         411             :
         412             :   !$acc parallel vector_length(VLENS) default(present)
         413             :   !$acc loop gang vector private(t_tmp,weight,i)
         414  8104683799 :   do j = 1, vlen
         415  7618454328 :      t_tmp = max(min(t(j),tmax)-tmin, 0._r8)   ! Number of table entries above tmin
         416  7618454328 :      i = int(t_tmp) + 1                     ! Corresponding index.
         417  7618454328 :      weight = t_tmp - aint(t_tmp, r8)       ! Fractional part of t_tmp (for interpolation).
         418  8104683799 :      es(j) = (1._r8 - weight)*estbl(i) + weight*estbl(i+1)
         419             :   end do
         420             :   !$acc end parallel
         421             :
         422             :   !$acc end data
         423   486229471 : end subroutine estblf_vect
         424             :
         425             : ! Get enthalpy based only on temperature
         426             : ! and specific humidity.
         427  2472156146 : elemental function tq_enthalpy(t, q, hltalt) result(enthalpy)
         428             :
         429             :   real(r8), intent(in) :: t      ! Temperature
         430             :   real(r8), intent(in) :: q      ! Specific humidity
         431             :   real(r8), intent(in) :: hltalt ! Modified hlat for T derivatives
         432             :
         433             :   real(r8) :: enthalpy
         434             :
         435  2472156146 :   enthalpy = cpair * t + hltalt * q
         436             :
         437  2472156146 : end function tq_enthalpy
         438             :
         439             : ! Get enthalpy based only on temperature
         440             : ! and specific humidity.
         441           0 : subroutine tq_enthalpy_vect(t, q, hltalt, enthalpy, vlen)
         442             :
         443             :   integer,                   intent(in)  :: vlen
         444             :   real(r8), dimension(vlen), intent(in)  :: t      ! Temperature
         445             :   real(r8), dimension(vlen), intent(in)  :: q      ! Specific humidity
         446             :   real(r8), dimension(vlen), intent(in)  :: hltalt ! Modified hlat for T derivatives
         447             :
         448             :   real(r8), dimension(vlen), intent(out) :: enthalpy
         449             :
         450             :   integer :: i
         451             :
         452             :   !$acc data present(t,q,hltalt,enthalpy)
         453             :
         454             :   !$acc parallel vector_length(VLENS) default(present)
         455             :   !$acc loop gang vector
         456           0 :   do i = 1, vlen
         457           0 :      enthalpy(i) = cpair * t(i) + hltalt(i) * q(i)
         458             :   end do
         459             :   !$acc end parallel
         460             :
         461             :   !$acc end data
         462           0 : end subroutine tq_enthalpy_vect
         463             :
         464             : !---------------------------------------------------------------------
         465             : ! LATENT HEAT OF VAPORIZATION CORRECTIONS
         466             : !---------------------------------------------------------------------
         467             :
         468  3692793746 : elemental subroutine no_ip_hltalt(t, hltalt)
         469             :   !------------------------------------------------------------------!
         470             :   ! Purpose:                                                         !
         471             :   !   Calculate latent heat of vaporization of pure liquid water at  !
         472             :   !   a given temperature.                                           !
         473             :   !------------------------------------------------------------------!
         474             :
         475             :   ! Inputs
         476             :   real(r8), intent(in) :: t        ! Temperature
         477             :   ! Outputs
         478             :   real(r8), intent(out) :: hltalt  ! Appropriately modified hlat
         479             :
         480  3692793746 :   hltalt = latvap
         481             :
         482             :   ! Account for change of latvap with t above freezing where
         483             :   ! constant slope is given by -2369 j/(kg c) = cpv - cw
         484  3692793746 :   if (t >= tmelt) then
         485   586171924 :      hltalt = hltalt - 2369.0_r8*(t-tmelt)
         486             :   end if
         487             :
         488  3692793746 : end subroutine no_ip_hltalt
         489             :
         490    76786567 : subroutine no_ip_hltalt_vect(t, hltalt, vlen)
         491             :   !------------------------------------------------------------------!
         492             :   ! Purpose:                                                         !
         493             :   !   Calculate latent heat of vaporization of pure liquid water at  !
         494             :   !   a given temperature.                                           !
         495             :   !------------------------------------------------------------------!
         496             :
         497             :   ! Inputs
         498             :   integer,                   intent(in) :: vlen
         499             :   real(r8), dimension(vlen), intent(in) :: t        ! Temperature
         500             :   ! Outputs
         501             :   real(r8), dimension(vlen), intent(out) :: hltalt  ! Appropriately modified hlat
         502             :
         503             :   integer :: i
         504             :
         505             :   !$acc data present(t,hltalt)
         506             :
         507             :   !$acc parallel vector_length(VLENS) default(present)
         508             :   !$acc loop gang vector
         509  1267939495 :   do i = 1, vlen
         510  1191152928 :      hltalt(i) = latvap
         511             :      ! Account for change of latvap with t above freezing where
         512             :      ! constant slope is given by -2369 j/(kg c) = cpv - cw
         513  1267939495 :      if (t(i) >= tmelt) then
         514   227409119 :         hltalt(i) = hltalt(i) - 2369.0_r8*(t(i)-tmelt)
         515             :      end if
         516             :   end do
         517             :   !$acc end parallel
         518             :
         519             :   !$acc end data
         520    76786567 : end subroutine no_ip_hltalt_vect
         521             :
         522  3692793746 : elemental subroutine calc_hltalt(t, hltalt, tterm)
         523             :   !------------------------------------------------------------------!
         524             :   ! Purpose:                                                         !
         525             :   !   Calculate latent heat of vaporization of water at a given      !
         526             :   !   temperature, taking into account the ice phase if temperature  !
         527             :   !   is below freezing.                                             !
         528             :   !   Optional argument also calculates a term used to calculate     !
         529             :   !   d(es)/dT within the water-ice transition range.                !
         530             :   !------------------------------------------------------------------!
         531             :
         532             :   ! Inputs
         533             :   real(r8), intent(in) :: t        ! Temperature
         534             :   ! Outputs
         535             :   real(r8), intent(out) :: hltalt  ! Appropriately modified hlat
         536             :   ! Term to account for d(es)/dT in transition region.
         537             :   real(r8), intent(out), optional :: tterm
         538             :
         539             :   ! Local variables
         540             :   real(r8) :: tc      ! Temperature in degrees C
         541             :   real(r8) :: weight  ! Weight for es transition from water to ice
         542             :   ! Loop iterator
         543             :   integer :: i
         544             :
         545  3692793746 :   if (present(tterm)) tterm = 0.0_r8
         546             :
         547  3692793746 :   call no_ip_hltalt(t,hltalt)
         548  3692793746 :   if (t < tmelt) then
         549             :      ! Weighting of hlat accounts for transition from water to ice.
         550  3106621822 :      tc = t - tmelt
         551             :
         552  3106621822 :      if (tc >= -ttrice) then
         553   577722468 :         weight = -tc/ttrice
         554             :
         555             :         ! polynomial expression approximates difference between es
         556             :         ! over water and es over ice from 0 to -ttrice (C) (max of
         557             :         ! ttrice is 40): required for accurate estimate of es
         558             :         ! derivative in transition range from ice to water
         559   577722468 :         if (present(tterm)) then
         560  2993778546 :            do i = size(pcf), 1, -1
         561  2993778546 :               tterm = pcf(i) + tc*tterm
         562             :            end do
         563   498963091 :            tterm = tterm/ttrice
         564             :         end if
         565             :
         566             :      else
         567             :         weight = 1.0_r8
         568             :      end if
         569             :
         570  3106621822 :      hltalt = hltalt + weight*latice
         571             :
         572             :   end if
         573             :
         574  3692793746 : end subroutine calc_hltalt
         575             :
         576    76786567 : subroutine calc_hltalt_vect(t, hltalt, vlen, tterm)
         577             :   !------------------------------------------------------------------!
         578             :   ! Purpose:                                                         !
         579             :   !   Calculate latent heat of vaporization of water at a given      !
         580             :   !   temperature, taking into account the ice phase if temperature  !
         581             :   !   is below freezing.                                             !
         582             :   !   Optional argument also calculates a term used to calculate     !
         583             :   !   d(es)/dT within the water-ice transition range.                !
         584             :   !------------------------------------------------------------------!
         585             :
         586             :   ! Inputs
         587             :   integer,                   intent(in) :: vlen
         588             :   real(r8), dimension(vlen), intent(in) :: t        ! Temperature
         589             :   ! Outputs
         590             :   real(r8), dimension(vlen), intent(out) :: hltalt  ! Appropriately modified hlat
         591             :   ! Term to account for d(es)/dT in transition region.
         592             :   real(r8), dimension(vlen), intent(out), optional :: tterm
         593             :
         594             :   ! Local variables
         595             :   real(r8) :: tc      ! Temperature in degrees C
         596             :   real(r8) :: weight  ! Weight for es transition from water to ice
         597             :   logical  :: present_tterm
         598             :   ! Loop iterator
         599             :   integer :: i, j, size_pcf
         600             :
         601    76786567 :   present_tterm = present(tterm)
         602    76786567 :   size_pcf      = size(pcf)
         603             :
         604             :   !$acc data present(t,hltalt,tterm)
         605             :
         606    76786567 :   if (present_tterm) then
         607             :      !$acc parallel vector_length(VLENS) default(present)
         608             :      !$acc loop gang vector
         609  1267939495 :      do i = 1, vlen
         610  1267939495 :         tterm(i) = 0.0_r8
         611             :      end do
         612             :      !$acc end parallel
         613             :   end if
         614             :
         615    76786567 :   call no_ip_hltalt_vect(t,hltalt,vlen)
         616             :
         617             :   !$acc parallel vector_length(VLENS) default(present)
         618             :   !$acc loop gang vector private(tc,weight)
         619  1267939495 :   do j = 1, vlen
         620  1267939495 :      if (t(j) < tmelt) then
         621             :         ! Weighting of hlat accounts for transition from water to ice.
         622   963743809 :         tc = t(j) - tmelt
         623             :
         624   963743809 :         if (tc >= -ttrice) then
         625   180643049 :            weight = -tc/ttrice
         626             :
         627             :            ! polynomial expression approximates difference between es
         628             :            ! over water and es over ice from 0 to -ttrice (C) (max of
         629             :            ! ttrice is 40): required for accurate estimate of es
         630             :            ! derivative in transition range from ice to water
         631   180643049 :            if (present_tterm) then
         632             :               !$acc loop seq
         633  1083858294 :               do i = size_pcf, 1, -1
         634  1083858294 :                  tterm(j) = pcf(i) + tc*tterm(j)
         635             :               end do
         636   180643049 :               tterm(j) = tterm(j)/ttrice
         637             :            end if
         638             :
         639             :         else
         640             :            weight = 1.0_r8
         641             :         end if
         642             :
         643   963743809 :         hltalt(j) = hltalt(j) + weight*latice
         644             :
         645             :      end if
         646             :   end do
         647             :   !$acc end parallel
         648             :
         649             :   !$acc end data
         650    76786567 : end subroutine calc_hltalt_vect
         651             :
         652             : !---------------------------------------------------------------------
         653             : ! OPTIONAL OUTPUTS
         654             : !---------------------------------------------------------------------
         655             :
         656             : ! Temperature derivative outputs, for qsat_*
         657  3082515268 : subroutine deriv_outputs_line(t, p, es, qs, hltalt, tterm, &
         658             :      gam, dqsdt)
         659             :
         660             :   ! Inputs
         661             :   real(r8), intent(in) :: t      ! Temperature
         662             :   real(r8), intent(in) :: p      ! Pressure
         663             :   real(r8), intent(in) :: es     ! Saturation vapor pressure
         664             :   real(r8), intent(in) :: qs     ! Saturation specific humidity
         665             :   real(r8), intent(in) :: hltalt ! Modified latent heat
         666             :   real(r8), intent(in) :: tterm  ! Extra term for d(es)/dT in
         667             :                                  ! transition region.
         668             :
         669             :   ! Outputs
         670             :   real(r8), intent(out), optional :: gam      ! (hltalt/cpair)*(d(qs)/dt)
         671             :   real(r8), intent(out), optional :: dqsdt    ! (d(qs)/dt)
         672             :
         673             :   ! Local variables
         674             :   real(r8) :: desdt        ! d(es)/dt
         675             :   real(r8) :: dqsdt_loc    ! local copy of dqsdt
         676             :
         677  3082515268 :   if (qs == 1.0_r8) then
         678             :      dqsdt_loc = 0._r8
         679             :   else
         680  3082515268 :      desdt = hltalt*es/(rh2o*t*t) + tterm
         681  3082515268 :      dqsdt_loc = qs*p*desdt/(es*(p-omeps*es))
         682             :   end if
         683             :
         684  3082515268 :   if (present(dqsdt)) dqsdt = dqsdt_loc
         685  3082515268 :   if (present(gam))   gam   = dqsdt_loc * (hltalt/cpair)
         686             :
         687  3082515268 : end subroutine deriv_outputs_line
         688             :
         689             : ! Temperature derivative outputs, for qsat_*
         690   153573134 : subroutine deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
         691    76786567 :      gam, dqsdt)
         692             :
         693             :   ! Inputs
         694             :   integer,                   intent(in) :: vlen
         695             :   real(r8), dimension(vlen), intent(in) :: t      ! Temperature
         696             :   real(r8), dimension(vlen), intent(in) :: p      ! Pressure
         697             :   real(r8), dimension(vlen), intent(in) :: es     ! Saturation vapor pressure
         698             :   real(r8), dimension(vlen), intent(in) :: qs     ! Saturation specific humidity
         699             :   real(r8), dimension(vlen), intent(in) :: hltalt ! Modified latent heat
         700             :   real(r8), dimension(vlen), intent(in) :: tterm  ! Extra term for d(es)/dT in
         701             :                                  ! transition region.
         702             :
         703             :   ! Outputs
         704             :   real(r8), dimension(vlen), intent(out), optional :: gam      ! (hltalt/cpair)*(d(qs)/dt)
         705             :   real(r8), dimension(vlen), intent(out), optional :: dqsdt    ! (d(qs)/dt)
         706             :
         707             :   ! Local variables
         708             :   real(r8) :: desdt        ! d(es)/dt
         709             :   real(r8) :: dqsdt_loc    ! local copy of dqsdt
         710             :   logical  :: present_dqsdt, present_gam
         711             :   integer  :: i
         712             :
         713    76786567 :   present_dqsdt = present(dqsdt)
         714    76786567 :   present_gam   = present(gam)
         715             :
         716             :   !$acc data present(t,p,es,qs,hltalt,tterm,gam,dqsdt)
         717             :
         718             :   !$acc parallel vector_length(VLENS) default(present)
         719             :   !$acc loop gang vector private(dqsdt_loc,desdt)
         720  1267939495 :   do i = 1, vlen
         721  1191152928 :      if (qs(i) == 1.0_r8) then
         722             :         dqsdt_loc = 0._r8
         723             :      else
         724  1191152928 :         desdt = hltalt(i)*es(i)/(rh2o*t(i)*t(i)) + tterm(i)
         725  1191152928 :         dqsdt_loc = qs(i)*p(i)*desdt/(es(i)*(p(i)-omeps*es(i)))
         726             :      end if
         727             :
         728  1191152928 :      if (present_dqsdt) dqsdt(i) = dqsdt_loc
         729  1267939495 :      if (present_gam)   gam(i)   = dqsdt_loc * (hltalt(i)/cpair)
         730             :   end do
         731             :   !$acc end parallel
         732             :
         733             :   !$acc end data
         734    76786567 : end subroutine deriv_outputs_vect
         735             :
         736             : !---------------------------------------------------------------------
         737             : ! QSAT (SPECIFIC HUMIDITY) PROCEDURES
         738             : !---------------------------------------------------------------------
         739             :
         740  3692834068 : subroutine qsat_line(t, p, es, qs, gam, dqsdt, enthalpy)
         741             :   !------------------------------------------------------------------!
         742             :   ! Purpose:                                                         !
         743             :   !   Look up and return saturation vapor pressure from precomputed  !
         744             :   !   table, then calculate and return saturation specific humidity. !
         745             :   !   Optionally return various temperature derivatives or enthalpy  !
         746             :   !   at saturation.                                                 !
         747             :   !------------------------------------------------------------------!
         748             :
         749             :   ! Inputs
         750             :   real(r8), intent(in) :: t    ! Temperature
         751             :   real(r8), intent(in) :: p    ! Pressure
         752             :   ! Outputs
         753             :   real(r8), intent(out) :: es  ! Saturation vapor pressure
         754             :   real(r8), intent(out) :: qs  ! Saturation specific humidity
         755             :
         756             :   real(r8), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
         757             :   real(r8), intent(out), optional :: dqsdt  ! (d(qs)/dt)
         758             :   real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
         759             :
         760             :   ! Local variables
         761             :   real(r8) :: hltalt       ! Modified latent heat for T derivatives
         762             :   real(r8) :: tterm        ! Account for d(es)/dT in transition region
         763             :
         764  3692834068 :   es = estblf(t)
         765             :
         766  3692834068 :   qs = svp_to_qsat(es, p)
         767             :
         768             :   ! Ensures returned es is consistent with limiters on qs.
         769  3692834068 :   es = min(es, p)
         770             :
         771             :   ! Calculate optional arguments.
         772  3692834068 :   if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
         773             :
         774             :      ! "generalized" analytic expression for t derivative of es
         775             :      ! accurate to within 1 percent for 173.16 < t < 373.16
         776  3082515268 :      call calc_hltalt(t, hltalt, tterm)
         777             :
         778  3082515268 :      if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
         779             :
         780             :      call deriv_outputs_line(t, p, es, qs, hltalt, tterm, &
         781  3082515268 :           gam=gam, dqsdt=dqsdt)
         782             :
         783             :   end if
         784             :
         785  3692834068 : end subroutine qsat_line
         786             :
         787   486229471 : subroutine qsat_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
         788             :   !------------------------------------------------------------------!
         789             :   ! Purpose:                                                         !
         790             :   !   Look up and return saturation vapor pressure from precomputed  !
         791             :   !   table, then calculate and return saturation specific humidity. !
         792             :   !   Optionally return various temperature derivatives or enthalpy  !
         793             :   !   at saturation.                                                 !
         794             :   !------------------------------------------------------------------!
         795             :
         796             :   ! Inputs
         797             :   integer,                   intent(in) :: vlen
         798             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
         799             :   real(r8), dimension(vlen), intent(in) :: p    ! Pressure
         800             :   ! Outputs
         801             :   real(r8), dimension(vlen), intent(out) :: es  ! Saturation vapor pressure
         802             :   real(r8), dimension(vlen), intent(out) :: qs  ! Saturation specific humidity
         803             :
         804             :   real(r8), dimension(vlen), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
         805             :   real(r8), dimension(vlen), intent(out), optional :: dqsdt  ! (d(qs)/dt)
         806             :   real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
         807             :
         808             :   ! Local variables
         809   972458942 :   real(r8), dimension(vlen) :: hltalt       ! Modified latent heat for T derivatives
         810   972458942 :   real(r8), dimension(vlen) :: tterm        ! Account for d(es)/dT in transition region
         811             :   integer                   :: i
         812             :   logical                   :: present_gam, present_dqsdt, present_enthalpy
         813             :
         814   486229471 :   present_gam = present(gam)
         815   486229471 :   present_dqsdt = present(dqsdt)
         816   486229471 :   present_enthalpy = present(enthalpy)
         817             :
         818             :   !$acc data copyin  (t,p) &
         819             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
         820             :   !$acc      create  (hltalt,tterm)
         821             :
         822   486229471 :   call estblf_vect(t, es, vlen)
         823             :
         824   486229471 :   call svp_to_qsat_vect(es, p, qs, vlen)
         825             :
         826             :   ! Ensures returned es is consistent with limiters on qs.
         827             :
         828             :   !$acc parallel vector_length(VLENS) default(present)
         829             :   !$acc loop gang vector
         830  8104683799 :   do i = 1, vlen
         831  8104683799 :      es(i) = min(es(i), p(i))
         832             :   end do
         833             :   !$acc end parallel
         834             :
         835             :   ! Calculate optional arguments.
         836   486229471 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
         837             :
         838             :      ! "generalized" analytic expression for t derivative of es
         839             :      ! accurate to within 1 percent for 173.16 < t < 373.16
         840    76786567 :      call calc_hltalt_vect(t, hltalt, vlen, tterm)
         841             :
         842    76786567 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
         843             :
         844             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
         845    76786567 :           gam=gam, dqsdt=dqsdt)
         846             :
         847             :   end if
         848             :
         849             :   !$acc end data
         850   486229471 : end subroutine qsat_vect
         851             :
         852           0 : subroutine qsat_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
         853             :   !------------------------------------------------------------------!
         854             :   ! Purpose:                                                         !
         855             :   !   Look up and return saturation vapor pressure from precomputed  !
         856             :   !   table, then calculate and return saturation specific humidity. !
         857             :   !   Optionally return various temperature derivatives or enthalpy  !
         858             :   !   at saturation.                                                 !
         859             :   !------------------------------------------------------------------!
         860             :
         861             :   ! Inputs
         862             :   integer,                        intent(in) :: dim1, dim2
         863             :   real(r8), dimension(dim1,dim2), intent(in) :: t    ! Temperature
         864             :   real(r8), dimension(dim1,dim2), intent(in) :: p    ! Pressure
         865             :   ! Outputs
         866             :   real(r8), dimension(dim1,dim2), intent(out) :: es  ! Saturation vapor pressure
         867             :   real(r8), dimension(dim1,dim2), intent(out) :: qs  ! Saturation specific humidity
         868             :
         869             :   real(r8), dimension(dim1,dim2), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
         870             :   real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt  ! (d(qs)/dt)
         871             :   real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
         872             :
         873             :   ! Local variables
         874           0 :   real(r8), dimension(dim1,dim2) :: hltalt       ! Modified latent heat for T derivatives
         875           0 :   real(r8), dimension(dim1,dim2) :: tterm        ! Account for d(es)/dT in transition region
         876             :   integer                        :: i, k, vlen
         877             :   logical                        :: present_gam, present_dqsdt, present_enthalpy
         878             :
         879           0 :   vlen = dim1 * dim2
         880           0 :   present_gam = present(gam)
         881           0 :   present_dqsdt = present(dqsdt)
         882           0 :   present_enthalpy = present(enthalpy)
         883             :
         884             :   !$acc data copyin  (t,p) &
         885             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
         886             :   !$acc      create  (hltalt,tterm)
         887             :
         888           0 :   call estblf_vect(t, es, vlen)
         889             :
         890           0 :   call svp_to_qsat_vect(es, p, qs, vlen)
         891             :
         892             :   ! Ensures returned es is consistent with limiters on qs.
         893             :
         894             :   !$acc parallel vector_length(VLENS) default(present)
         895             :   !$acc loop gang vector collapse(2)
         896           0 :   do k = 1, dim2
         897           0 :      do i = 1, dim1
         898           0 :         es(i,k) = min(es(i,k), p(i,k))
         899             :      end do
         900             :   end do
         901             :   !$acc end parallel
         902             :
         903             :   ! Calculate optional arguments.
         904           0 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
         905             :
         906             :      ! "generalized" analytic expression for t derivative of es
         907             :      ! accurate to within 1 percent for 173.16 < t < 373.16
         908           0 :      call calc_hltalt_vect(t, hltalt, vlen, tterm)
         909             :
         910           0 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
         911             :
         912             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
         913           0 :           gam=gam, dqsdt=dqsdt)
         914             :
         915             :   end if
         916             :
         917             :   !$acc end data
         918           0 : end subroutine qsat_2D
         919             :
         920 12128668192 : subroutine qsat_water_line(t, p, es, qs, gam, dqsdt, enthalpy)
         921             :   !------------------------------------------------------------------!
         922             :   ! Purpose:                                                         !
         923             :   !   Calculate SVP over water at a given temperature, and then      !
         924             :   !   calculate and return saturation specific humidity.             !
         925             :   !   Optionally return various temperature derivatives or enthalpy  !
         926             :   !   at saturation.                                                 !
         927             :   !------------------------------------------------------------------!
         928             :
         929             :   use wv_sat_methods, only: wv_sat_qsat_water
         930             :
         931             :   ! Inputs
         932             :   real(r8), intent(in) :: t    ! Temperature
         933             :   real(r8), intent(in) :: p    ! Pressure
         934             :   ! Outputs
         935             :   real(r8), intent(out) :: es  ! Saturation vapor pressure
         936             :   real(r8), intent(out) :: qs  ! Saturation specific humidity
         937             :
         938             :   real(r8), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
         939             :   real(r8), intent(out), optional :: dqsdt  ! (d(qs)/dt)
         940             :   real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
         941             :
         942             :   ! Local variables
         943             :   real(r8) :: hltalt       ! Modified latent heat for T derivatives
         944             :
         945 12128668192 :   call wv_sat_qsat_water(t, p, es, qs)
         946             :
         947 12128668192 :   if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
         948             :
         949             :      ! "generalized" analytic expression for t derivative of es
         950             :      ! accurate to within 1 percent for 173.16 < t < 373.16
         951           0 :      call no_ip_hltalt(t, hltalt)
         952             :
         953           0 :      if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
         954             :
         955             :      ! For pure water/ice transition term is 0.
         956             :      call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, &
         957           0 :           gam=gam, dqsdt=dqsdt)
         958             :
         959             :   end if
         960             :
         961 12128668192 : end subroutine qsat_water_line
         962             :
         963           0 : subroutine qsat_water_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
         964             :   !------------------------------------------------------------------!
         965             :   ! Purpose:                                                         !
         966             :   !   Calculate SVP over water at a given temperature, and then      !
         967             :   !   calculate and return saturation specific humidity.             !
         968             :   !   Optionally return various temperature derivatives or enthalpy  !
         969             :   !   at saturation.                                                 !
         970             :   !------------------------------------------------------------------!
         971             :
         972             :   use wv_sat_methods, only: wv_sat_qsat_water_vect
         973             :
         974             :   ! Inputs
         975             :   integer,                   intent(in) :: vlen
         976             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
         977             :   real(r8), dimension(vlen), intent(in) :: p    ! Pressure
         978             :   ! Outputs
         979             :   real(r8), dimension(vlen), intent(out) :: es  ! Saturation vapor pressure
         980             :   real(r8), dimension(vlen), intent(out) :: qs  ! Saturation specific humidity
         981             :
         982             :   real(r8), dimension(vlen), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
         983             :   real(r8), dimension(vlen), intent(out), optional :: dqsdt  ! (d(qs)/dt)
         984             :   real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
         985             :
         986             :   ! Local variables
         987           0 :   real(r8), dimension(vlen) :: hltalt       ! Modified latent heat for T derivatives
         988           0 :   real(r8), dimension(vlen) :: tterm
         989             :   integer                   :: i
         990             :   logical                   :: present_gam, present_dqsdt, present_enthalpy
         991             :
         992           0 :   present_gam      = present(gam)
         993           0 :   present_dqsdt    = present(dqsdt)
         994           0 :   present_enthalpy = present(enthalpy)
         995             :
         996             :   !$acc data copyin  (t,p) &
         997             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
         998             :   !$acc      create  (tterm,hltalt)
         999             :
        1000             :   !$acc parallel vector_length(VLENS) default(present)
        1001             :   !$acc loop gang vector
        1002           0 :   do i = 1, vlen
        1003           0 :      tterm(i) = 0._r8
        1004             :   end do
        1005             :   !$acc end parallel
        1006             :
        1007           0 :   call wv_sat_qsat_water_vect(t, p, es, qs, vlen)
        1008             :
        1009           0 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
        1010             :
        1011             :      ! "generalized" analytic expression for t derivative of es
        1012             :      ! accurate to within 1 percent for 173.16 < t < 373.16
        1013           0 :      call no_ip_hltalt_vect(t, hltalt, vlen)
        1014             :
        1015           0 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
        1016             :
        1017             :      ! For pure water/ice transition term is 0.
        1018             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
        1019           0 :           gam=gam, dqsdt=dqsdt)
        1020             :
        1021             :   end if
        1022             :
        1023             :   !$acc end data
        1024           0 : end subroutine qsat_water_vect
        1025             :
        1026           0 : subroutine qsat_water_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
        1027             :   !------------------------------------------------------------------!
        1028             :   ! Purpose:                                                         !
        1029             :   !   Calculate SVP over water at a given temperature, and then      !
        1030             :   !   calculate and return saturation specific humidity.             !
        1031             :   !   Optionally return various temperature derivatives or enthalpy  !
        1032             :   !   at saturation.                                                 !
        1033             :   !------------------------------------------------------------------!
        1034             :
        1035             :   use wv_sat_methods, only: wv_sat_qsat_water_vect
        1036             :
        1037             :   ! Inputs
        1038             :   integer,                        intent(in) :: dim1, dim2
        1039             :   real(r8), dimension(dim1,dim2), intent(in) :: t    ! Temperature
        1040             :   real(r8), dimension(dim1,dim2), intent(in) :: p    ! Pressure
        1041             :   ! Outputs
        1042             :   real(r8), dimension(dim1,dim2), intent(out) :: es  ! Saturation vapor pressure
        1043             :   real(r8), dimension(dim1,dim2), intent(out) :: qs  ! Saturation specific humidity
        1044             :
        1045             :   real(r8), dimension(dim1,dim2), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
        1046             :   real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt  ! (d(qs)/dt)
        1047             :   real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
        1048             :
        1049             :   ! Local variables
        1050           0 :   real(r8), dimension(dim1,dim2) :: hltalt       ! Modified latent heat for T derivatives
        1051           0 :   real(r8), dimension(dim1,dim2) :: tterm
        1052             :   integer                        :: i, k, vlen
        1053             :   logical                        :: present_gam, present_dqsdt, present_enthalpy
        1054             :
        1055           0 :   vlen = dim1 * dim2
        1056           0 :   present_gam = present(gam)
        1057           0 :   present_dqsdt = present(dqsdt)
        1058           0 :   present_enthalpy = present(enthalpy)
        1059             :
        1060             :   !$acc data copyin  (t,p) &
        1061             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
        1062             :   !$acc      create  (hltalt,tterm)
        1063             :
        1064             :   !$acc parallel vector_length(VLENS) default(present)
        1065             :   !$acc loop gang vector collapse(2)
        1066           0 :   do k = 1, dim2
        1067           0 :      do i = 1, dim1
        1068           0 :         tterm(i,k) = 0._r8
        1069             :      end do
        1070             :   end do
        1071             :   !$acc end parallel
        1072             :
        1073           0 :   call wv_sat_qsat_water_vect(t, p, es, qs, vlen)
        1074             :
        1075           0 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
        1076             :
        1077             :      ! "generalized" analytic expression for t derivative of es
        1078             :      ! accurate to within 1 percent for 173.16 < t < 373.16
        1079           0 :      call no_ip_hltalt_vect(t, hltalt, vlen)
        1080             :
        1081           0 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
        1082             :
        1083             :      ! For pure water/ice transition term is 0.
        1084             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
        1085           0 :           gam=gam, dqsdt=dqsdt)
        1086             :
        1087             :   end if
        1088             :
        1089             :   !$acc end data
        1090           0 : end subroutine qsat_water_2D
        1091             :
        1092           0 : subroutine qsat_ice_line(t, p, es, qs, gam, dqsdt, enthalpy)
        1093             :   !------------------------------------------------------------------!
        1094             :   ! Purpose:                                                         !
        1095             :   !   Calculate SVP over ice at a given temperature, and then        !
        1096             :   !   calculate and return saturation specific humidity.             !
        1097             :   !   Optionally return various temperature derivatives or enthalpy  !
        1098             :   !   at saturation.                                                 !
        1099             :   !------------------------------------------------------------------!
        1100             :
        1101             :   use wv_sat_methods, only: wv_sat_qsat_ice
        1102             :
        1103             :   ! Inputs
        1104             :   real(r8), intent(in) :: t    ! Temperature
        1105             :   real(r8), intent(in) :: p    ! Pressure
        1106             :   ! Outputs
        1107             :   real(r8), intent(out) :: es  ! Saturation vapor pressure
        1108             :   real(r8), intent(out) :: qs  ! Saturation specific humidity
        1109             :
        1110             :   real(r8), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
        1111             :   real(r8), intent(out), optional :: dqsdt  ! (d(qs)/dt)
        1112             :   real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
        1113             :
        1114             :   ! Local variables
        1115             :   real(r8) :: hltalt       ! Modified latent heat for T derivatives
        1116             :
        1117           0 :   call wv_sat_qsat_ice(t, p, es, qs)
        1118             :
        1119           0 :   if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then
        1120             :
        1121             :      ! For pure ice, just add latent heats.
        1122           0 :      hltalt = latvap + latice
        1123             :
        1124           0 :      if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt)
        1125             :
        1126             :      ! For pure water/ice transition term is 0.
        1127             :      call deriv_outputs_line(t, p, es, qs, hltalt, 0._r8, &
        1128           0 :           gam=gam, dqsdt=dqsdt)
        1129             :
        1130             :   end if
        1131             :
        1132           0 : end subroutine qsat_ice_line
        1133             :
        1134           0 : subroutine qsat_ice_vect(t, p, es, qs, vlen, gam, dqsdt, enthalpy)
        1135             :   !------------------------------------------------------------------!
        1136             :   ! Purpose:                                                         !
        1137             :   !   Calculate SVP over ice at a given temperature, and then        !
        1138             :   !   calculate and return saturation specific humidity.             !
        1139             :   !   Optionally return various temperature derivatives or enthalpy  !
        1140             :   !   at saturation.                                                 !
        1141             :   !------------------------------------------------------------------!
        1142             :
        1143             :   use wv_sat_methods, only: wv_sat_qsat_ice_vect
        1144             :
        1145             :   ! Inputs
        1146             :   integer,                   intent(in) :: vlen
        1147             :   real(r8), dimension(vlen), intent(in) :: t    ! Temperature
        1148             :   real(r8), dimension(vlen), intent(in) :: p    ! Pressure
        1149             :   ! Outputs
        1150             :   real(r8), dimension(vlen), intent(out) :: es  ! Saturation vapor pressure
        1151             :   real(r8), dimension(vlen), intent(out) :: qs  ! Saturation specific humidity
        1152             :
        1153             :   real(r8), dimension(vlen), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
        1154             :   real(r8), dimension(vlen), intent(out), optional :: dqsdt  ! (d(qs)/dt)
        1155             :   real(r8), dimension(vlen), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
        1156             :
        1157             :   ! Local variables
        1158           0 :   real(r8), dimension(vlen) :: hltalt       ! Modified latent heat for T derivatives
        1159           0 :   real(r8), dimension(vlen) :: tterm
        1160             :   integer                   :: i
        1161             :   logical                   :: present_gam, present_dqsdt, present_enthalpy
        1162             :
        1163           0 :   present_gam = present(gam)
        1164           0 :   present_dqsdt = present(dqsdt)
        1165           0 :   present_enthalpy = present(enthalpy)
        1166             :
        1167             :   !$acc data copyin  (t,p) &
        1168             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
        1169             :   !$acc      create  (hltalt,tterm)
        1170             :
        1171             :   !$acc parallel vector_length(VLENS) default(present)
        1172             :   !$acc loop gang vector
        1173           0 :   do i = 1, vlen
        1174           0 :      tterm(i) = 0._r8
        1175             :   end do
        1176             :   !$acc end parallel
        1177             :
        1178           0 :   call wv_sat_qsat_ice_vect(t, p, es, qs, vlen)
        1179             :
        1180           0 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
        1181             :
        1182             :      !$acc parallel vector_length(VLENS) default(present)
        1183             :      !$acc loop gang vector
        1184           0 :      do i = 1, vlen
        1185             :         ! For pure ice, just add latent heats.
        1186           0 :         hltalt(i) = latvap + latice
        1187             :      end do
        1188             :      !$acc end parallel
        1189             :
        1190           0 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
        1191             :
        1192             :      ! For pure water/ice transition term is 0.
        1193             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
        1194           0 :           gam=gam, dqsdt=dqsdt)
        1195             :
        1196             :   end if
        1197             :
        1198             :   !$acc end data
        1199           0 : end subroutine qsat_ice_vect
        1200             :
        1201           0 : subroutine qsat_ice_2D(t, p, es, qs, dim1, dim2, gam, dqsdt, enthalpy)
        1202             :   !------------------------------------------------------------------!
        1203             :   ! Purpose:                                                         !
        1204             :   !   Calculate SVP over ice at a given temperature, and then        !
        1205             :   !   calculate and return saturation specific humidity.             !
        1206             :   !   Optionally return various temperature derivatives or enthalpy  !
        1207             :   !   at saturation.                                                 !
        1208             :   !------------------------------------------------------------------!
        1209             :
        1210             :   use wv_sat_methods, only: wv_sat_qsat_ice_vect
        1211             :
        1212             :   ! Inputs
        1213             :   integer,                        intent(in) :: dim1, dim2
        1214             :   real(r8), dimension(dim1,dim2), intent(in) :: t    ! Temperature
        1215             :   real(r8), dimension(dim1,dim2), intent(in) :: p    ! Pressure
        1216             :   ! Outputs
        1217             :   real(r8), dimension(dim1,dim2), intent(out) :: es  ! Saturation vapor pressure
        1218             :   real(r8), dimension(dim1,dim2), intent(out) :: qs  ! Saturation specific humidity
        1219             :
        1220             :   real(r8), dimension(dim1,dim2), intent(out), optional :: gam    ! (l/cpair)*(d(qs)/dt)
        1221             :   real(r8), dimension(dim1,dim2), intent(out), optional :: dqsdt  ! (d(qs)/dt)
        1222             :   real(r8), dimension(dim1,dim2), intent(out), optional :: enthalpy ! cpair*t + hltalt*q
        1223             :
        1224             :   ! Local variables
        1225           0 :   real(r8), dimension(dim1,dim2) :: hltalt       ! Modified latent heat for T derivatives
        1226           0 :   real(r8), dimension(dim1,dim2) :: tterm
        1227             :   integer                        :: i, k, vlen
        1228             :   logical                        :: present_gam, present_dqsdt, present_enthalpy
        1229             :
        1230           0 :   vlen = dim1 * dim2
        1231           0 :   present_gam = present(gam)
        1232           0 :   present_dqsdt = present(dqsdt)
        1233           0 :   present_enthalpy = present(enthalpy)
        1234             :
        1235             :   !$acc data copyin  (t,p) &
        1236             :   !$acc      copyout (es,qs,gam,dqsdt,enthalpy) &
        1237             :   !$acc      create  (hltalt,tterm)
        1238             :
        1239             :   !$acc parallel vector_length(VLENS) default(present)
        1240             :   !$acc loop gang vector collapse(2)
        1241           0 :   do k = 1, dim2
        1242           0 :      do i = 1, dim1
        1243           0 :         tterm(i,k) = 0._r8
        1244             :      end do
        1245             :   end do
        1246             :   !$acc end parallel
        1247             :
        1248           0 :   call wv_sat_qsat_ice_vect(t, p, es, qs, vlen)
        1249             :
        1250           0 :   if (present_gam .or. present_dqsdt .or. present_enthalpy) then
        1251             :
        1252             :      !$acc parallel vector_length(VLENS) default(present)
        1253             :      !$acc loop gang vector collapse(2)
        1254           0 :      do k = 1, dim2
        1255           0 :         do i = 1, dim1
        1256             :            ! For pure ice, just add latent heats.
        1257           0 :            hltalt(i,k) = latvap + latice
        1258             :         end do
        1259             :      end do
        1260             :      !$acc end parallel
        1261             :
        1262           0 :      if (present_enthalpy) call tq_enthalpy_vect(t, qs, hltalt, enthalpy, vlen)
        1263             :
        1264             :      ! For pure water/ice transition term is 0.
        1265             :      call deriv_outputs_vect(t, p, es, qs, hltalt, tterm, vlen, &
        1266           0 :           gam=gam, dqsdt=dqsdt)
        1267             :
        1268             :   end if
        1269             :
        1270             :   !$acc end data
        1271           0 : end subroutine qsat_ice_2D
        1272             :
        1273             : !---------------------------------------------------------------------
        1274             : ! FINDSP (WET BULB TEMPERATURE) PROCEDURES
        1275             : !---------------------------------------------------------------------
        1276             :
        1277    38879568 : subroutine findsp_vc(q, t, p, use_ice, tsp, qsp)
        1278             :
        1279             :   use cam_logfile,      only: iulog
        1280             :   use cam_abortutils,   only: endrun
        1281             :
        1282             :   ! Wrapper for findsp which is 1D and handles the output status.
        1283             :   ! Changing findsp to elemental restricted debugging output.
        1284             :   ! If that output is needed again, it's preferable *not* to copy findsp,
        1285             :   ! but to change the existing version.
        1286             :
        1287             :   ! input arguments
        1288             :   real(r8), intent(in) :: q(:)        ! water vapor (kg/kg)
        1289             :   real(r8), intent(in) :: t(:)        ! temperature (K)
        1290             :   real(r8), intent(in) :: p(:)        ! pressure    (Pa)
        1291             :   logical,  intent(in) :: use_ice     ! flag to include ice phase in calculations
        1292             :
        1293             :   ! output arguments
        1294             :   real(r8), intent(out) :: tsp(:)     ! saturation temp (K)
        1295             :   real(r8), intent(out) :: qsp(:)     ! saturation mixing ratio (kg/kg)
        1296             :
        1297    77759136 :   integer :: status(size(q))   ! flag representing state of output
        1298             :                                ! 0 => Successful convergence
        1299             :                                ! 1 => No calculation done: pressure or specific
        1300             :                                !      humidity not within usable range
        1301             :                                ! 2 => Run failed to converge
        1302             :                                ! 4 => Temperature fell below minimum
        1303             :                                ! 8 => Enthalpy not conserved
        1304             :
        1305             :   integer :: n, i
        1306             :
        1307    38879568 :   n = size(q)
        1308             :
        1309             :   ! Currently, only 2 and 8 seem to be treated as fatal errors.
        1310   649198368 :   do i = 1,n
        1311   610318800 :      call findsp(q(i), t(i), p(i), use_ice, tsp(i), qsp(i), status(i))
        1312   649198368 :      if (status(i) == 2) then
        1313           0 :         write(iulog,*) ' findsp not converging at i = ', i
        1314           0 :         write(iulog,*) ' t, q, p ', t(i), q(i), p(i)
        1315           0 :         write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i)
        1316           0 :         call endrun ('wv_saturation::FINDSP -- not converging')
        1317   610318800 :      else if (status(i) == 8) then
        1318           0 :         write(iulog,*) ' the enthalpy is not conserved at i = ', i
        1319           0 :         write(iulog,*) ' t, q, p ', t(i), q(i), p(i)
        1320           0 :         write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i)
        1321           0 :         call endrun ('wv_saturation::FINDSP -- enthalpy is not conserved')
        1322             :      endif
        1323             :   end do
        1324             :
        1325    38879568 : end subroutine findsp_vc
        1326             :
        1327   610318800 : subroutine findsp (q, t, p, use_ice, tsp, qsp, status)
        1328             : !-----------------------------------------------------------------------
        1329             : !
        1330             : ! Purpose:
        1331             : !     find the wet bulb temperature for a given t and q
        1332             : !     in a longitude height section
        1333             : !     wet bulb temp is the temperature and spec humidity that is
        1334             : !     just saturated and has the same enthalpy
        1335             : !     if q > qs(t) then tsp > t and qsp = qs(tsp) < q
        1336             : !     if q < qs(t) then tsp < t and qsp = qs(tsp) > q
        1337             : !
        1338             : ! Method:
        1339             : ! a Newton method is used
        1340             : ! first guess uses an algorithm provided by John Petch from the UKMO
        1341             : ! we exclude points where the physical situation is unrealistic
        1342             : ! e.g. where the temperature is outside the range of validity for the
        1343             : !      saturation vapor pressure, or where the water vapor pressure
        1344             : !      exceeds the ambient pressure, or the saturation specific humidity is
        1345             : !      unrealistic
        1346             : !
        1347             : ! Author: P. Rasch
        1348             : !
        1349             : !-----------------------------------------------------------------------
        1350             : !
        1351             : !     input arguments
        1352             : !
        1353             :
        1354             :   real(r8), intent(in) :: q        ! water vapor (kg/kg)
        1355             :   real(r8), intent(in) :: t        ! temperature (K)
        1356             :   real(r8), intent(in) :: p        ! pressure    (Pa)
        1357             :   logical,  intent(in) :: use_ice  ! flag to include ice phase in calculations
        1358             : !
        1359             : ! output arguments
        1360             : !
        1361             :   real(r8), intent(out) :: tsp      ! saturation temp (K)
        1362             :   real(r8), intent(out) :: qsp      ! saturation mixing ratio (kg/kg)
        1363             :   integer,  intent(out) :: status   ! flag representing state of output
        1364             :                                     ! 0 => Successful convergence
        1365             :                                     ! 1 => No calculation done: pressure or specific
        1366             :                                     !      humidity not within usable range
        1367             :                                     ! 2 => Run failed to converge
        1368             :                                     ! 4 => Temperature fell below minimum
        1369             :                                     ! 8 => Enthalpy not conserved
        1370             : !
        1371             : ! local variables
        1372             : !
        1373             :   integer, parameter :: iter = 8    ! max number of times to iterate the calculation
        1374             :   integer :: l                      ! iterator
        1375             :
        1376             :   real(r8) es                   ! sat. vapor pressure
        1377             :   real(r8) gam                  ! change in sat spec. hum. wrt temperature (times hltalt/cpair)
        1378             :   real(r8) dgdt                 ! work variable
        1379             :   real(r8) g                    ! work variable
        1380             :   real(r8) hltalt               ! lat. heat. of vap.
        1381             :   real(r8) qs                   ! spec. hum. of water vapor
        1382             :
        1383             : ! work variables
        1384             :   real(r8) t1, q1, dt, dq
        1385             :   real(r8) qvd
        1386             :   real(r8) r1b, c1, c2
        1387             :   real(r8), parameter :: dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
        1388             :   real(r8), parameter :: dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
        1389             :   real(r8) enin, enout
        1390             :
        1391             :   ! Saturation specific humidity at this temperature
        1392   610318800 :   if (use_ice) then
        1393   610318800 :      call qsat(t, p, es, qs)
        1394             :   else
        1395           0 :      call qsat_water(t, p, es, qs)
        1396             :   end if
        1397             :
        1398             :   ! make sure a meaningful calculation is possible
        1399             :   if (p <= 5._r8*es .or. qs <= 0._r8 .or. qs >= 0.5_r8 &
        1400   610318800 :        .or. t < tmin .or. t > tmax) then
        1401       40322 :      status = 1
        1402             :      ! Keep initial parameters when conditions aren't suitable
        1403       40322 :      tsp = t
        1404       40322 :      qsp = q
        1405       40322 :      enin = 1._r8
        1406       40322 :      enout = 1._r8
        1407             :
        1408       40322 :      return
        1409             :   end if
        1410             :
        1411             :   ! Prepare to iterate
        1412   610278478 :   status = 2
        1413             :
        1414             :   ! Get initial enthalpy
        1415   610278478 :   if (use_ice) then
        1416   610278478 :      call calc_hltalt(t,hltalt)
        1417             :   else
        1418           0 :      call no_ip_hltalt(t,hltalt)
        1419             :   end if
        1420   610278478 :   enin = tq_enthalpy(t, q, hltalt)
        1421             :
        1422             :   ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
        1423   610278478 :   c1 = hltalt*c3
        1424   610278478 :   c2 = (t + 36._r8)**2
        1425   610278478 :   r1b = c2/(c2 + c1*qs)
        1426   610278478 :   qvd = r1b * (q - qs)
        1427   610278478 :   tsp = t + ((hltalt/cpair)*qvd)
        1428             :
        1429             :   ! Generate qsp, gam, and enout from tsp.
        1430   610278478 :   if (use_ice) then
        1431   610278478 :      call qsat(tsp, p, es, qsp, gam=gam, enthalpy=enout)
        1432             :   else
        1433           0 :      call qsat_water(tsp, p, es, qsp, gam=gam, enthalpy=enout)
        1434             :   end if
        1435             :
        1436             :   ! iterate on first guess
        1437  1251599190 :   do l = 1, iter
        1438             :
        1439  1251599190 :      g = enin - enout
        1440  1251599190 :      dgdt = -cpair * (1 + gam)
        1441             :
        1442             :      ! New tsp
        1443  1251599190 :      t1 = tsp - g/dgdt
        1444  1251599190 :      dt = abs(t1 - tsp)/t1
        1445  1251599190 :      tsp = t1
        1446             :
        1447             :      ! bail out if past end of temperature range
        1448  1251599190 :      if ( tsp < tmin ) then
        1449           0 :         tsp = tmin
        1450             :         ! Get latent heat and set qsp to a value
        1451             :         ! that preserves enthalpy.
        1452           0 :         if (use_ice) then
        1453           0 :            call calc_hltalt(tsp,hltalt)
        1454             :         else
        1455           0 :            call no_ip_hltalt(tsp,hltalt)
        1456             :         end if
        1457           0 :         qsp = (enin - cpair*tsp)/hltalt
        1458           0 :         enout = tq_enthalpy(tsp, qsp, hltalt)
        1459           0 :         status = 4
        1460           0 :         exit
        1461             :      end if
        1462             :
        1463             :      ! Re-generate qsp, gam, and enout from new tsp.
        1464  1251599190 :      if (use_ice) then
        1465  1251599190 :         call qsat(tsp, p, es, q1, gam=gam, enthalpy=enout)
        1466             :      else
        1467           0 :         call qsat_water(tsp, p, es, q1, gam=gam, enthalpy=enout)
        1468             :      end if
        1469  1251599190 :      dq = abs(q1 - qsp)/max(q1,1.e-12_r8)
        1470  1251599190 :      qsp = q1
        1471             :
        1472             :      ! if converged at this point, exclude it from more iterations
        1473  1251599190 :      if (dt < dttol .and. dq < dqtol) then
        1474   610278478 :         status = 0
        1475   610278478 :         exit
        1476             :      endif
        1477             :   end do
        1478             :
        1479             :   ! Test for enthalpy conservation
        1480   610278478 :   if (abs((enin-enout)/(enin+enout)) > 1.e-4_r8) status = 8
        1481             :
        1482             : end subroutine findsp
        1483             :
        1484             : end module wv_saturation

Generated by: LCOV version 1.14