1 2! KGEN-generated Fortran source file 3! 4! Filename : wv_sat_methods.F90 5! Generated at: 2015-03-31 09:44:41 6! KGEN version: 0.4.5 7 8 9 10 MODULE wv_sat_methods 11 ! This portable module contains all 1 methods for estimating 12 ! the saturation vapor pressure of water. 13 ! 14 ! wv_saturation provides 1-specific interfaces and utilities 15 ! based on these formulae. 16 ! 17 ! Typical usage of this module: 18 ! 19 ! Init: 20 ! call wv_sat_methods_init(r8, <constants>, errstring) 21 ! 22 ! Get scheme index from a name string: 23 ! scheme_idx = wv_sat_get_scheme_idx(scheme_name) 24 ! if (.not. wv_sat_valid_idx(scheme_idx)) <throw some error> 25 ! 26 ! Get pressures: 27 ! es = wv_sat_svp_water(t, scheme_idx) 28 ! es = wv_sat_svp_ice(t, scheme_idx) 29 ! 30 ! Use ice/water transition range: 31 ! es = wv_sat_svp_trice(t, ttrice, scheme_idx) 32 ! 33 ! Note that elemental functions cannot be pointed to, nor passed 34 ! as arguments. If you need to do either, it is recommended to 35 ! wrap the function so that it can be given an explicit (non- 36 ! elemental) interface. 37 IMPLICIT NONE 38 PRIVATE 39 INTEGER, parameter :: r8 = selected_real_kind(12) ! 8 byte real 40 REAL(KIND=r8) :: tmelt ! Melting point of water at 1 atm (K) 41 REAL(KIND=r8) :: h2otrip ! Triple point temperature of water (K) 42 REAL(KIND=r8) :: tboil ! Boiling point of water at 1 atm (K) 43 ! Ice-water transition range 44 REAL(KIND=r8) :: epsilo ! Ice-water transition range 45 REAL(KIND=r8) :: omeps ! 1._r8 - epsilo 46 ! Indices representing individual schemes 47 INTEGER, parameter :: oldgoffgratch_idx = 0 48 INTEGER, parameter :: goffgratch_idx = 1 49 INTEGER, parameter :: murphykoop_idx = 2 50 INTEGER, parameter :: bolton_idx = 3 51 ! Index representing the current default scheme. 52 INTEGER, parameter :: initial_default_idx = goffgratch_idx 53 INTEGER :: default_idx = initial_default_idx 54 PUBLIC wv_sat_svp_water 55 PUBLIC wv_sat_svp_ice 56 ! pressure -> humidity conversion 57 PUBLIC wv_sat_svp_to_qsat 58 ! Combined qsat operations 59 PUBLIC wv_sat_qsat_water 60 PUBLIC wv_sat_qsat_ice 61 PUBLIC kgen_read_externs_wv_sat_methods 62 CONTAINS 63 64 ! write subroutines 65 ! No subroutines 66 67 ! module extern variables 68 69 SUBROUTINE kgen_read_externs_wv_sat_methods(kgen_unit) 70 INTEGER, INTENT(IN) :: kgen_unit 71 READ(UNIT=kgen_unit) tmelt 72 READ(UNIT=kgen_unit) h2otrip 73 READ(UNIT=kgen_unit) tboil 74 READ(UNIT=kgen_unit) epsilo 75 READ(UNIT=kgen_unit) omeps 76 READ(UNIT=kgen_unit) default_idx 77 END SUBROUTINE kgen_read_externs_wv_sat_methods 78 79 !--------------------------------------------------------------------- 80 ! ADMINISTRATIVE FUNCTIONS 81 !--------------------------------------------------------------------- 82 ! Get physical constants 83 84 ! Look up index by name. 85 86 ! Check validity of an index from the above routine. 87 88 ! Set default scheme (otherwise, Goff & Gratch is default) 89 ! Returns a logical representing success (.true.) or 90 ! failure (.false.). 91 92 ! Reset default scheme to initial value. 93 ! The same thing can be accomplished with wv_sat_set_default; 94 ! the real reason to provide this routine is to reset the 95 ! module for testing purposes. 96 97 !--------------------------------------------------------------------- 98 ! UTILITIES 99 !--------------------------------------------------------------------- 100 ! Get saturation specific humidity given pressure and SVP. 101 ! Specific humidity is limited to range 0-1. 102 103 elemental FUNCTION wv_sat_svp_to_qsat(es, p) RESULT ( qs ) 104 REAL(KIND=r8), intent(in) :: es ! SVP 105 REAL(KIND=r8), intent(in) :: p ! Current pressure. 106 REAL(KIND=r8) :: qs 107 ! If pressure is less than SVP, set qs to maximum of 1. 108 IF ((p - es) <= 0._r8) THEN 109 qs = 1.0_r8 110 ELSE 111 qs = epsilo*es / (p - omeps*es) 112 END IF 113 END FUNCTION wv_sat_svp_to_qsat 114 115 elemental SUBROUTINE wv_sat_qsat_water(t, p, es, qs, idx) 116 !------------------------------------------------------------------! 117 ! Purpose: ! 118 ! Calculate SVP over water at a given temperature, and then ! 119 ! calculate and return saturation specific humidity. ! 120 !------------------------------------------------------------------! 121 ! Inputs 122 REAL(KIND=r8), intent(in) :: t ! Temperature 123 REAL(KIND=r8), intent(in) :: p ! Pressure 124 ! Outputs 125 REAL(KIND=r8), intent(out) :: es ! Saturation vapor pressure 126 REAL(KIND=r8), intent(out) :: qs ! Saturation specific humidity 127 INTEGER, intent(in), optional :: idx ! Scheme index 128 es = wv_sat_svp_water(t, idx) 129 qs = wv_sat_svp_to_qsat(es, p) 130 ! Ensures returned es is consistent with limiters on qs. 131 es = min(es, p) 132 END SUBROUTINE wv_sat_qsat_water 133 134 elemental SUBROUTINE wv_sat_qsat_ice(t, p, es, qs, idx) 135 !------------------------------------------------------------------! 136 ! Purpose: ! 137 ! Calculate SVP over ice at a given temperature, and then ! 138 ! calculate and return saturation specific humidity. ! 139 !------------------------------------------------------------------! 140 ! Inputs 141 REAL(KIND=r8), intent(in) :: t ! Temperature 142 REAL(KIND=r8), intent(in) :: p ! Pressure 143 ! Outputs 144 REAL(KIND=r8), intent(out) :: es ! Saturation vapor pressure 145 REAL(KIND=r8), intent(out) :: qs ! Saturation specific humidity 146 INTEGER, intent(in), optional :: idx ! Scheme index 147 es = wv_sat_svp_ice(t, idx) 148 qs = wv_sat_svp_to_qsat(es, p) 149 ! Ensures returned es is consistent with limiters on qs. 150 es = min(es, p) 151 END SUBROUTINE wv_sat_qsat_ice 152 153 !--------------------------------------------------------------------- 154 ! SVP INTERFACE FUNCTIONS 155 !--------------------------------------------------------------------- 156 157 elemental FUNCTION wv_sat_svp_water(t, idx) RESULT ( es ) 158 REAL(KIND=r8), intent(in) :: t 159 INTEGER, intent(in), optional :: idx 160 REAL(KIND=r8) :: es 161 INTEGER :: use_idx 162 IF (present(idx)) THEN 163 use_idx = idx 164 ELSE 165 use_idx = default_idx 166 END IF 167 SELECT CASE ( use_idx ) 168 CASE ( goffgratch_idx ) 169 es = goffgratch_svp_water(t) 170 CASE ( murphykoop_idx ) 171 es = murphykoop_svp_water(t) 172 CASE ( oldgoffgratch_idx ) 173 es = oldgoffgratch_svp_water(t) 174 CASE ( bolton_idx ) 175 es = bolton_svp_water(t) 176 END SELECT 177 END FUNCTION wv_sat_svp_water 178 179 elemental FUNCTION wv_sat_svp_ice(t, idx) RESULT ( es ) 180 REAL(KIND=r8), intent(in) :: t 181 INTEGER, intent(in), optional :: idx 182 REAL(KIND=r8) :: es 183 INTEGER :: use_idx 184 IF (present(idx)) THEN 185 use_idx = idx 186 ELSE 187 use_idx = default_idx 188 END IF 189 SELECT CASE ( use_idx ) 190 CASE ( goffgratch_idx ) 191 es = goffgratch_svp_ice(t) 192 CASE ( murphykoop_idx ) 193 es = murphykoop_svp_ice(t) 194 CASE ( oldgoffgratch_idx ) 195 es = oldgoffgratch_svp_ice(t) 196 CASE ( bolton_idx ) 197 es = bolton_svp_water(t) 198 END SELECT 199 END FUNCTION wv_sat_svp_ice 200 201 !--------------------------------------------------------------------- 202 ! SVP METHODS 203 !--------------------------------------------------------------------- 204 ! Goff & Gratch (1946) 205 206 elemental FUNCTION goffgratch_svp_water(t) RESULT ( es ) 207 REAL(KIND=r8), intent(in) :: t ! Temperature in Kelvin 208 REAL(KIND=r8) :: es ! SVP in Pa 209 ! uncertain below -70 C 210 es = 10._r8**(-7.90298_r8*(tboil/t-1._r8)+ 5.02808_r8*log10(tboil/t)- 1.3816e-7_r8*(10._r8**(11.344_r8*(& 211 1._r8-t/tboil))-1._r8)+ 8.1328e-3_r8*(10._r8**(-3.49149_r8*(tboil/t-1._r8))-1._r8)+ log10(1013.246_r8))*100._r8 212 END FUNCTION goffgratch_svp_water 213 214 elemental FUNCTION goffgratch_svp_ice(t) RESULT ( es ) 215 REAL(KIND=r8), intent(in) :: t ! Temperature in Kelvin 216 REAL(KIND=r8) :: es ! SVP in Pa 217 ! good down to -100 C 218 es = 10._r8**(-9.09718_r8*(h2otrip/t-1._r8)-3.56654_r8* log10(h2otrip/t)+0.876793_r8*(1._r8-t/h2otrip)+ & 219 log10(6.1071_r8))*100._r8 220 END FUNCTION goffgratch_svp_ice 221 ! Murphy & Koop (2005) 222 223 elemental FUNCTION murphykoop_svp_water(t) RESULT ( es ) 224 REAL(KIND=r8), intent(in) :: t ! Temperature in Kelvin 225 REAL(KIND=r8) :: es ! SVP in Pa 226 ! (good for 123 < T < 332 K) 227 es = exp(54.842763_r8 - (6763.22_r8 / t) - (4.210_r8 * log(t)) + (0.000367_r8 * t) + (tanh(0.0415_r8 * (t - & 228 218.8_r8)) * (53.878_r8 - (1331.22_r8 / t) - (9.44523_r8 * log(t)) + 0.014025_r8 * t))) 229 END FUNCTION murphykoop_svp_water 230 231 elemental FUNCTION murphykoop_svp_ice(t) RESULT ( es ) 232 REAL(KIND=r8), intent(in) :: t ! Temperature in Kelvin 233 REAL(KIND=r8) :: es ! SVP in Pa 234 ! (good down to 110 K) 235 es = exp(9.550426_r8 - (5723.265_r8 / t) + (3.53068_r8 * log(t)) - (0.00728332_r8 * t)) 236 END FUNCTION murphykoop_svp_ice 237 ! Old 1 implementation, also labelled Goff & Gratch (1946) 238 ! The water formula differs only due to compiler-dependent order of 239 ! operations, so differences are roundoff level, usually 0. 240 ! The ice formula gives fairly close answers to the current 241 ! implementation, but has been rearranged, and uses the 242 ! 1 atm melting point of water as the triple point. 243 ! Differences are thus small but above roundoff. 244 ! A curious fact: although using the melting point of water was 245 ! probably a mistake, it mildly improves accuracy for ice svp, 246 ! since it compensates for a systematic error in Goff & Gratch. 247 248 elemental FUNCTION oldgoffgratch_svp_water(t) RESULT ( es ) 249 REAL(KIND=r8), intent(in) :: t 250 REAL(KIND=r8) :: es 251 REAL(KIND=r8) :: ps 252 REAL(KIND=r8) :: e1 253 REAL(KIND=r8) :: e2 254 REAL(KIND=r8) :: f1 255 REAL(KIND=r8) :: f2 256 REAL(KIND=r8) :: f3 257 REAL(KIND=r8) :: f4 258 REAL(KIND=r8) :: f5 259 REAL(KIND=r8) :: f 260 ps = 1013.246_r8 261 e1 = 11.344_r8*(1.0_r8 - t/tboil) 262 e2 = -3.49149_r8*(tboil/t - 1.0_r8) 263 f1 = -7.90298_r8*(tboil/t - 1.0_r8) 264 f2 = 5.02808_r8*log10(tboil/t) 265 f3 = -1.3816_r8*(10.0_r8**e1 - 1.0_r8)/10000000.0_r8 266 f4 = 8.1328_r8*(10.0_r8**e2 - 1.0_r8)/1000.0_r8 267 f5 = log10(ps) 268 f = f1 + f2 + f3 + f4 + f5 269 es = (10.0_r8**f)*100.0_r8 270 END FUNCTION oldgoffgratch_svp_water 271 272 elemental FUNCTION oldgoffgratch_svp_ice(t) RESULT ( es ) 273 REAL(KIND=r8), intent(in) :: t 274 REAL(KIND=r8) :: es 275 REAL(KIND=r8) :: term1 276 REAL(KIND=r8) :: term2 277 REAL(KIND=r8) :: term3 278 term1 = 2.01889049_r8/(tmelt/t) 279 term2 = 3.56654_r8*log(tmelt/t) 280 term3 = 20.947031_r8*(tmelt/t) 281 es = 575.185606e10_r8*exp(-(term1 + term2 + term3)) 282 END FUNCTION oldgoffgratch_svp_ice 283 ! Bolton (1980) 284 ! zm_conv deep convection scheme contained this SVP calculation. 285 ! It appears to be from D. Bolton, 1980, Monthly Weather Review. 286 ! Unlike the other schemes, no distinct ice formula is associated 287 ! with it. (However, a Bolton ice formula exists in CLUBB.) 288 ! The original formula used degrees C, but this function 289 ! takes Kelvin and internally converts. 290 291 elemental FUNCTION bolton_svp_water(t) RESULT ( es ) 292 REAL(KIND=r8), parameter :: c1 = 611.2_r8 293 REAL(KIND=r8), parameter :: c2 = 17.67_r8 294 REAL(KIND=r8), parameter :: c3 = 243.5_r8 295 REAL(KIND=r8), intent(in) :: t ! Temperature in Kelvin 296 REAL(KIND=r8) :: es ! SVP in Pa 297 es = c1*exp((c2*(t - tmelt))/((t - tmelt)+c3)) 298 END FUNCTION bolton_svp_water 299 END MODULE wv_sat_methods 300