1 2! KGEN-generated Fortran source file 3! 4! Filename : mo_rrtm_coeffs.f90 5! Generated at: 2015-02-19 15:30:32 6! KGEN version: 0.4.4 7 8 9 10 MODULE mo_rrtm_coeffs 11 USE mo_kind, ONLY: wp 12 USE mo_rrtm_params, ONLY: preflog 13 USE mo_rrtm_params, ONLY: tref 14 USE rrlw_planck, ONLY: chi_mls 15 IMPLICIT NONE 16 REAL(KIND=wp), parameter :: stpfac = 296._wp/1013._wp 17 CONTAINS 18 19 ! read subroutines 20 ! -------------------------------------------------------------------------------------------- 21 22 SUBROUTINE lrtm_coeffs(kproma, kbdim, klev, play, tlay, coldry, wkl, wbroad, laytrop, jp, jt, jt1, colh2o, colco2, colo3, & 23 coln2o, colco, colch4, colo2, colbrd, fac00, fac01, fac10, fac11, rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, & 24 rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, selffac, selffrac, & 25 indself, forfac, forfrac, indfor, minorfrac, scaleminor, scaleminorn2, indminor) 26 INTEGER, intent(in) :: kbdim 27 INTEGER, intent(in) :: klev 28 INTEGER, intent(in) :: kproma 29 ! number of columns 30 ! maximum number of column as first dim is declared in calling (sub)prog. 31 ! total number of layers 32 REAL(KIND=wp), intent(in) :: wkl(:,:,:) 33 REAL(KIND=wp), intent(in) :: play(kbdim,klev) 34 REAL(KIND=wp), intent(in) :: tlay(kbdim,klev) 35 REAL(KIND=wp), intent(in) :: coldry(kbdim,klev) 36 REAL(KIND=wp), intent(in) :: wbroad(kbdim,klev) 37 ! layer pressures (mb) 38 ! layer temperatures (K) 39 ! dry air column density (mol/cm2) 40 ! broadening gas column density (mol/cm2) 41 !< molecular amounts (mol/cm-2) (mxmol,klev) 42 ! 43 ! Output Dimensions kproma, klev unless otherwise specified 44 ! 45 INTEGER, intent(out) :: laytrop(kbdim) 46 INTEGER, intent(out) :: jp(kbdim,klev) 47 INTEGER, intent(out) :: jt(kbdim,klev) 48 INTEGER, intent(out) :: jt1(kbdim,klev) 49 INTEGER, intent(out) :: indfor(kbdim,klev) 50 INTEGER, intent(out) :: indself(kbdim,klev) 51 INTEGER, intent(out) :: indminor(kbdim,klev) 52 !< tropopause layer index 53 ! 54 ! 55 ! 56 ! 57 ! 58 ! 59 REAL(KIND=wp), intent(out) :: forfac(kbdim,klev) 60 REAL(KIND=wp), intent(out) :: colch4(kbdim,klev) 61 REAL(KIND=wp), intent(out) :: colco2(kbdim,klev) 62 REAL(KIND=wp), intent(out) :: colh2o(kbdim,klev) 63 REAL(KIND=wp), intent(out) :: coln2o(kbdim,klev) 64 REAL(KIND=wp), intent(out) :: forfrac(kbdim,klev) 65 REAL(KIND=wp), intent(out) :: selffrac(kbdim,klev) 66 REAL(KIND=wp), intent(out) :: colo3(kbdim,klev) 67 REAL(KIND=wp), intent(out) :: fac00(kbdim,klev) 68 REAL(KIND=wp), intent(out) :: colo2(kbdim,klev) 69 REAL(KIND=wp), intent(out) :: fac01(kbdim,klev) 70 REAL(KIND=wp), intent(out) :: fac10(kbdim,klev) 71 REAL(KIND=wp), intent(out) :: fac11(kbdim,klev) 72 REAL(KIND=wp), intent(out) :: selffac(kbdim,klev) 73 REAL(KIND=wp), intent(out) :: colbrd(kbdim,klev) 74 REAL(KIND=wp), intent(out) :: colco(kbdim,klev) 75 REAL(KIND=wp), intent(out) :: rat_h2oco2(kbdim,klev) 76 REAL(KIND=wp), intent(out) :: rat_h2oco2_1(kbdim,klev) 77 REAL(KIND=wp), intent(out) :: rat_h2oo3(kbdim,klev) 78 REAL(KIND=wp), intent(out) :: rat_h2oo3_1(kbdim,klev) 79 REAL(KIND=wp), intent(out) :: rat_h2on2o(kbdim,klev) 80 REAL(KIND=wp), intent(out) :: rat_h2on2o_1(kbdim,klev) 81 REAL(KIND=wp), intent(out) :: rat_h2och4(kbdim,klev) 82 REAL(KIND=wp), intent(out) :: rat_h2och4_1(kbdim,klev) 83 REAL(KIND=wp), intent(out) :: rat_n2oco2(kbdim,klev) 84 REAL(KIND=wp), intent(out) :: rat_n2oco2_1(kbdim,klev) 85 REAL(KIND=wp), intent(out) :: rat_o3co2(kbdim,klev) 86 REAL(KIND=wp), intent(out) :: rat_o3co2_1(kbdim,klev) 87 REAL(KIND=wp), intent(out) :: scaleminor(kbdim,klev) 88 REAL(KIND=wp), intent(out) :: scaleminorn2(kbdim,klev) 89 REAL(KIND=wp), intent(out) :: minorfrac(kbdim,klev) 90 !< column amount (h2o) 91 !< column amount (co2) 92 !< column amount (o3) 93 !< column amount (n2o) 94 !< column amount (co) 95 !< column amount (ch4) 96 !< column amount (o2) 97 !< column amount (broadening gases) 98 !< 99 !< 100 !< 101 !< 102 !< 103 INTEGER :: jk 104 REAL(KIND=wp) :: colmol(kbdim,klev) 105 REAL(KIND=wp) :: factor(kbdim,klev) 106 ! ------------------------------------------------ 107 CALL srtm_coeffs(kproma, kbdim, klev, play, tlay, coldry, wkl, laytrop, jp, jt, jt1, colch4, colco2, colh2o, colmol, & 108 coln2o, colo2, colo3, fac00, fac01, fac10, fac11, selffac, selffrac, indself, forfac, forfrac, indfor) 109 colbrd(1:kproma,1:klev) = 1.e-20_wp * wbroad(1:kproma,1:klev) 110 colco(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,5,1:klev), 1.e-32_wp * & 111 coldry(1:kproma,1:klev), wkl(1:kproma,5,1:klev) > 0._wp) 112 ! 113 ! Water vapor continuum broadening factors are used differently in LW and SW? 114 ! 115 forfac(1:kproma,1:klev) = forfac(1:kproma,1:klev) * colh2o(1:kproma,1:klev) 116 selffac(1:kproma,1:klev) = selffac(1:kproma,1:klev) * colh2o(1:kproma,1:klev) 117 ! 118 ! Setup reference ratio to be used in calculation of binary species parameter. 119 ! 120 DO jk = 1, klev 121 rat_h2oco2(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk))/chi_mls(2,jp(1:kproma, jk)) 122 rat_h2oco2_1(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk)+1)/chi_mls(2,jp(1:kproma, jk)+1) 123 ! 124 ! Needed only in lower atmos (plog > 4.56_wp) 125 ! 126 rat_h2oo3(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk))/chi_mls(3,jp(1:kproma, jk)) 127 rat_h2oo3_1(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk)+1)/chi_mls(3,jp(1:kproma, jk)+1) 128 rat_h2on2o(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk))/chi_mls(4,jp(1:kproma, jk)) 129 rat_h2on2o_1(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk)+1)/chi_mls(4,jp(1:kproma, jk)+1) 130 rat_h2och4(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk))/chi_mls(6,jp(1:kproma, jk)) 131 rat_h2och4_1(1:kproma, jk) = chi_mls(1,jp(1:kproma, jk)+1)/chi_mls(6,jp(1:kproma, jk)+1) 132 rat_n2oco2(1:kproma, jk) = chi_mls(4,jp(1:kproma, jk))/chi_mls(2,jp(1:kproma, jk)) 133 rat_n2oco2_1(1:kproma, jk) = chi_mls(4,jp(1:kproma, jk)+1)/chi_mls(2,jp(1:kproma, jk)+1) 134 ! 135 ! Needed only in upper atmos (plog <= 4.56_wp) 136 ! 137 rat_o3co2(1:kproma, jk) = chi_mls(3,jp(1:kproma, jk))/chi_mls(2,jp(1:kproma, jk)) 138 rat_o3co2_1(1:kproma, jk) = chi_mls(3,jp(1:kproma, jk)+1)/chi_mls(2,jp(1:kproma, jk)+1) 139 END DO 140 ! 141 ! Set up factors needed to separately include the minor gases 142 ! in the calculation of absorption coefficient 143 ! 144 scaleminor(1:kproma,1:klev) = play(1:kproma,1:klev)/tlay(1:kproma,1:klev) 145 scaleminorn2(1:kproma,1:klev) = scaleminor(1:kproma,1:klev) * (wbroad(1:kproma,1:klev)/(& 146 coldry(1:kproma,1:klev)+wkl(1:kproma,1,1:klev))) 147 factor(1:kproma,1:klev) = (tlay(1:kproma,1:klev)-180.8_wp)/7.2_wp 148 indminor(1:kproma,1:klev) = min(18, max(1, int(factor(1:kproma,1:klev)))) 149 minorfrac(1:kproma,1:klev) = (tlay(1:kproma,1:klev)-180.8_wp)/7.2_wp - float(indminor(1:kproma,1:klev)) 150 END SUBROUTINE lrtm_coeffs 151 ! -------------------------------------------------------------------------------------------- 152 153 SUBROUTINE srtm_coeffs(kproma, kbdim, klev, play, tlay, coldry, wkl, laytrop, jp, jt, jt1, colch4, colco2, colh2o, colmol,& 154 coln2o, colo2, colo3, fac00, fac01, fac10, fac11, selffac, selffrac, indself, forfac, forfrac, indfor) 155 INTEGER, intent(in) :: kbdim 156 INTEGER, intent(in) :: klev 157 INTEGER, intent(in) :: kproma 158 ! number of columns 159 ! maximum number of col. as declared in calling (sub)programs 160 ! total number of layers 161 REAL(KIND=wp), intent(in) :: play(kbdim,klev) 162 REAL(KIND=wp), intent(in) :: tlay(kbdim,klev) 163 REAL(KIND=wp), intent(in) :: wkl(:,:,:) 164 REAL(KIND=wp), intent(in) :: coldry(kbdim,klev) 165 ! layer pressures (mb) 166 ! layer temperatures (K) 167 ! dry air column density (mol/cm2) 168 !< molecular amounts (mol/cm-2) (mxmol,klev) 169 ! 170 ! Output Dimensions kproma, klev unless otherwise specified 171 ! 172 INTEGER, intent(out) :: jp(kbdim,klev) 173 INTEGER, intent(out) :: jt(kbdim,klev) 174 INTEGER, intent(out) :: jt1(kbdim,klev) 175 INTEGER, intent(out) :: laytrop(kbdim) 176 INTEGER, intent(out) :: indfor(kbdim,klev) 177 INTEGER, intent(out) :: indself(kbdim,klev) 178 !< tropopause layer index 179 ! 180 ! 181 ! 182 ! 183 REAL(KIND=wp), intent(out) :: fac10(kbdim,klev) 184 REAL(KIND=wp), intent(out) :: fac00(kbdim,klev) 185 REAL(KIND=wp), intent(out) :: fac11(kbdim,klev) 186 REAL(KIND=wp), intent(out) :: fac01(kbdim,klev) 187 REAL(KIND=wp), intent(out) :: colh2o(kbdim,klev) 188 REAL(KIND=wp), intent(out) :: colco2(kbdim,klev) 189 REAL(KIND=wp), intent(out) :: colo3(kbdim,klev) 190 REAL(KIND=wp), intent(out) :: coln2o(kbdim,klev) 191 REAL(KIND=wp), intent(out) :: colch4(kbdim,klev) 192 REAL(KIND=wp), intent(out) :: colo2(kbdim,klev) 193 REAL(KIND=wp), intent(out) :: colmol(kbdim,klev) 194 REAL(KIND=wp), intent(out) :: forfac(kbdim,klev) 195 REAL(KIND=wp), intent(out) :: selffac(kbdim,klev) 196 REAL(KIND=wp), intent(out) :: forfrac(kbdim,klev) 197 REAL(KIND=wp), intent(out) :: selffrac(kbdim,klev) 198 !< column amount (h2o) 199 !< column amount (co2) 200 !< column amount (o3) 201 !< column amount (n2o) 202 !< column amount (ch4) 203 !< column amount (o2) 204 !< 205 !< 206 !< 207 !< 208 !< 209 INTEGER :: jp1(kbdim,klev) 210 INTEGER :: jk 211 REAL(KIND=wp) :: plog (kbdim,klev) 212 REAL(KIND=wp) :: fp (kbdim,klev) 213 REAL(KIND=wp) :: ft (kbdim,klev) 214 REAL(KIND=wp) :: ft1 (kbdim,klev) 215 REAL(KIND=wp) :: water (kbdim,klev) 216 REAL(KIND=wp) :: scalefac(kbdim,klev) 217 REAL(KIND=wp) :: compfp(kbdim,klev) 218 REAL(KIND=wp) :: factor (kbdim,klev) 219 ! ------------------------------------------------------------------------- 220 ! 221 ! Find the two reference pressures on either side of the 222 ! layer pressure. Store them in JP and JP1. Store in FP the 223 ! fraction of the difference (in ln(pressure)) between these 224 ! two values that the layer pressure lies. 225 ! 226 plog(1:kproma,1:klev) = log(play(1:kproma,1:klev)) 227 jp(1:kproma,1:klev) = min(58,max(1,int(36._wp - 5*(plog(1:kproma,1:klev)+0.04_wp)))) 228 jp1(1:kproma,1:klev) = jp(1:kproma,1:klev) + 1 229 DO jk = 1, klev 230 fp(1:kproma,jk) = 5._wp *(preflog(jp(1:kproma,jk)) - plog(1:kproma,jk)) 231 END DO 232 ! 233 ! Determine, for each reference pressure (JP and JP1), which 234 ! reference temperature (these are different for each 235 ! reference pressure) is nearest the layer temperature but does 236 ! not exceed it. Store these indices in JT and JT1, resp. 237 ! Store in FT (resp. FT1) the fraction of the way between JT 238 ! (JT1) and the next highest reference temperature that the 239 ! layer temperature falls. 240 ! 241 DO jk = 1, klev 242 jt(1:kproma,jk) = min(4,max(1,int(3._wp + (tlay(1:kproma,jk) - tref(& 243 jp (1:kproma,jk)))/15._wp))) 244 jt1(1:kproma,jk) = min(4,max(1,int(3._wp + (tlay(1:kproma,jk) - & 245 tref(jp1(1:kproma,jk)))/15._wp))) 246 END DO 247 DO jk = 1, klev 248 ft(1:kproma,jk) = ((tlay(1:kproma,jk)-tref(jp (1:kproma,jk)))/15._wp) - float(jt (& 249 1:kproma,jk)-3) 250 ft1(1:kproma,jk) = ((tlay(1:kproma,jk)-tref(jp1(1:kproma,jk)))/15._wp) - float(jt1(& 251 1:kproma,jk)-3) 252 END DO 253 water(1:kproma,1:klev) = wkl(1:kproma,1,1:klev)/coldry(1:kproma,1:klev) 254 scalefac(1:kproma,1:klev) = play(1:kproma,1:klev) * stpfac / tlay(1:kproma,1:klev) 255 ! 256 ! We have now isolated the layer ln pressure and temperature, 257 ! between two reference pressures and two reference temperatures 258 ! (for each reference pressure). We multiply the pressure 259 ! fraction FP with the appropriate temperature fractions to get 260 ! the factors that will be needed for the interpolation that yields 261 ! the optical depths (performed in routines TAUGBn for band n).` 262 ! 263 compfp(1:kproma,1:klev) = 1. - fp(1:kproma,1:klev) 264 fac10(1:kproma,1:klev) = compfp(1:kproma,1:klev) * ft(1:kproma,1:klev) 265 fac00(1:kproma,1:klev) = compfp(1:kproma,1:klev) * (1._wp - ft(1:kproma,1:klev)) 266 fac11(1:kproma,1:klev) = fp(1:kproma,1:klev) * ft1(1:kproma,1:klev) 267 fac01(1:kproma,1:klev) = fp(1:kproma,1:klev) * (1._wp - ft1(1:kproma,1:klev)) 268 ! Tropopause defined in terms of pressure (~100 hPa) 269 ! We're looking for the first layer (counted from the bottom) at which the pressure reaches 270 ! or falls below this value 271 ! 272 laytrop(1:kproma) = count(plog(1:kproma,1:klev) > 4.56_wp, dim = 2) 273 ! 274 ! Calculate needed column amounts. 275 ! Only a few ratios are used in the upper atmosphere but masking may be less efficient 276 ! 277 colh2o(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,1,1:klev), 1.e-32_wp * & 278 coldry(1:kproma,1:klev), wkl(1:kproma,1,1:klev) > 0._wp) 279 colco2(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,2,1:klev), 1.e-32_wp * & 280 coldry(1:kproma,1:klev), wkl(1:kproma,2,1:klev) > 0._wp) 281 colo3(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,3,1:klev), 1.e-32_wp * & 282 coldry(1:kproma,1:klev), wkl(1:kproma,3,1:klev) > 0._wp) 283 coln2o(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,4,1:klev), 1.e-32_wp * & 284 coldry(1:kproma,1:klev), wkl(1:kproma,4,1:klev) > 0._wp) 285 colch4(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,6,1:klev), 1.e-32_wp * & 286 coldry(1:kproma,1:klev), wkl(1:kproma,6,1:klev) > 0._wp) 287 colo2(1:kproma,1:klev) = merge(1.e-20_wp * wkl(1:kproma,7,1:klev), 1.e-32_wp * & 288 coldry(1:kproma,1:klev), wkl(1:kproma,7,1:klev) > 0._wp) 289 colmol(1:kproma,1:klev) = 1.e-20_wp * coldry(1:kproma,1:klev) + colh2o(1:kproma,1:klev) 290 ! ------------------------------------------ 291 ! Interpolation coefficients 292 ! 293 forfac(1:kproma,1:klev) = scalefac(1:kproma,1:klev) / (1._wp+water(1:kproma,1:klev)) 294 ! 295 ! Set up factors needed to separately include the water vapor 296 ! self-continuum in the calculation of absorption coefficient. 297 ! 298 selffac(1:kproma,1:klev) = water(1:kproma,1:klev) * forfac(1:kproma,1:klev) 299 ! 300 ! If the pressure is less than ~100mb, perform a different set of species 301 ! interpolations. 302 ! 303 factor(1:kproma,1:klev) = (332.0_wp-tlay(1:kproma,1:klev))/36.0_wp 304 indfor(1:kproma,1:klev) = merge(3, min(2, max(1, int(factor(1:kproma,& 305 1:klev)))), plog(1:kproma,1:klev) <= 4.56_wp) 306 forfrac(1:kproma,1:klev) = merge((tlay(1:kproma,1:klev)-188.0_wp)/36.0_wp - 1.0_wp, factor(1:kproma,& 307 1:klev) - float(indfor(1:kproma,1:klev)), plog(1:kproma,1:klev) <= 4.56_wp) 308 ! In RRTMG code, this calculation is done only in the lower atmosphere (plog > 4.56) 309 ! 310 factor(1:kproma,1:klev) = (tlay(1:kproma,1:klev)-188.0_wp)/7.2_wp 311 indself(1:kproma,1:klev) = min(9, max(1, int(factor(1:kproma,1:klev))-7)) 312 selffrac(1:kproma,1:klev) = factor(1:kproma,1:klev) - float(indself(1:kproma,1:klev) + 7) 313 END SUBROUTINE srtm_coeffs 314 END MODULE mo_rrtm_coeffs 315