LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - saturation.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 37 | 116 | 31.9 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 10 | 30.0 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------
           2             : !$Id$
           3             : !===============================================================================
           4             : module saturation
           5             :
           6             : ! Description:
           7             : !   Contains functions that compute saturation with respect
           8             : !   to liquid or ice.
           9             : !-----------------------------------------------------------------------
          10             :
          11             :   use model_flags, only: &  ! h1g, 2010-06-18
          12             :      I_sat_sphum
          13             :
          14             :   use clubb_precision, only: &
          15             :     core_rknd ! Variable(s)
          16             :
          17             :   use model_flags, only: &
          18             :     saturation_bolton, &
          19             :     saturation_gfdl, &
          20             :     saturation_flatau, &
          21             :     saturation_lookup
          22             :
          23             :   implicit none
          24             :
          25             :   private  ! Change default so all items private
          26             :
          27             :   public   :: sat_mixrat_liq, sat_mixrat_ice, rcm_sat_adj, sat_vapor_press_liq
          28             :
          29             :   private  :: sat_vapor_press_ice_flatau, sat_vapor_press_ice_bolton
          30             :
          31             :   interface sat_mixrat_liq
          32             :     module procedure sat_mixrat_liq_k   ! Works over a single vertical level
          33             :     module procedure sat_mixrat_liq_2D  ! Works over all vertical levels and columns
          34             :   end interface sat_mixrat_liq
          35             :
          36             :   interface sat_mixrat_ice
          37             :     module procedure sat_mixrat_ice_k   ! Works over a single vertical level
          38             :     module procedure sat_mixrat_ice_2D  ! Works over all vertical levels and columns
          39             :   end interface sat_mixrat_ice
          40             :
          41             :   ! Lookup table of values for saturation
          42             :   real( kind = core_rknd ), parameter, private, dimension(188:343) :: &
          43             :     svp_liq_lookup_table = &
          44             :     (/ 0.049560547_core_rknd, 0.059753418_core_rknd, 0.070129395_core_rknd, 0.083618164_core_rknd, &
          45             :     0.09814453_core_rknd, 0.11444092_core_rknd, 0.13446045_core_rknd, 0.15686035_core_rknd,     &
          46             :     0.18218994_core_rknd, 0.21240234_core_rknd, 0.24725342_core_rknd, 0.28668213_core_rknd,     &
          47             :     0.33184814_core_rknd, 0.3826294_core_rknd, 0.4416504_core_rknd, 0.50775146_core_rknd,       &
          48             :     0.58343506_core_rknd, 0.6694946_core_rknd, 0.7668457_core_rknd, 0.87750244_core_rknd,       &
          49             :     1.0023804_core_rknd, 1.1434937_core_rknd, 1.3028564_core_rknd, 1.482544_core_rknd,          &
          50             :     1.6847534_core_rknd, 1.9118042_core_rknd, 2.1671143_core_rknd, 2.4535522_core_rknd,         &
          51             :     2.774231_core_rknd, 3.1330566_core_rknd, 3.5343628_core_rknd, 3.9819336_core_rknd,          &
          52             :     4.480713_core_rknd, 5.036072_core_rknd, 5.6540527_core_rknd, 6.340088_core_rknd,            &
          53             :     7.1015015_core_rknd, 7.9450684_core_rknd, 8.8793335_core_rknd, 9.91217_core_rknd,           &
          54             :     11.053528_core_rknd, 12.313049_core_rknd, 13.70166_core_rknd, 15.231018_core_rknd,          &
          55             :     16.91394_core_rknd, 18.764038_core_rknd, 20.795898_core_rknd, 23.025574_core_rknd,          &
          56             :     25.470093_core_rknd, 28.147766_core_rknd, 31.078003_core_rknd, 34.282043_core_rknd,         &
          57             :     37.782593_core_rknd, 41.60382_core_rknd, 45.771606_core_rknd, 50.31366_core_rknd,           &
          58             :     55.259644_core_rknd, 60.641174_core_rknd, 66.492004_core_rknd, 72.84802_core_rknd,          &
          59             :     79.74756_core_rknd, 87.23126_core_rknd, 95.34259_core_rknd, 104.12747_core_rknd,            &
          60             :     113.634796_core_rknd, 123.91641_core_rknd, 135.02725_core_rknd, 147.02563_core_rknd,        &
          61             :     159.97308_core_rknd, 173.93488_core_rknd, 188.97995_core_rknd, 205.18109_core_rknd,         &
          62             :     222.61517_core_rknd, 241.36334_core_rknd, 261.51108_core_rknd, 283.14853_core_rknd,         &
          63             :     306.37054_core_rknd, 331.27698_core_rknd, 357.97278_core_rknd, 386.56842_core_rknd,         &
          64             :     417.17978_core_rknd, 449.9286_core_rknd, 484.94254_core_rknd, 522.3556_core_rknd,           &
          65             :     562.30804_core_rknd, 604.947_core_rknd, 650.42645_core_rknd, 698.9074_core_rknd,            &
          66             :     750.55835_core_rknd, 805.55554_core_rknd, 864.0828_core_rknd, 926.3325_core_rknd,           &
          67             :     992.5052_core_rknd, 1062.8102_core_rknd, 1137.4657_core_rknd, 1216.6995_core_rknd,          &
          68             :     1300.7483_core_rknd, 1389.8594_core_rknd, 1484.2896_core_rknd, 1584.3064_core_rknd,         &
          69             :     1690.1881_core_rknd, 1802.224_core_rknd, 1920.7146_core_rknd, 2045.9724_core_rknd,          &
          70             :     2178.3218_core_rknd, 2318.099_core_rknd, 2465.654_core_rknd, 2621.3489_core_rknd,           &
          71             :     2785.5596_core_rknd, 2958.6758_core_rknd, 3141.101_core_rknd, 3333.2534_core_rknd,          &
          72             :     3535.5657_core_rknd, 3748.4863_core_rknd, 3972.4792_core_rknd, 4208.024_core_rknd,          &
          73             :     4455.616_core_rknd, 4715.7686_core_rknd, 4989.0127_core_rknd, 5275.8945_core_rknd,          &
          74             :     5576.9795_core_rknd, 5892.8535_core_rknd, 6224.116_core_rknd, 6571.3926_core_rknd,          &
          75             :     6935.3213_core_rknd, 7316.5674_core_rknd, 7715.8105_core_rknd, 8133.755_core_rknd,          &
          76             :     8571.125_core_rknd, 9028.667_core_rknd, 9507.15_core_rknd, 10007.367_core_rknd,             &
          77             :     10530.132_core_rknd, 11076.282_core_rknd, 11646.683_core_rknd, 12242.221_core_rknd,         &
          78             :     12863.808_core_rknd, 13512.384_core_rknd, 14188.913_core_rknd, 14894.385_core_rknd,         &
          79             :     15629.823_core_rknd, 16396.268_core_rknd, 17194.799_core_rknd, 18026.516_core_rknd,         &
          80             :     18892.55_core_rknd, 19794.07_core_rknd, 20732.262_core_rknd, 21708.352_core_rknd,           &
          81             :     22723.592_core_rknd, 23779.273_core_rknd, 24876.709_core_rknd, 26017.258_core_rknd,         &
          82             :     27202.3_core_rknd, 28433.256_core_rknd, 29711.578_core_rknd, 31038.766_core_rknd /)
          83             : !$acc declare create( svp_liq_lookup_table )
          84             :
          85             :
          86             :   contains
          87             :
          88             :   !-------------------------------------------------------------------------
          89             :   ! Wrapped in interface sat_mixrat_liq
          90  6368645631 :   function sat_mixrat_liq_k( p_in_Pa, T_in_K, &
          91             :                              saturation_formula )
          92             : !$acc routine seq
          93             :   ! Description:
          94             :   !   Used to compute the saturation mixing ratio of liquid water.
          95             :
          96             :   ! References:
          97             :   !   Formula from Emanuel 1994, 4.4.14
          98             :   !-------------------------------------------------------------------------
          99             :
         100             :     use constants_clubb, only: &
         101             :         ep
         102             :
         103             :     use constants_clubb, only: &
         104             :         T_freeze_K
         105             :
         106             :     implicit none
         107             :
         108             :     ! -------------------- Input Variables --------------------
         109             :     real( kind = core_rknd ), intent(in) ::  &
         110             :       p_in_Pa,  & ! Pressure    [Pa]
         111             :       T_in_K      ! Temperature [K]
         112             :
         113             :     integer, intent(in) :: &
         114             :       saturation_formula
         115             :
         116             :     ! -------------------- Output Variables --------------------
         117             :     real( kind = core_rknd ) ::  &
         118             :       sat_mixrat_liq_k
         119             :
         120             :     ! -------------------- Local Variables --------------------
         121             :     real( kind = core_rknd ) :: &
         122             :         T_in_C, &
         123             :         T_in_C_sqd, &
         124             :         T_in_K_clipped
         125             :
         126             :       integer :: &
         127             :         T_in_K_int
         128             :
         129             :     ! Constant parameters
         130             :
         131             :     ! Relative error norm expansion (-50 to 50 deg_C) from
         132             :     ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
         133             :     ! (The 100 coefficient converts from mb to Pa)
         134             :     !   real, dimension(7), parameter :: a = &
         135             :     !   100.* (/ 6.11176750,      0.443986062,     0.143053301E-01, &
         136             :     !            0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
         137             :     !            0.638780966E-10 /)
         138             :
         139             :     ! Relative error norm expansion (-85 to 70 deg_C) from
         140             :     ! Table 4 of pp. 1511 of Flatau et al.
         141             :     !real( kind = core_rknd ), dimension(9), parameter :: a = &
         142             :     !100._core_rknd * &
         143             :     !  Commented out because the form has been redone, causing these number to no longer be needed,
         144             :     !  leaving them in for now for reference.
         145             :     !         (/ 6.11583699_core_rknd,      0.444606896_core_rknd,     0.143177157E-01_core_rknd, &
         146             :     !         0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
         147             :     !         0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
         148             :
         149             :     real( kind = core_rknd ), parameter :: &
         150             :       min_T_in_C = -85._core_rknd,  & ! [deg_C]
         151             :       min_T_in_K = 173.15_core_rknd   ! Lowest temperature at which Goff-Gratch is valid [K]
         152             :
         153             :     ! ---------------------- Output Variables ----------------------
         154             :     real( kind = core_rknd ) :: &
         155             :       esat  ! Saturation vapor pressure over water [Pa]
         156             :
         157             :     ! -------------------- Begin Code --------------------
         158             :
         159             :     ! Calculate the SVP for water vapor.
         160  6368645631 :     select case ( saturation_formula )
         161             :     case ( saturation_flatau )
         162             :
         163             :       ! Using the Flatau, et al. polynomial approximation for SVP over vapor
         164             :
         165             :       ! Determine deg K - 273.15
         166           0 :       T_in_C = T_in_K - T_freeze_K
         167             :
         168             :       ! Since this approximation is only good out to -85 degrees Celsius we
         169             :       ! truncate the result here (Flatau, et al. 1992)
         170           0 :       T_in_C = max( T_in_C, min_T_in_C )
         171             :
         172             :       ! Polynomial approx. (Flatau, et al. 1992)
         173             :
         174             :       ! This is the generalized formula but is not computationally efficient.
         175             :       ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
         176             :       ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
         177             :
         178             :       ! esat = a(1)
         179             :
         180             :       ! do i = 2, size( a ) , 1
         181             :       !   esat = esat + a(i) * ( T_in_C )**(i-1)
         182             :       ! end do
         183             :
         184             :       ! The 8th order polynomial fit.  When running deep
         185             :       ! convective cases I noticed that absolute temperature often dips below
         186             :       ! -50 deg_C at higher altitudes, where the 6th order approximation is
         187             :       ! not accurate.  -dschanen 20 Nov 2008
         188             :       !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
         189             :       !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
         190             :
         191             :
         192             :       ! Factoring the polynomial above and changing it into this form allows the cpu
         193             :       ! to complete the calculations out of order. This is because modern cpus can complete
         194             :       ! multiple instructions at once if they do not depend on eachother, in the above case
         195             :       ! each instruction relies on the result of the last. In this version however, the terms
         196             :       ! in the parentheses could potentially be calculated in parallel by different execution
         197             :       ! units in the cpu, then only when those terms are being multiplied together do the
         198             :       ! instructions need to be done one at a time. See clubb issue 834 for more info.
         199             :       !   - Gunther Huebler, Aug 2018
         200           0 :       T_in_C_sqd = T_in_C**2
         201             :
         202             :       esat = &
         203             :         - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
         204             :           * ( T_in_C + 90.72381630364440_core_rknd ) &
         205             :           * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
         206             :           * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
         207           0 :           * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
         208             :
         209             :     case ( saturation_bolton )
         210             :
         211             :       ! Using the Bolton 1980 approximations for SVP over vapor
         212             :       ! Generally this more computationally expensive than the Flatau polnomial expansion
         213             :       esat = 611.2_core_rknd &
         214             :               * exp( (17.67_core_rknd *(T_in_K-T_freeze_K))  &
         215           0 :                      / (T_in_K-29.65_core_rknd) ) ! Known magic number
         216             :
         217             : ! ---> h1g
         218             :     case ( saturation_gfdl )
         219             :
         220             :       ! Using GFDL polynomial approximation for SVP with respect to liquid
         221             :       ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
         222             :       !   we threshold the temperature.  This will yield a minimal saturation at
         223             :       !   cold temperatures.
         224  6368645631 :       T_in_K_clipped = max( min_T_in_K, T_in_K )
         225             :
         226             :       ! Goff Gratch equation, uncertain below -70 C
         227             :
         228             :       esat = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
         229             :            5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
         230             :            1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
         231             :              *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
         232             :            8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
         233             :              *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
         234  6368645631 :            log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
         235             :
         236             : ! <--- h1g
         237             :
         238             :     case ( saturation_lookup )
         239             :
         240           0 :       T_in_K_int = int( anint( T_in_K ) )
         241             :
         242             :       ! Since this approximation is only good out to -85 degrees Celsius we
         243             :       ! truncate the result here
         244           0 :       T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
         245             :
         246             :       ! Use the lookup table to determine the saturation vapor pressure.
         247           0 :       esat = svp_liq_lookup_table( T_in_K_int )
         248             :
         249             :     case default
         250             :
         251             :       ! Undefined approximation
         252  6368645631 :       esat = -99999.999_core_rknd
         253             :
         254             :     end select
         255             :
         256             :     ! If esat exceeds the air pressure, then assume esat~=0.5*pressure
         257             :     !   and set rsat = ep = 0.622
         258  6368645631 :     if ( p_in_Pa-esat < 1.0_core_rknd ) then
         259             :       sat_mixrat_liq_k = ep
         260             :     else
         261             :
         262  6368645631 :       if ( I_sat_sphum )  then   ! h1g, 2010-06-18 begin mod
         263             :
         264             :         ! GFDL uses specific humidity
         265             :         ! Formula for Saturation Specific Humidity
         266             :         sat_mixrat_liq_k = ep * ( esat / ( p_in_Pa &
         267             :                                              - (1.0_core_rknd-ep) * esat ) )
         268             :       else
         269             :
         270             :         ! Formula for Saturation Mixing Ratio:
         271             :         !
         272             :         ! rs = (epsilon) * [ esat / ( p - esat ) ];
         273             :         ! where epsilon = R_d / R_v
         274  6368645631 :         sat_mixrat_liq_k = ep * esat / ( p_in_Pa - esat )
         275             :
         276             :       endif                     ! h1g, 2010-06-18 end mod
         277             :
         278             :     end if
         279             :
         280             :     return
         281             :
         282             :   end function sat_mixrat_liq_k
         283             :
         284             :   !-------------------------------------------------------------------------
         285             :   !
         286    80415504 :   function sat_mixrat_liq_2D( nz, ngrdcol, p_in_Pa, T_in_K, &
         287             :                               saturation_formula, &
         288             :                               start_index_in )
         289             :   !
         290             :   ! Description:
         291             :   !   Used to compute the saturation mixing ratio of liquid water.
         292             :   !
         293             :   ! References:
         294             :   !   Formula from Emanuel 1994, 4.4.14
         295             :   !-------------------------------------------------------------------------
         296             :
         297             :     use constants_clubb, only: &
         298             :       ep
         299             :
         300             :     use constants_clubb, only: T_freeze_K
         301             :
         302             :     implicit none
         303             :
         304             :     ! -------------------- Input Variables --------------------
         305             :     integer, intent(in) :: &
         306             :       nz, &
         307             :       ngrdcol
         308             :
         309             :     integer, intent(in) :: &
         310             :       saturation_formula
         311             :
         312             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         313             :       p_in_Pa,  & ! Pressure    [Pa]
         314             :       T_in_K      ! Temperature [K]
         315             :
         316             :     ! -------------------- Output Variables --------------------
         317             :     real( kind = core_rknd ), dimension(ngrdcol,nz) ::  &
         318             :       sat_mixrat_liq_2D
         319             :
         320             :     ! -------------------- Optional Input --------------------
         321             :     integer, intent(in), optional :: &
         322             :       start_index_in
         323             :
         324             :     ! -------------------- Local Variables --------------------
         325             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         326    80415504 :       esat
         327             :
         328             :     real( kind = core_rknd ) :: &
         329             :       T_in_C, &
         330             :       T_in_C_sqd, &
         331             :       T_in_K_clipped
         332             :
         333             :     integer :: &
         334             :       T_in_K_int
         335             :
         336             :     integer :: i,k
         337             :
         338             :     ! Constant parameters
         339             :
         340             :     ! Relative error norm expansion (-50 to 50 deg_C) from
         341             :     ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
         342             :     ! (The 100 coefficient converts from mb to Pa)
         343             :     !   real, dimension(7), parameter :: a = &
         344             :     !   100.* (/ 6.11176750,      0.443986062,     0.143053301E-01, &
         345             :     !            0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
         346             :     !            0.638780966E-10 /)
         347             :
         348             :     ! Relative error norm expansion (-85 to 70 deg_C) from
         349             :     ! Table 4 of pp. 1511 of Flatau et al.
         350             :     !real( kind = core_rknd ), dimension(9), parameter :: a = &
         351             :     !100._core_rknd * &
         352             :     !  Commented out because the form has been redone, causing these number to no longer be needed,
         353             :     !  leaving them in for now for reference.
         354             :     !         (/ 6.11583699_core_rknd,      0.444606896_core_rknd,     0.143177157E-01_core_rknd, &
         355             :     !         0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
         356             :     !         0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
         357             :
         358             :     real( kind = core_rknd ), parameter :: &
         359             :       min_T_in_C = -85._core_rknd,  & ! [deg_C]
         360             :       min_T_in_K = 173.15_core_rknd   ! Lowest temperature at which Goff-Gratch is valid [K]
         361             :
         362             :     integer :: &
         363             :       start_index
         364             :
         365             :     ! -------------------- Begin Code --------------------
         366             :
         367             :     !$acc data create(esat) &
         368             :     !$acc      copyin(p_in_Pa,T_in_K), &
         369             :     !$acc      copyout(sat_mixrat_liq_2D)
         370             :
         371             :     ! start_index is an optional argument and
         372             :     ! used for choosing the sub-arrays
         373    80415504 :     if ( present(start_index_in) ) then
         374    17870112 :       start_index = start_index_in
         375             :     else
         376             :       start_index = 1
         377             :     end if
         378             :
         379    80415504 :     select case ( saturation_formula )
         380             :     case ( saturation_flatau )
         381             :
         382             :       !$acc parallel loop gang vector collapse(2) default(present)
         383           0 :       do k = start_index, nz
         384           0 :         do i = 1, ngrdcol
         385             :
         386             :           ! Determine deg K - 273.15
         387           0 :           T_in_C = T_in_K(i,k) - T_freeze_K
         388             :
         389             :           ! Since this approximation is only good out to -85 degrees Celsius we
         390             :           ! truncate the result here (Flatau, et al. 1992)
         391           0 :           T_in_C = max( T_in_C, min_T_in_C )
         392             :
         393             :           ! Polynomial approx. (Flatau, et al. 1992)
         394             :
         395             :           ! This is the generalized formula but is not computationally efficient.
         396             :           ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
         397             :           ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
         398             :
         399             :           ! esat = a(1)
         400             :
         401             :           ! do i = 2, size( a ) , 1
         402             :           !   esat = esat + a(i) * ( T_in_C )**(i-1)
         403             :           ! end do
         404             :
         405             :           ! The 8th order polynomial fit.  When running deep
         406             :           ! convective cases I noticed that absolute temperature often dips below
         407             :           ! -50 deg_C at higher altitudes, where the 6th order approximation is
         408             :           ! not accurate.  -dschanen 20 Nov 2008
         409             :           !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
         410             :           !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
         411             :
         412             :
         413             :           ! Factoring the polynomial above and changing it into this form allows the cpu
         414             :           ! to complete the calculations out of order. This is because modern cpus can complete
         415             :           ! multiple instructions at once if they do not depend on eachother, in the above case
         416             :           ! each instruction relies on the result of the last. In this version however, the terms
         417             :           ! in the parentheses could potentially be calculated in parallel by different execution
         418             :           ! units in the cpu, then only when those terms are being multiplied together do the
         419             :           ! instructions need to be done one at a time. See clubb issue 834 for more info.
         420             :           !   - Gunther Huebler, Aug 2018
         421           0 :           T_in_C_sqd = T_in_C**2
         422             :
         423             :           esat(i,k) = &
         424             :           - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
         425             :             * ( T_in_C + 90.72381630364440_core_rknd ) &
         426             :             * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
         427             :             * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
         428           0 :             * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
         429             :         end do
         430             :       end do
         431             :       !$acc end parallel loop
         432             :
         433             :     case ( saturation_bolton )
         434             :
         435             :       ! Using the Bolton 1980 approximations for SVP over vapor
         436             :       ! Generally this more computationally expensive than the Flatau polnomial expansion
         437             :       !$acc parallel loop gang vector collapse(2) default(present)
         438           0 :       do k = start_index, nz
         439           0 :         do i = 1, ngrdcol
         440           0 :           esat(i,k) = 611.2_core_rknd &
         441             :                       * exp( (17.67_core_rknd *(T_in_K(i,k)-T_freeze_K))  &
         442           0 :                              / (T_in_K(i,k)-29.65_core_rknd) ) ! Known magic number
         443             :         end do
         444             :       end do
         445             :       !$acc end parallel loop
         446             :
         447             : ! ---> h1g
         448             :     case ( saturation_gfdl )
         449             :
         450             :       ! Using GFDL polynomial approximation for SVP with respect to liquid
         451             :       !$acc parallel loop gang vector collapse(2) default(present)
         452  6888928176 :       do k = start_index, nz
         453 >11376*10^7 :         do i = 1, ngrdcol
         454             :
         455             :           ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
         456             :           !   we threshold the temperature.  This will yield a minimal saturation at
         457             :           !   cold temperatures.
         458 >10687*10^7 :           T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
         459             :
         460             :           ! Goff Gratch equation, uncertain below -70 C
         461             :
         462             :           esat(i,k) = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
         463             :                5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
         464             :                1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
         465             :                  *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
         466             :                8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
         467             :                  *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
         468 >11368*10^7 :                log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
         469             :         end do
         470             :       end do
         471             :       !$acc end parallel loop
         472             :
         473             : ! <--- h1g
         474             :
         475             :     case ( saturation_lookup )
         476             :
         477             :       !$acc parallel loop gang vector collapse(2) default(present)
         478           0 :       do k = start_index, nz
         479           0 :         do i = 1, ngrdcol
         480           0 :           T_in_K_int = int( anint( T_in_K(i,k) ) )
         481             :
         482             :           ! Since this approximation is only good out to -85 degrees Celsius we
         483             :           ! truncate the result here
         484           0 :           T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
         485             :
         486             :           ! Use the lookup table to determine the saturation vapor pressure.
         487           0 :           esat(i,k) = svp_liq_lookup_table( T_in_K_int )
         488             :         end do
         489             :       end do
         490             :       !$acc end parallel loop
         491             :
         492             :     case default
         493             :
         494             :       ! Undefined approximation
         495    80415504 :       esat = -99999.999_core_rknd
         496             :
         497             :     end select
         498             :
         499             :     !$acc parallel loop gang vector collapse(2) default(present)
         500  6888928176 :     do k = start_index, nz
         501 >11376*10^7 :       do i = 1, ngrdcol
         502             :
         503             :         ! If esat exceeds the air pressure, then assume esat~=0.5*pressure
         504             :         !   and set rsat = ep = 0.622
         505 >11368*10^7 :         if ( p_in_Pa(i,k)-esat(i,k) < 1.0_core_rknd ) then
         506  2040744207 :           sat_mixrat_liq_2D(i,k) = ep
         507             :         else
         508             :
         509             :           if ( I_sat_sphum )  then   ! h1g, 2010-06-18 begin mod
         510             :
         511             :             ! GFDL uses specific humidity
         512             :             ! Formula for Saturation Specific Humidity
         513             :             sat_mixrat_liq_2D(i,k) = ep * ( esat(i,k) &
         514             :                                      / ( p_in_Pa(i,k) - (1.0_core_rknd-ep * esat(i,k) ) ) )
         515             :           else
         516             :
         517             :             ! Formula for Saturation Mixing Ratio:
         518             :             !
         519             :             ! rs = (epsilon) * [ esat / ( p - esat ) ];
         520             :             ! where epsilon = R_d / R_v
         521 >10483*10^7 :             sat_mixrat_liq_2D(i,k) = ep * esat(i,k) / ( p_in_Pa(i,k) - esat(i,k) )
         522             :
         523             :           endif                     ! h1g, 2010-06-18 end mod
         524             :
         525             :         end if
         526             :
         527             :       end do
         528             :     end do
         529             :     !$acc end parallel loop
         530             :
         531             :     !$acc end data
         532             :
         533             :   end function sat_mixrat_liq_2D
         534             :
         535             :   !-----------------------------------------------------------------
         536             :   ! Wrapped in interface sat_vapor_press_liq
         537           0 :   subroutine sat_vapor_press_liq( T_in_K, &
         538             :                                   saturation_formula, &
         539             :                                   esat )
         540             :
         541             :   ! Description:
         542             :   !   Computes SVP for water vapor. Calls one of the other functions
         543             :   !   that calculate an approximation to SVP.
         544             :
         545             :   ! References:
         546             :   !   None
         547             :   !-----------------------------------------------------------------
         548             :
         549             :     use clubb_precision, only: &
         550             :         core_rknd ! Variable(s)
         551             :
         552             :     implicit none
         553             :
         554             :     ! ------------------------ Input Variables ------------------------
         555             :     real( kind = core_rknd ), intent(in) :: &
         556             :       T_in_K     ! Temperature                          [K]
         557             :
         558             :     integer, intent(in) :: &
         559             :       saturation_formula
         560             :
         561             :     ! ------------------------ Output Variables ------------------------
         562             :     real( kind = core_rknd ), intent(out) :: &
         563             :       esat      ! Saturation Vapor Pressure over Water [Pa]
         564             :
         565             :     ! ------------------------ Being Code ------------------------
         566             :
         567             :     ! Saturation Vapor Pressure, esat, can be found to be approximated
         568             :     ! in many different ways.
         569           0 :     select case ( saturation_formula )
         570             :     case ( saturation_bolton )
         571             :
         572             :       ! Using the Bolton 1980 approximations for SVP over vapor
         573             :       call sat_vapor_press_liq_bolton( T_in_K, &
         574           0 :                                        esat )
         575             :
         576             :     !Earthworks case
         577             :     case ( saturation_flatau )
         578             :
         579             :       ! Using the Flatau, et al. polynomial approximation for SVP over vapor
         580             :       call sat_vapor_press_liq_flatau( T_in_K, &
         581           0 :                                        esat )
         582             :
         583             : ! ---> h1g
         584             :     case ( saturation_gfdl )
         585             :
         586             :       ! Using GFDL polynomial approximation for SVP with respect to liquid
         587             :       call sat_vapor_press_liq_gfdl( T_in_K, &
         588           0 :                                      esat )
         589             :
         590             : ! <--- h1g
         591             :     case ( saturation_lookup )
         592             :
         593             :       ! Use the lookup table to determine the saturation vapor pressure.
         594           0 :       esat = sat_vapor_press_liq_lookup( T_in_K )
         595             :
         596             :     case default
         597             :
         598             :       ! Undefined approximation
         599           0 :       esat = -99999.999_core_rknd
         600             :
         601             :     end select
         602             :
         603           0 :     return
         604             :
         605             :   end subroutine sat_vapor_press_liq
         606             :
         607             :   !------------------------------------------------------------------------
         608           0 :   subroutine sat_vapor_press_liq_flatau( T_in_K, &
         609             :                                          esat )
         610             :
         611             :   ! Description:
         612             :   !   Computes SVP for water vapor.
         613             :
         614             :   ! References:
         615             :   !   ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
         616             :   !     and Cotton.  (1992)  Journal of Applied Meteorology, Vol. 31,
         617             :   !     pp. 1507--1513
         618             :   !------------------------------------------------------------------------
         619             :
         620             :     use constants_clubb, only: T_freeze_K
         621             :
         622             :     use clubb_precision, only: &
         623             :         core_rknd ! Variable(s)
         624             :
         625             :     implicit none
         626             :
         627             :     ! Constant parameters
         628             :
         629             :     ! Relative error norm expansion (-50 to 50 deg_C) from
         630             :     ! Table 3 of pp. 1510 of Flatau et al. 1992 (Water Vapor)
         631             :     ! (The 100 coefficient converts from mb to Pa)
         632             : !   real, dimension(7), parameter :: a = &
         633             : !   100.* (/ 6.11176750,      0.443986062,     0.143053301E-01, &
         634             : !            0.265027242E-03, 0.302246994E-05, 0.203886313E-07, &
         635             : !            0.638780966E-10 /)
         636             :
         637             :     ! Relative error norm expansion (-85 to 70 deg_C) from
         638             :     ! Table 4 of pp. 1511 of Flatau et al.
         639             :     !real( kind = core_rknd ), dimension(9), parameter :: a = &
         640             :     !100._core_rknd * &
         641             :     !  Commented out because the form has been redone, causing these number to no longer be needed,
         642             :     !  leaving them in for now for reference.
         643             :     !         (/ 6.11583699_core_rknd,      0.444606896_core_rknd,     0.143177157E-01_core_rknd, &
         644             :     !         0.264224321E-03_core_rknd, 0.299291081E-05_core_rknd, 0.203154182E-07_core_rknd, &
         645             :     !         0.702620698E-10_core_rknd, 0.379534310E-13_core_rknd,-0.321582393E-15_core_rknd /)
         646             :
         647             :     real( kind = core_rknd ), parameter :: min_T_in_C = -85._core_rknd ! [deg_C]
         648             :
         649             :     ! ---------------------- Input Variables ----------------------
         650             :     real( kind = core_rknd ), intent(in) :: &
         651             :       T_in_K   ! Temperature   [K]
         652             :
         653             :     ! ---------------------- Output Variables ----------------------
         654             :     real( kind = core_rknd ), intent(out) :: &
         655             :       esat  ! Saturation vapor pressure over water [Pa]
         656             :
         657             :     ! ---------------------- Local Variables ----------------------
         658             :     real( kind = core_rknd ) :: T_in_C, T_in_C_sqd
         659             :
         660             :     ! ---------------------- Begin Code ----------------------
         661             :
         662             :     ! Determine deg K - 273.15
         663           0 :     T_in_C = T_in_K - T_freeze_K
         664             :
         665             :     ! Since this approximation is only good out to -85 degrees Celsius we
         666             :     ! truncate the result here (Flatau, et al. 1992)
         667           0 :     T_in_C = max( T_in_C, min_T_in_C )
         668             :
         669             :     ! Polynomial approx. (Flatau, et al. 1992)
         670             :
         671             :     ! This is the generalized formula but is not computationally efficient.
         672             :     ! Based on Wexler's expressions(2.1)-(2.4) (See Flatau et al. p 1508)
         673             :     ! e_{sat} = a_1 + a_2 ( T - T_0 ) + ... + a_{n+1} ( T - T_0 )^n
         674             :
         675             :     ! esat = a(1)
         676             :
         677             :     ! do i = 2, size( a ) , 1
         678             :     !   esat = esat + a(i) * ( T_in_C )**(i-1)
         679             :     ! end do
         680             :
         681             :     ! The 8th order polynomial fit.  When running deep
         682             :     ! convective cases I noticed that absolute temperature often dips below
         683             :     ! -50 deg_C at higher altitudes, where the 6th order approximation is
         684             :     ! not accurate.  -dschanen 20 Nov 2008
         685             :     !esat = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
         686             :     !*( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*( a(9) ) ) ) ) ) ) ) )
         687             :
         688             :
         689             :     ! Factoring the polynomial above and changing it into this form allows the cpu
         690             :     ! to complete the calculations out of order. This is because modern cpus can complete
         691             :     ! multiple instructions at once if they do not depend on eachother, in the above case
         692             :     ! each instruction relies on the result of the last. In this version however, the terms
         693             :     ! in the parentheses could potentially be calculated in parallel by different execution
         694             :     ! units in the cpu, then only when those terms are being multiplied together do the
         695             :     ! instructions need to be done one at a time. See clubb issue 834 for more info.
         696             :     !   - Gunther Huebler, Aug 2018
         697           0 :     T_in_C_sqd = T_in_C**2
         698             :
         699             :     esat = &
         700             :      - 3.21582393e-14_core_rknd * ( T_in_C - 646.5835252598777_core_rknd ) &
         701             :        * ( T_in_C + 90.72381630364440_core_rknd ) &
         702             :        * ( T_in_C_sqd + 111.0976961559954_core_rknd * T_in_C + 6459.629194243118_core_rknd ) &
         703             :        * ( T_in_C_sqd + 152.3131930092453_core_rknd * T_in_C + 6499.774954705265_core_rknd ) &
         704           0 :        * ( T_in_C_sqd + 174.4279584934021_core_rknd * T_in_C + 7721.679732114084_core_rknd )
         705             :
         706           0 :     return
         707             :   end subroutine sat_vapor_press_liq_flatau
         708             :
         709             :
         710             :   !------------------------------------------------------------------------
         711           0 :   subroutine sat_vapor_press_liq_bolton( T_in_K, &
         712             :                                          esat )
         713             :   ! Description:
         714             :   !   Computes SVP for water vapor.
         715             :   ! References:
         716             :   !   Bolton 1980
         717             :   !------------------------------------------------------------------------
         718             :
         719             :     use constants_clubb, only: T_freeze_K
         720             :
         721             :     use clubb_precision, only: &
         722             :         core_rknd ! Variable(s)
         723             :
         724             :     implicit none
         725             :
         726             :     ! --------------------- Input Variables ---------------------
         727             :     real( kind = core_rknd ), intent(in) :: &
         728             :       T_in_K   ! Temperature   [K]
         729             :
         730             :     ! --------------------- Output Variables ---------------------
         731             :     real( kind = core_rknd ), intent(out) :: &
         732             :       esat  ! Saturation vapor pressure over water [Pa]
         733             :
         734             :
         735             :     ! --------------------------- Begin Code ---------------------------
         736             :
         737             :     ! (Bolton 1980) approx.
         738             :     ! Generally this more computationally expensive than the Flatau polnomial expansion
         739             :     esat = 611.2_core_rknd &
         740             :                 * exp( (17.67_core_rknd *(T_in_K-T_freeze_K))  &
         741           0 :                        / (T_in_K-29.65_core_rknd) ) ! Known magic number
         742             :
         743           0 :     return
         744             :   end subroutine sat_vapor_press_liq_bolton
         745             :
         746             :
         747             :   ! ---> h1g, 2010-06-16
         748             :   !------------------------------------------------------------------------
         749           0 :   subroutine sat_vapor_press_liq_gfdl( T_in_K, &
         750             :                                        esat )
         751             :   ! Description:
         752             :   ! copy from "GFDL polysvp.F90"
         753             :   !  Compute saturation vapor pressure with respect to liquid  by using
         754             :   ! function from Goff and Gratch (1946)
         755             :
         756             :   !  Polysvp returned in units of pa.
         757             :   !  T_in_K  is input in units of K.
         758             :   !------------------------------------------------------------------------
         759             :
         760             :     use clubb_precision, only: &
         761             :         core_rknd ! Variable(s)
         762             :
         763             :     implicit none
         764             :
         765             :     ! --------------------------- Input Variables ---------------------------
         766             :     real( kind = core_rknd ), intent(in) :: &
         767             :       T_in_K   ! Absolute temperature   [K]
         768             :
         769             :     ! --------------------------- Output Variables ---------------------------
         770             :     real( kind = core_rknd ), intent(out) :: &
         771             :       esat  ! Saturation vapor pressure over water [Pa]
         772             :
         773             :     ! --------------------------- Local Variables ---------------------------
         774             :     real( kind = core_rknd ), parameter :: &
         775             :        min_T_in_K = 203.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
         776             :
         777             :     real( kind = core_rknd ) :: &
         778             :        T_in_K_clipped        ! Absolute temperature with minimum threshold applied [K]
         779             :
         780             :     ! --------------------------- Begin Code ---------------------------
         781             :
         782             :     ! Since the Goff-Gratch approximation is valid only down to -70 degrees Celsius,
         783             :     !   we threshold the temperature.  This will yield a minimal saturation at
         784             :     !   cold temperatures.
         785           0 :     T_in_K_clipped = max( min_T_in_K, T_in_K )
         786             :
         787             :     ! Goff Gratch equation, uncertain below -70 C
         788             :
         789             :     esat = 10._core_rknd**(-7.90298_core_rknd*(373.16_core_rknd/T_in_K_clipped-1._core_rknd)+ &
         790             :          5.02808_core_rknd*log10(373.16_core_rknd/T_in_K_clipped)- &
         791             :          1.3816e-7_core_rknd*(10._core_rknd**(11.344_core_rknd &
         792             :            *(1._core_rknd-T_in_K_clipped/373.16_core_rknd))-1._core_rknd)+ &
         793             :          8.1328e-3_core_rknd*(10._core_rknd**(-3.49149_core_rknd &
         794             :            *(373.16_core_rknd/T_in_K_clipped-1._core_rknd))-1._core_rknd)+ &
         795           0 :          log10(1013.246_core_rknd))*100._core_rknd ! Known magic number
         796             :
         797           0 :     return
         798             :   end subroutine sat_vapor_press_liq_gfdl
         799             : ! <--- h1g, 2010-06-16
         800             :
         801             :   !------------------------------------------------------------------------
         802           0 :   elemental function sat_vapor_press_liq_lookup( T_in_K ) result ( esat )
         803             :
         804             : ! Description:
         805             : !   Computes SVP for water vapor, using a lookup table.
         806             : !
         807             : !   The lookup table was constructed using the Flatau approximation.
         808             :
         809             : ! References:
         810             : !   ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
         811             : !     and Cotton.  (1992)  Journal of Applied Meteorology, Vol. 31,
         812             : !     pp. 1507--1513
         813             : !------------------------------------------------------------------------
         814             :
         815             :     implicit none
         816             :
         817             :     ! External
         818             :     intrinsic :: max, min, int, anint
         819             :
         820             :     ! Input Variables
         821             :     real( kind = core_rknd ), intent(in) :: T_in_K   ! Temperature   [K]
         822             :
         823             :     ! Output Variables
         824             :     real( kind = core_rknd ) :: esat  ! Saturation vapor pressure over water [Pa]
         825             :
         826             :     ! Local Variables
         827             :     integer :: T_in_K_int
         828             :
         829             :     ! ---- Begin Code ----
         830             :
         831           0 :     T_in_K_int = int( anint( T_in_K ) )
         832             :
         833             :     ! Since this approximation is only good out to -85 degrees Celsius we
         834             :     ! truncate the result here
         835           0 :     T_in_K_int = min( max( T_in_K_int, 188 ), 343 )
         836             :
         837             :     ! Use the lookup table to determine the saturation vapor pressure.
         838           0 :     esat = svp_liq_lookup_table( T_in_K_int )
         839             :
         840             :     return
         841             :   end function sat_vapor_press_liq_lookup
         842             :
         843             :   !------------------------------------------------------------------------
         844             :   ! Wrapped in interface sat_mixrat_ice
         845           0 :   function sat_mixrat_ice_k( p_in_Pa, T_in_K, saturation_formula )
         846             :
         847             :   ! Description:
         848             :   !   Used to compute the saturation mixing ratio of ice.
         849             :
         850             :   ! References:
         851             :   !   Formula from Emanuel 1994, 4.4.15
         852             :   !-------------------------------------------------------------------------
         853             :
         854             :     use constants_clubb, only: &
         855             :         ep ! Variable(s)
         856             :
         857             :     use clubb_precision, only: &
         858             :         core_rknd ! Variable(s)
         859             :
         860             :     implicit none
         861             :
         862             :     ! ------------------------ Input Variables ------------------------
         863             :     real( kind = core_rknd ), intent(in) :: &
         864             :       p_in_Pa, &          ! Pressure [Pa]
         865             :       T_in_K              ! Temperature [K]
         866             :
         867             :     integer, intent(in) :: &
         868             :       saturation_formula
         869             :
         870             :     ! ------------------------ Output Variables ------------------------
         871             :     real( kind = core_rknd ) :: &
         872             :       sat_mixrat_ice_k
         873             :
         874             :     ! ------------------------ Local Variables ------------------------
         875             :     real( kind = core_rknd ), dimension(1,1) ::  &
         876             :       p_in_Pa_col,  &
         877             :       T_in_K_col
         878             :
         879             :     real( kind = core_rknd ), dimension(1,1) ::  &
         880             :       sat_mixrat_ice_col
         881             :
         882             :     integer :: i, k  ! Loop indices
         883             :
         884             :     ! ------------------------ Begin Code ------------------------
         885             :
         886             :     ! Copy inputs to 2D arrays
         887           0 :     p_in_Pa_col(1,1) = p_in_Pa
         888           0 :     T_in_K_col(1,1) = T_in_K
         889             :
         890             :     ! Call 2D version
         891           0 :     sat_mixrat_ice_col = sat_mixrat_ice_2D( 1, 1, p_in_Pa_col, T_in_K_col, saturation_formula )
         892             :
         893             :     ! Copy 2D result into output
         894           0 :     sat_mixrat_ice_k = sat_mixrat_ice_col( 1, 1 )
         895             :
         896             :     return
         897             :   end function sat_mixrat_ice_k
         898             :
         899             :   !------------------------------------------------------------------------
         900             :   ! Wrapped in interface sat_mixrat_ice
         901    35740224 :   function sat_mixrat_ice_2D( nz, ngrdcol, p_in_Pa, T_in_K, saturation_formula )
         902             :
         903             :   ! Description:
         904             :   !   Used to compute the saturation mixing ratio of ice.
         905             :
         906             :   ! References:
         907             :   !   Formula from Emanuel 1994, 4.4.15
         908             :   !-------------------------------------------------------------------------
         909             :
         910             :     use constants_clubb, only: &
         911             :         T_freeze_K, & ! Variable(s)
         912             :         ep
         913             :
         914             :     use clubb_precision, only: &
         915             :         core_rknd ! Variable(s)
         916             :
         917             :     implicit none
         918             :
         919             :     ! ------------------------ Input Variables ------------------------
         920             :     integer, intent(in) :: &
         921             :       nz, &
         922             :       ngrdcol
         923             :
         924             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         925             :       p_in_Pa, &          ! Pressure [Pa]
         926             :       T_in_K              ! Temperature [K]
         927             :
         928             :     integer, intent(in) :: &
         929             :       saturation_formula
         930             :
         931             :     ! ------------------------ Output Variables ------------------------
         932             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         933             :       sat_mixrat_ice_2D
         934             :
         935             :     ! ------------------------ Local Variables ------------------------
         936             :     ! Relative error norm expansion (-90 to 0 deg_C) from
         937             :     ! Table 4 of pp. 1511 of Flatau et al. 1992 (Ice)
         938             :     real( kind = core_rknd ), dimension(9), parameter :: a = &
         939             :       100._core_rknd * (/ 6.09868993_core_rknd, 0.499320233_core_rknd, 0.184672631E-01_core_rknd, &
         940             :                 0.402737184E-03_core_rknd, 0.565392987E-05_core_rknd, 0.521693933E-07_core_rknd, &
         941             :                 0.307839583E-09_core_rknd, 0.105785160E-11_core_rknd, 0.161444444E-14_core_rknd /)
         942             :
         943             :     real( kind = core_rknd ), parameter :: &
         944             :       min_T_in_C = -90._core_rknd,  & ! [deg_C]
         945             :       min_T_in_K = 173.15_core_rknd   ! Lowest temperature at which Goff-Gratch is valid [K]
         946             :
         947             :     real( kind = core_rknd ) :: &
         948             :       T_in_C,         & ! Temperature [deg_C]
         949             :       T_in_K_clipped    ! Absolute temperature with minimum threshold applied [K]
         950             :
         951             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
         952    35740224 :       esat_ice
         953             :
         954             :     integer :: i, k  ! Loop indices
         955             :
         956             :     ! ------------------------ Begin Code ------------------------
         957             :
         958             :     !$acc data create( esat_ice ) &
         959             :     !$acc      copyin( p_in_Pa, T_in_K ) &
         960             :     !$acc      copyout( sat_mixrat_ice_2D )
         961             :
         962             :     ! Determine the SVP for the given temperature
         963    35740224 :     select case ( saturation_formula )
         964             :     case ( saturation_bolton )
         965             :
         966             :       ! Using the Bolton 1980 approximations for SVP over ice
         967             :       !$acc parallel loop gang vector collapse(2) default(present)
         968           0 :       do i = 1, ngrdcol
         969           0 :         do k = 1, nz
         970             :
         971             :           ! Exponential approx.
         972           0 :           esat_ice(i,k) = 100.0_core_rknd * exp( 23.33086_core_rknd - &
         973           0 :             (6111.72784_core_rknd/T_in_K(i,k)) + (0.15215_core_rknd*log( T_in_K(i,k) )) )
         974             :
         975             :         end do
         976             :       end do
         977             :       !$acc end parallel loop
         978             :
         979             :     case ( saturation_flatau )
         980             :
         981             :       ! Using the Flatau, et al. polynomial approximation for SVP over ice
         982             :       !$acc parallel loop gang vector collapse(2) default(present)
         983           0 :       do i = 1, ngrdcol
         984           0 :         do k = 1, nz
         985             :
         986             :           ! Determine deg K - 273.15
         987           0 :           T_in_C = T_in_K(i,k) - T_freeze_K
         988             :
         989             :           ! Since this approximation is only good out to -90 degrees Celsius we
         990             :           ! truncate the result here (Flatau, et al. 1992)
         991           0 :           T_in_C = max( T_in_C, min_T_in_C )
         992             :
         993             :           ! Polynomial approx. (Flatau, et al. 1992)
         994             :           !   esati = a(1)
         995             :
         996             :           !   do i = 2, size( a ), 1
         997             :           !     esati = esati + a(i) * ( T_in_C )**(i-1)
         998             :           !   end do
         999             :
        1000             :           esat_ice(i,k) = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
        1001           0 :           *( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*a(9) ) ) ) ) ) ) )
        1002             :
        1003             :         end do
        1004             :       end do
        1005             :       !$acc end parallel loop
        1006             :
        1007             : ! ---> h1g, 2010-06-16
        1008             :     case ( saturation_gfdl )
        1009             :
        1010             :       ! Using GFDL polynomial approximation for SVP with respect to ice
        1011             :       !$acc parallel loop gang vector collapse(2) default(present)
        1012   596778624 :       do i = 1, ngrdcol
        1013 48285042624 :         do k = 1, nz
        1014             :
        1015             :           ! Since the Goff-Gratch ice approximation is valid only down to -100 degrees Celsius,
        1016             :           !   we threshold the temperature.  This will yield a minimal saturation at
        1017             :           !   cold temperatures.
        1018 47688264000 :           T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
        1019             :
        1020             :           ! Goff Gratch equation (good down to -100 C)
        1021             :
        1022             :           esat_ice(i,k) = 10._core_rknd**(-9.09718_core_rknd* &
        1023             :                       (273.16_core_rknd/T_in_K_clipped-1._core_rknd)-3.56654_core_rknd* &
        1024             :                     log10(273.16_core_rknd/T_in_K_clipped)+0.876793_core_rknd* &
        1025             :                       (1._core_rknd-T_in_K_clipped/273.16_core_rknd)+ &
        1026 48249302400 :                     log10(6.1071_core_rknd))*100._core_rknd ! Known magic number
        1027             :         end do
        1028             :       end do
        1029             :       !$acc end parallel loop
        1030             :
        1031             : ! <--- h1g, 2010-06-16
        1032             :
        1033             :     case default
        1034             :
        1035             :       ! Undefined approximation
        1036    35740224 :       esat_ice = -99999.999_core_rknd
        1037             :
        1038             :     end select
        1039             :
        1040             :
        1041             :     !$acc parallel loop gang vector collapse(2) default(present)
        1042  3073659264 :     do k = 1, nz
        1043 50761923264 :       do i = 1, ngrdcol
        1044             :
        1045             :         ! If esat_ice exceeds the air pressure, then assume esat_ice~=0.5*pressure
        1046             :         !   and set rsat = ep = 0.622
        1047 50726183040 :         if ( p_in_Pa(i,k)-esat_ice(i,k) < 1.0_core_rknd ) then
        1048   890022639 :           sat_mixrat_ice_2D(i,k) = ep
        1049             :         else
        1050             :
        1051             :           if ( I_sat_sphum )  then   ! h1g, 2010-06-18 begin mod
        1052             :
        1053             :             ! GFDL uses specific humidity
        1054             :             ! Formula for Saturation Specific Humidity
        1055             :             sat_mixrat_ice_2D(i,k) = ep * ( esat_ice(i,k) &
        1056             :                                   / ( p_in_Pa(i,k) - (1.0_core_rknd-ep) * esat_ice(i,k) ) )
        1057             :           else
        1058             :
        1059             :             ! Formula for Saturation Mixing Ratio:
        1060             :             !
        1061             :             ! rs = (epsilon) * [ esat / ( p - esat ) ];
        1062             :             ! where epsilon = R_d / R_v
        1063 46798241361 :             sat_mixrat_ice_2D(i,k) = ep * esat_ice(i,k) / ( p_in_Pa(i,k) - esat_ice(i,k) )
        1064             :
        1065             :           endif                     ! h1g, 2010-06-18 end mod
        1066             :         end if
        1067             :
        1068             :       end do
        1069             :     end do
        1070             :     !$acc end parallel loop
        1071             :
        1072             :     !$acc end data
        1073             :
        1074    35740224 :     return
        1075             :   end function sat_mixrat_ice_2D
        1076             :
        1077             :   !------------------------------------------------------------------------
        1078             :   subroutine sat_vapor_press_ice( nz, ngrdcol, T_in_K, saturation_formula, &
        1079             :                                   esat_ice )
        1080             :   !
        1081             :   ! Description:
        1082             :   !   Computes SVP for ice, using one of the various approximations.
        1083             :   !
        1084             :   ! References:
        1085             :   !   None
        1086             :   !------------------------------------------------------------------------
        1087             :
        1088             :
        1089             :     use clubb_precision, only: &
        1090             :         core_rknd ! Variable(s)
        1091             :
        1092             :     implicit none
        1093             :
        1094             :     ! ---------------------- Input Variable ----------------------
        1095             :     integer, intent(in) :: &
        1096             :       nz, &
        1097             :       ngrdcol
        1098             :
        1099             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1100             :       T_in_K      ! Temperature     [K]
        1101             :
        1102             :     integer, intent(in) :: &
        1103             :       saturation_formula
        1104             :
        1105             :     ! ---------------------- Output Variable ----------------------
        1106             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1107             :       esat_ice    ! Saturation Vapor Pressure over Ice [Pa]
        1108             :
        1109             :     ! ---------------------- Begin Code ----------------------
        1110             :
        1111             :     select case ( saturation_formula )
        1112             :     case ( saturation_bolton )
        1113             :
        1114             :       ! Using the Bolton 1980 approximations for SVP over ice
        1115             :       call sat_vapor_press_ice_bolton( nz, ngrdcol, T_in_K, &
        1116             :                                        esat_ice )
        1117             :
        1118             :     case ( saturation_flatau )
        1119             :
        1120             :       ! Using the Flatau, et al. polynomial approximation for SVP over ice
        1121             :       call sat_vapor_press_ice_flatau( nz, ngrdcol, T_in_K, &
        1122             :                                        esat_ice )
        1123             :
        1124             : ! ---> h1g, 2010-06-16
        1125             :     case ( saturation_gfdl )
        1126             :
        1127             :       ! Using GFDL polynomial approximation for SVP with respect to ice
        1128             :       call sat_vapor_press_ice_gfdl( nz, ngrdcol, T_in_K, &
        1129             :                                      esat_ice )
        1130             : ! <--- h1g, 2010-06-16
        1131             :
        1132             :     case default
        1133             :
        1134             :       ! Undefined approximation
        1135             :       esat_ice = -99999.999_core_rknd
        1136             :
        1137             :     end select
        1138             :
        1139             :     return
        1140             :
        1141             :   end subroutine sat_vapor_press_ice
        1142             :
        1143             :   !------------------------------------------------------------------------
        1144             :   subroutine sat_vapor_press_ice_flatau( nz, ngrdcol, T_in_K, &
        1145             :                                          esati )
        1146             :   !
        1147             :   ! Description:
        1148             :   !   Computes SVP for ice.
        1149             :   !
        1150             :   ! References:
        1151             :   !   ``Polynomial Fits to Saturation Vapor Pressure'' Falatau, Walko,
        1152             :   !     and Cotton.  (1992)  Journal of Applied Meteorology, Vol. 31,
        1153             :   !     pp. 1507--1513
        1154             :   !------------------------------------------------------------------------
        1155             :     use constants_clubb, only: T_freeze_K
        1156             :
        1157             :     use clubb_precision, only: &
        1158             :         core_rknd ! Variable(s)
        1159             :
        1160             :     implicit none
        1161             :
        1162             :     ! ------------------------ Input Variables ------------------------
        1163             :     integer, intent(in) :: &
        1164             :       nz, &
        1165             :       ngrdcol
        1166             :
        1167             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1168             :       T_in_K   ! Temperature   [deg_K]
        1169             :
        1170             :     ! ------------------------ Output Variables ------------------------
        1171             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1172             :       esati  ! Saturation vapor pressure over ice [Pa]
        1173             :
        1174             :     ! ------------------------ Local Variables ------------------------
        1175             :     ! Relative error norm expansion (-90 to 0 deg_C) from
        1176             :     ! Table 4 of pp. 1511 of Flatau et al. 1992 (Ice)
        1177             :     real( kind = core_rknd ), dimension(9), parameter :: a = &
        1178             :     100._core_rknd * (/ 6.09868993_core_rknd, 0.499320233_core_rknd, 0.184672631E-01_core_rknd, &
        1179             :               0.402737184E-03_core_rknd, 0.565392987E-05_core_rknd, 0.521693933E-07_core_rknd, &
        1180             :               0.307839583E-09_core_rknd, 0.105785160E-11_core_rknd, 0.161444444E-14_core_rknd /)
        1181             :
        1182             :     real( kind = core_rknd ), parameter :: min_T_in_C = -90._core_rknd ! [deg_C]
        1183             :
        1184             :     real( kind = core_rknd ) :: &
        1185             :       T_in_C ! Temperature [deg_C]
        1186             :
        1187             :     integer :: i, k
        1188             :
        1189             :     ! ------------------------ Begin Code ------------------------
        1190             :
        1191             :     do i = 1, ngrdcol
        1192             :       do k = 1, nz
        1193             :
        1194             :         ! Determine deg K - 273.15
        1195             :         T_in_C = T_in_K(i,k) - T_freeze_K
        1196             :
        1197             :         ! Since this approximation is only good out to -90 degrees Celsius we
        1198             :         ! truncate the result here (Flatau, et al. 1992)
        1199             :         T_in_C = max( T_in_C, min_T_in_C )
        1200             :
        1201             :         ! Polynomial approx. (Flatau, et al. 1992)
        1202             :         !   esati = a(1)
        1203             :
        1204             :         !   do i = 2, size( a ), 1
        1205             :         !     esati = esati + a(i) * ( T_in_C )**(i-1)
        1206             :         !   end do
        1207             :
        1208             :         esati(i,k) = a(1) + T_in_C*( a(2) + T_in_C*( a(3) + T_in_C*( a(4) + T_in_C &
        1209             :         *( a(5) + T_in_C*( a(6) + T_in_C*( a(7) + T_in_C*( a(8) + T_in_C*a(9) ) ) ) ) ) ) )
        1210             :
        1211             :       end do
        1212             :     end do
        1213             :
        1214             :     return
        1215             :
        1216             :   end subroutine sat_vapor_press_ice_flatau
        1217             :
        1218             :   !------------------------------------------------------------------------
        1219             :   subroutine sat_vapor_press_ice_bolton( nz, ngrdcol, T_in_K, &
        1220             :                                          esati )
        1221             :   !
        1222             :   ! Description:
        1223             :   !   Computes SVP for ice.
        1224             :   !
        1225             :   ! References:
        1226             :   !   Bolton 1980
        1227             :   !------------------------------------------------------------------------
        1228             :
        1229             :     use clubb_precision, only: &
        1230             :         core_rknd ! Variable(s)
        1231             :
        1232             :     implicit none
        1233             :
        1234             :     ! ------------------------ Input Variables ------------------------
        1235             :     integer, intent(in) :: &
        1236             :       nz, &
        1237             :       ngrdcol
        1238             :
        1239             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1240             :       T_in_K   ! Temperature   [deg_K]
        1241             :
        1242             :     ! ------------------------ Output Variables ------------------------
        1243             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1244             :       esati  ! Saturation vapor pressure over ice [Pa]
        1245             :
        1246             :     ! ------------------------ Local Variables ------------------------
        1247             :     integer :: i, k
        1248             :
        1249             :     ! ------------------------ Begin Code ------------------------
        1250             :
        1251             :     do i = 1, ngrdcol
        1252             :       do k = 1, nz
        1253             :
        1254             :         ! Exponential approx.
        1255             :         esati(i,k) = 100.0_core_rknd * exp( 23.33086_core_rknd - &
        1256             :           (6111.72784_core_rknd/T_in_K(i,k)) + (0.15215_core_rknd*log( T_in_K(i,k) )) )
        1257             :
        1258             :       end do
        1259             :     end do
        1260             :
        1261             :     return
        1262             :
        1263             :   end subroutine sat_vapor_press_ice_bolton
        1264             :
        1265             :
        1266             :   ! ---> h1g, 2010-06-16
        1267             :   !------------------------------------------------------------------------
        1268             :   subroutine sat_vapor_press_ice_gfdl( nz, ngrdcol, T_in_K, &
        1269             :                                        esati )
        1270             :   ! Description:
        1271             :   ! copy from "GFDL polysvp.F90"
        1272             :   !  Compute saturation vapor pressure with respect to liquid by using
        1273             :   ! function from Goff and Gratch (1946)
        1274             :   !
        1275             :   !  Polysvp returned in units of pa.
        1276             :   !  T_in_K is input in units of K.
        1277             :   !------------------------------------------------------------------------
        1278             :
        1279             :     use clubb_precision, only: &
        1280             :         core_rknd ! Variable(s)
        1281             :
        1282             :     implicit none
        1283             :
        1284             :     ! ------------------------ Input Variables ------------------------
        1285             :     integer, intent(in) :: &
        1286             :       nz, &
        1287             :       ngrdcol
        1288             :
        1289             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
        1290             :       T_in_K   ! Temperature   [deg_K]
        1291             :
        1292             :     ! ------------------------ Output Variables ------------------------
        1293             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
        1294             :       esati  ! Saturation vapor pressure over ice [Pa]
        1295             :
        1296             :     ! ------------------------ Local Variables ------------------------
        1297             :     real( kind = core_rknd ), parameter :: &
        1298             :        min_T_in_K = 173.15_core_rknd ! Lowest temperature at which Goff-Gratch is valid [K]
        1299             :
        1300             :     real( kind = core_rknd ) :: &
        1301             :        T_in_K_clipped        ! Absolute temperature with minimum threshold applied [K]
        1302             :
        1303             :     integer :: i, k
        1304             :
        1305             :     ! ------------------------ Begin Code ------------------------
        1306             :
        1307             :     do i = 1, ngrdcol
        1308             :       do k = 1, nz
        1309             :
        1310             :         ! Since the Goff-Gratch ice approximation is valid only down to -100 degrees Celsius,
        1311             :         !   we threshold the temperature.  This will yield a minimal saturation at
        1312             :         !   cold temperatures.
        1313             :         T_in_K_clipped = max( min_T_in_K, T_in_K(i,k) )
        1314             :
        1315             :         ! Goff Gratch equation (good down to -100 C)
        1316             :
        1317             :         esati(i,k) = 10._core_rknd**(-9.09718_core_rknd* &
        1318             :                     (273.16_core_rknd/T_in_K_clipped-1._core_rknd)-3.56654_core_rknd* &
        1319             :                   log10(273.16_core_rknd/T_in_K_clipped)+0.876793_core_rknd* &
        1320             :                     (1._core_rknd-T_in_K_clipped/273.16_core_rknd)+ &
        1321             :                   log10(6.1071_core_rknd))*100._core_rknd ! Known magic number
        1322             :       end do
        1323             :     end do
        1324             :
        1325             :     return
        1326             :
        1327             :   end subroutine sat_vapor_press_ice_gfdl
        1328             : ! <--- h1g, 2010-06-16
        1329             :
        1330             : !-------------------------------------------------------------------------
        1331           0 :   function rcm_sat_adj( thlm, rtm, p_in_Pa, exner, saturation_formula ) result ( rcm )
        1332             :
        1333             :     ! Description:
        1334             :     !
        1335             :     !   This function uses an iterative method to find the value of rcm
        1336             :     !   from an initial profile that has saturation at some point.
        1337             :     !
        1338             :     ! References:
        1339             :     !   None
        1340             :     !-------------------------------------------------------------------------
        1341             :
        1342             :     use clubb_precision, only: &
        1343             :         core_rknd ! Variable(s)
        1344             :
        1345             :     use constants_clubb, only: &
        1346             :         Cp,            & ! Variable(s)
        1347             :         Lv,            &
        1348             :         zero_threshold
        1349             :
        1350             :     implicit none
        1351             :
        1352             :     ! Local Constant(s)
        1353             :     real( kind = core_rknd ), parameter :: &
        1354             :       tolerance = 0.001_core_rknd ! Tolerance on theta calculation [K]
        1355             :
        1356             :     integer, parameter :: &
        1357             :       itermax = 1000000 ! Maximum interations
        1358             :
        1359             :     ! External
        1360             :     intrinsic :: max, abs
        1361             :
        1362             :     ! Input Variable(s)
        1363             :     real( kind = core_rknd ), intent(in) :: &
        1364             :       thlm,    & ! Liquid Water Potential Temperature [K]
        1365             :       rtm,     & ! Total Water Mixing Ratio       [kg/kg]
        1366             :       p_in_Pa, & ! Pressure                          [Pa]
        1367             :       exner      ! Exner function                     [-]
        1368             :
        1369             :     integer, intent(in) :: &
        1370             :       saturation_formula
        1371             :
        1372             :     ! Output Variable(s)
        1373             :     real( kind = core_rknd ) :: rcm ! Cloud water mixing ratio      [kg/kg]
        1374             :
        1375             :     ! Local Variable(s)
        1376             :     real( kind = core_rknd ) :: &
        1377             :       theta, answer, too_low, too_high, & ! [K]
        1378             :       rsat
        1379             :
        1380             :     integer :: iteration
        1381             :
        1382             :     ! ----- Begin Code -----
        1383             :
        1384             :     ! Default initialization
        1385           0 :     theta = thlm
        1386           0 :     too_high = 0.0_core_rknd
        1387           0 :     too_low  = 0.0_core_rknd
        1388             :
        1389           0 :     do iteration = 1, itermax, 1
        1390             :
        1391             :       answer = theta - (Lv/(Cp*exner)) &
        1392           0 :                        *(MAX( rtm - sat_mixrat_liq(p_in_Pa,theta*exner, saturation_formula), zero_threshold ))
        1393             :
        1394           0 :       if ( ABS(answer - thlm) <= tolerance ) then
        1395             :         exit
        1396           0 :       else if ( answer - thlm > tolerance ) then
        1397             :         too_high = theta
        1398           0 :       else if ( thlm - answer > tolerance ) THEN
        1399           0 :         too_low = theta
        1400             :       end if
        1401             :
        1402             :       ! For the first timestep, be sure to set a "too_high"
        1403             :       ! that is "way too high."
        1404           0 :       if ( iteration == 1 ) then
        1405           0 :         too_high = theta + 20.0_core_rknd
        1406             :       end if
        1407             :
        1408           0 :       theta = (too_low + too_high)/2.0_core_rknd
        1409             :
        1410             :     end do ! 1..itermax
        1411             :
        1412           0 :     if ( iteration == itermax ) then
        1413             :       ! Magic Eric Raut added to remove compiler warning (clearly this value is not used)
        1414           0 :       rcm = 0.0_core_rknd
        1415             :
        1416           0 :       error stop "Error in rcm_sat_adj: could not determine rcm"
        1417             :     else
        1418           0 :       rcm = MAX( rtm - sat_mixrat_liq( p_in_Pa, theta*exner, saturation_formula), zero_threshold )
        1419             :       return
        1420             :     end if
        1421             :
        1422             :   end function rcm_sat_adj
        1423             :
        1424             : end module saturation

Generated by: LCOV version 1.14