1 2! KGEN-generated Fortran source file 3! 4! Filename : mo_lrtm_driver.f90 5! Generated at: 2015-02-19 15:30:29 6! KGEN version: 0.4.4 7 8 9 10 MODULE mo_lrtm_driver 11 USE mo_kind, ONLY: wp 12 USE mo_physical_constants, ONLY: amw 13 USE mo_physical_constants, ONLY: amd 14 USE mo_physical_constants, ONLY: grav 15 USE mo_rrtm_params, ONLY: nbndlw 16 USE mo_rrtm_params, ONLY: ngptlw 17 USE mo_radiation_parameters, ONLY: do_gpoint 18 USE mo_radiation_parameters, ONLY: i_overlap 19 USE mo_radiation_parameters, ONLY: l_do_sep_clear_sky 20 USE mo_radiation_parameters, ONLY: rad_undef 21 USE mo_lrtm_setup, ONLY: ngb 22 USE mo_lrtm_setup, ONLY: ngs 23 USE mo_lrtm_setup, ONLY: delwave 24 USE mo_lrtm_setup, ONLY: nspa 25 USE mo_lrtm_setup, ONLY: nspb 26 USE rrlw_planck, ONLY: totplanck 27 USE mo_rrtm_coeffs, ONLY: lrtm_coeffs 28 USE mo_lrtm_gas_optics, ONLY: gas_optics_lw 29 USE mo_lrtm_solver, ONLY: find_secdiff 30 USE mo_lrtm_solver, ONLY: lrtm_solver 31 USE mo_cld_sampling, ONLY: sample_cld_state 32 USE mo_spec_sampling, ONLY: spec_sampling_strategy 33 USE mo_spec_sampling, ONLY: get_gpoint_set 34 USE mo_taumol03, ONLY: taumol03_lwr,taumol03_upr 35 USE mo_taumol04, ONLY: taumol04_lwr,taumol04_upr 36 USE rrlw_planck, ONLY: chi_mls 37 USE rrlw_kg03, ONLY: selfref 38 USE rrlw_kg03, ONLY: forref 39 USE rrlw_kg03, ONLY: ka_mn2o 40 USE rrlw_kg03, ONLY: absa 41 USE rrlw_kg03, ONLY: fracrefa 42 USE rrlw_kg03, ONLY: kb_mn2o 43 USE rrlw_kg03, ONLY: absb 44 USE rrlw_kg03, ONLY: fracrefb 45 IMPLICIT NONE 46 PRIVATE 47 PUBLIC lrtm 48 CONTAINS 49 50 ! read subroutines 51 !----------------------------------------------------------------------------- 52 !> 53 !! @brief Prepares information for radiation call 54 !! 55 !! @remarks: This program is the driver subroutine for the longwave radiative 56 !! transfer routine. This routine is adapted from the AER LW RRTMG_LW model 57 !! that itself has been adapted from RRTM_LW for improved efficiency. Our 58 !! routine does the spectral integration externally (the solver is explicitly 59 !! called for each g-point, so as to facilitate sampling of g-points 60 !! This routine: 61 !! 1) calls INATM to read in the atmospheric profile from GCM; 62 !! all layering in RRTMG is ordered from surface to toa. 63 !! 2) calls COEFFS to calculate various quantities needed for 64 !! the radiative transfer algorithm. This routine is called only once for 65 !! any given thermodynamic state, i.e., it does not change if clouds chanege 66 !! 3) calls TAUMOL to calculate gaseous optical depths for each 67 !! of the 16 spectral bands, this is updated band by band. 68 !! 4) calls SOLVER (for both clear and cloudy profiles) to perform the 69 !! radiative transfer calculation with a maximum-random cloud 70 !! overlap method, or calls RTRN to use random cloud overlap. 71 !! 5) passes the necessary fluxes and cooling rates back to GCM 72 !! 73 ! 74 75 SUBROUTINE lrtm(kproma, kbdim, klev, play, psfc, tlay, tlev, tsfc, wkl, wx, coldry, emis, cldfr, taucld, tauaer, rnseeds, & 76 strategy, n_gpts_ts, uflx, dflx, uflxc, dflxc) 77 INTEGER, intent(in) :: klev 78 INTEGER, intent(in) :: kbdim 79 INTEGER, intent(in) :: kproma 80 !< Maximum block length 81 !< Number of horizontal columns 82 !< Number of model layers 83 REAL(KIND=wp), intent(in) :: wx(:,:,:) 84 REAL(KIND=wp), intent(in) :: cldfr(kbdim,klev) 85 REAL(KIND=wp), intent(in) :: taucld(kbdim,klev,nbndlw) 86 REAL(KIND=wp), intent(in) :: wkl(:,:,:) 87 REAL(KIND=wp), intent(in) :: coldry(kbdim,klev) 88 REAL(KIND=wp), intent(in) :: play(kbdim,klev) 89 REAL(KIND=wp), intent(in) :: tlay(kbdim,klev) 90 REAL(KIND=wp), intent(in) :: tauaer(kbdim,klev,nbndlw) 91 REAL(KIND=wp), intent(in) :: tlev(kbdim,klev+1) 92 REAL(KIND=wp), intent(in) :: tsfc(kbdim) 93 REAL(KIND=wp), intent(in) :: psfc(kbdim) 94 REAL(KIND=wp), intent(in) :: emis(kbdim,nbndlw) 95 !< Layer pressures [hPa, mb] (kbdim,klev) 96 !< Surface pressure [hPa, mb] (kbdim) 97 !< Layer temperatures [K] (kbdim,klev) 98 !< Interface temperatures [K] (kbdim,klev+1) 99 !< Surface temperature [K] (kbdim) 100 !< Gas volume mixing ratios 101 !< CFC type gas volume mixing ratios 102 !< Column dry amount 103 !< Surface emissivity (kbdim,nbndlw) 104 !< Cloud fraction (kbdim,klev) 105 !< Coud optical depth (kbdim,klev,nbndlw) 106 !< Aerosol optical depth (kbdim,klev,nbndlw) 107 ! Variables for sampling cloud state and spectral points 108 INTEGER, intent(inout) :: rnseeds(:, :) !< Seeds for random number generator (kbdim,:) 109 TYPE(spec_sampling_strategy), intent(in) :: strategy 110 INTEGER, intent(in ) :: n_gpts_ts 111 REAL(KIND=wp), intent(out) :: uflx (kbdim,0:klev) 112 REAL(KIND=wp), intent(out) :: dflx (kbdim,0:klev) 113 REAL(KIND=wp), intent(out) :: uflxc(kbdim,0:klev) 114 REAL(KIND=wp), intent(out) :: dflxc(kbdim,0:klev) 115 !< Tot sky longwave upward flux [W/m2], (kbdim,0:klev) 116 !< Tot sky longwave downward flux [W/m2], (kbdim,0:klev) 117 !< Clr sky longwave upward flux [W/m2], (kbdim,0:klev) 118 !< Clr sky longwave downward flux [W/m2], (kbdim,0:klev) 119 REAL(KIND=wp) :: taug(klev) !< Properties for one column at a time >! gas optical depth 120 REAL(KIND=wp) :: rrpk_taug03(kbdim,klev) 121 REAL(KIND=wp) :: rrpk_taug04(kbdim,klev) 122 REAL(KIND=wp) :: fracs(kbdim,klev,n_gpts_ts) 123 REAL(KIND=wp) :: taut (kbdim,klev,n_gpts_ts) 124 REAL(KIND=wp) :: tautot(kbdim,klev,n_gpts_ts) 125 REAL(KIND=wp) :: pwvcm(kbdim) 126 REAL(KIND=wp) :: secdiff(kbdim) 127 !< Planck fraction per g-point 128 !< precipitable water vapor [cm] 129 !< diffusivity angle for RT calculation 130 !< gaseous + aerosol optical depths for all columns 131 !< cloud + gaseous + aerosol optical depths for all columns 132 REAL(KIND=wp) :: planklay(kbdim, klev,nbndlw) 133 REAL(KIND=wp) :: planklev(kbdim,0:klev,nbndlw) 134 REAL(KIND=wp) :: plankbnd(kbdim, nbndlw) ! Properties for all bands 135 ! Planck function at mid-layer 136 ! Planck function at level interfaces 137 ! Planck function at surface 138 REAL(KIND=wp) :: layplnk(kbdim, klev) 139 REAL(KIND=wp) :: levplnk(kbdim,0:klev) 140 REAL(KIND=wp) :: bndplnk(kbdim) 141 REAL(KIND=wp) :: srfemis(kbdim) ! Properties for a single set of columns/g-points 142 ! Planck function at mid-layer 143 ! Planck function at level interfaces 144 ! Planck function at surface 145 ! Surface emission 146 REAL(KIND=wp) :: zgpfd(kbdim,0:klev) 147 REAL(KIND=wp) :: zgpfu(kbdim,0:klev) 148 REAL(KIND=wp) :: zgpcu(kbdim,0:klev) 149 REAL(KIND=wp) :: zgpcd(kbdim,0:klev) 150 ! < gpoint clearsky downward flux 151 ! < gpoint clearsky downward flux 152 ! < gpoint fullsky downward flux 153 ! < gpoint fullsky downward flux 154 ! ----------------- 155 ! Variables for gas optics calculations 156 INTEGER :: jt1 (kbdim,klev) 157 INTEGER :: indfor (kbdim,klev) 158 INTEGER :: indself (kbdim,klev) 159 INTEGER :: indminor(kbdim,klev) 160 INTEGER :: laytrop (kbdim ) 161 INTEGER :: jp (kbdim,klev) 162 INTEGER :: rrpk_jp (klev,kbdim) 163 INTEGER :: jt (kbdim,klev) 164 INTEGER :: rrpk_jt (kbdim,0:1,klev) 165 !< tropopause layer index 166 !< lookup table index 167 !< lookup table index 168 !< lookup table index 169 REAL(KIND=wp) :: wbrodl (kbdim,klev) 170 REAL(KIND=wp) :: selffac (kbdim,klev) 171 REAL(KIND=wp) :: colh2o (kbdim,klev) 172 REAL(KIND=wp) :: colo3 (kbdim,klev) 173 REAL(KIND=wp) :: coln2o (kbdim,klev) 174 REAL(KIND=wp) :: colco (kbdim,klev) 175 REAL(KIND=wp) :: selffrac (kbdim,klev) 176 REAL(KIND=wp) :: colch4 (kbdim,klev) 177 REAL(KIND=wp) :: colo2 (kbdim,klev) 178 REAL(KIND=wp) :: colbrd (kbdim,klev) 179 REAL(KIND=wp) :: minorfrac (kbdim,klev) 180 REAL(KIND=wp) :: scaleminorn2(kbdim,klev) 181 REAL(KIND=wp) :: scaleminor (kbdim,klev) 182 REAL(KIND=wp) :: forfac (kbdim,klev) 183 REAL(KIND=wp) :: colco2 (kbdim,klev) 184 REAL(KIND=wp) :: forfrac (kbdim,klev) 185 !< column amount (h2o) 186 !< column amount (co2) 187 !< column amount (o3) 188 !< column amount (n2o) 189 !< column amount (co) 190 !< column amount (ch4) 191 !< column amount (o2) 192 !< column amount (broadening gases) 193 REAL(KIND=wp) :: wx_loc(size(wx, 2), size(wx, 3)) 194 !< Normalized CFC amounts (molecules/cm^2) 195 REAL(KIND=wp) :: fac00(kbdim,klev) 196 REAL(KIND=wp) :: fac01(kbdim,klev) 197 REAL(KIND=wp) :: fac10(kbdim,klev) 198 REAL(KIND=wp) :: fac11(kbdim,klev) 199 REAL(KIND=wp) :: rrpk_fac0(kbdim,0:1,klev) 200 REAL(KIND=wp) :: rrpk_fac1(kbdim,0:1,klev) 201 REAL(KIND=wp) :: rat_n2oco2 (kbdim,klev) 202 REAL(KIND=wp) :: rat_o3co2 (kbdim,klev) 203 REAL(KIND=wp) :: rat_h2on2o (kbdim,klev) 204 REAL(KIND=wp) :: rat_n2oco2_1(kbdim,klev) 205 REAL(KIND=wp) :: rat_h2on2o_1(kbdim,klev) 206 REAL(KIND=wp) :: rat_h2oco2_1(kbdim,klev) 207 REAL(KIND=wp) :: rat_h2oo3 (kbdim,klev) 208 REAL(KIND=wp) :: rat_h2och4 (kbdim,klev) 209 REAL(KIND=wp) :: rat_h2oco2 (kbdim,klev) 210 REAL(KIND=wp) :: rrpk_rat_h2oco2 (kbdim,0:1,klev) 211 REAL(KIND=wp) :: rrpk_rat_o3co2 (kbdim,0:1,klev) 212 REAL(KIND=wp) :: rat_h2oo3_1 (kbdim,klev) 213 REAL(KIND=wp) :: rat_o3co2_1 (kbdim,klev) 214 REAL(KIND=wp) :: rat_h2och4_1(kbdim,klev) 215 ! ----------------- 216 INTEGER :: jl,jlBegin,simdStep=96 217 INTEGER :: ig 218 INTEGER :: jk ! loop indicies 219 INTEGER :: igs(kbdim, n_gpts_ts) 220 INTEGER :: ibs(kbdim, n_gpts_ts) 221 INTEGER :: ib 222 INTEGER :: igpt 223 INTEGER*8 :: start_clock,stop_clock,rate_clock 224 REAL :: overall_time=0 225 ! minimum val for clouds 226 ! Variables for sampling strategy 227 REAL(KIND=wp) :: gpt_scaling 228 REAL(KIND=wp) :: clrsky_scaling(1:kbdim) 229 REAL(KIND=wp) :: smp_tau(kbdim, klev, n_gpts_ts) 230 LOGICAL :: cldmask(kbdim, klev, n_gpts_ts) 231 LOGICAL :: colcldmask(kbdim, n_gpts_ts) !< cloud mask in each cell 232 !< cloud mask for each column 233 ! 234 ! -------------------------------- 235 ! 236 ! 1.0 Choose a set of g-points to do consistent with the spectral sampling strategy 237 ! 238 ! -------------------------------- 239 gpt_scaling = real(ngptlw,kind=wp)/real(n_gpts_ts,kind=wp) 240 ! Standalone logic 241 IF (do_gpoint == 0) THEN 242 igs(1:kproma,1:n_gpts_ts) = get_gpoint_set(kproma, kbdim, strategy, rnseeds) 243 ELSE IF (n_gpts_ts == 1) THEN ! Standalone logic 244 IF (do_gpoint > ngptlw) RETURN 245 igs(:, 1:n_gpts_ts) = do_gpoint 246 ELSE 247 PRINT *, "Asking for gpoint fluxes for too many gpoints!" 248 STOP 249 END IF 250 ! Save the band nunber associated with each gpoint 251 DO jl = 1, kproma 252 DO ig = 1, n_gpts_ts 253 ibs(jl, ig) = ngb(igs(jl, ig)) 254 END DO 255 END DO 256 ! 257 ! --- 2.0 Optical properties 258 ! 259 ! --- 2.1 Cloud optical properties. 260 ! -------------------------------- 261 ! Cloud optical depth is only saved for the band associated with this g-point 262 ! We sample clouds first because we may want to adjust water vapor based 263 ! on presence/absence of clouds 264 ! 265 CALL sample_cld_state(kproma, kbdim, klev, n_gpts_ts, rnseeds(:,:), i_overlap, cldfr(:,:), cldmask(:,:,:)) 266 !IBM* ASSERT(NODEPS) 267 DO ig = 1, n_gpts_ts 268 DO jl = 1, kproma 269 smp_tau(jl,:,ig) = merge(taucld(jl,1:klev,ibs(jl,ig)), 0._wp, cldmask(jl,:,ig)) 270 END DO 271 END DO ! Loop over samples - done with cloud optical depth calculations 272 ! 273 ! Cloud masks for sorting out clear skies - by cell and by column 274 ! 275 IF (.not. l_do_sep_clear_sky) THEN 276 ! 277 ! Are any layers cloudy? 278 ! 279 colcldmask(1:kproma, 1:n_gpts_ts) = any(cldmask(1:kproma,1:klev,1:n_gpts_ts), dim=2) 280 ! 281 ! Clear-sky scaling is gpt_scaling/frac_clr or 0 if all samples are cloudy 282 ! 283 clrsky_scaling(1:kproma) = gpt_scaling * & 284 merge(real(n_gpts_ts,kind=wp) / (real(n_gpts_ts - count(& 285 colcldmask(1:kproma,:),dim=2),kind=wp)), & 286 0._wp, any(.not. colcldmask(1:kproma,:),dim=2)) 287 END IF 288 ! 289 ! --- 2.2. Gas optical depth calculations 290 ! 291 ! -------------------------------- 292 ! 293 ! 2.2.1 Calculate information needed by the radiative transfer routine 294 ! that is specific to this atmosphere, especially some of the 295 ! coefficients and indices needed to compute the optical depths 296 ! by interpolating data from stored reference atmospheres. 297 ! The coefficients are functions of temperature and pressure and remain the same 298 ! for all g-point samples. 299 ! If gas concentrations, temperatures, or pressures vary with sample (ig) 300 ! the coefficients need to be calculated inside the loop over samples 301 ! -------------------------------- 302 ! 303 ! Broadening gases -- the number of molecules per cm^2 of all gases not specified explicitly 304 ! (water is excluded) 305 wbrodl(1:kproma,1:klev) = coldry(1:kproma,1:klev) - sum(wkl(1:kproma,2:,1:klev), dim=2) 306 CALL lrtm_coeffs(kproma, kbdim, klev, play, tlay, coldry, wkl, wbrodl, laytrop, jp, jt, jt1, colh2o, colco2, colo3, & 307 coln2o, colco, colch4, colo2, colbrd, fac00, fac01, fac10, fac11, rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, & 308 rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, selffac, & 309 selffrac, indself, forfac, forfrac, indfor, minorfrac, scaleminor, scaleminorn2, indminor) 310 ! 311 ! 2.2.2 Loop over g-points calculating gas optical properties. 312 ! 313 ! -------------------------------- 314 !IBM* ASSERT(NODEPS) 315 !CALL system_clock(start_clock,rate_clock) 316 rrpk_rat_h2oco2(:,0,:) = rat_h2oco2 317 rrpk_rat_h2oco2(:,1,:) = (rat_h2oco2_1) 318 rrpk_rat_o3co2(:,0,:) = rat_o3co2 319 rrpk_rat_o3co2(:,1,:) = (rat_o3co2_1) 320 rrpk_fac0(:,0,:) = fac00 321 rrpk_fac0(:,1,:) = fac01 322 rrpk_fac1(:,0,:) = fac10 323 rrpk_fac1(:,1,:) = fac11 324 rrpk_jt(:,0,:) = jt 325 rrpk_jt(:,1,:) = jt1 326 !CALL system_clock(stop_clock,rate_clock) 327 !overall_time=overall_time + (stop_clock-start_clock)/REAL(rate_clock) 328 !print *,n_gpts_ts 329 !print *,"===",kproma 330 DO ig = 1, n_gpts_ts 331 igpt=igs(1,ig) 332 IF(ngb(igpt) == 3) Then 333 CALL system_clock(start_clock, rate_clock) 334 jl=kproma 335 DO jlBegin = 1,kproma,simdStep 336 jl = jlBegin+simdStep-1 337 call taumol03_lwr(jl,jlBegin,laytrop(1), klev, & 338 rrpk_rat_h2oco2(jlBegin:jl,:,:), colco2(jlBegin:jl,:), colh2o(jlBegin:jl,:), coln2o(jlBegin:jl,:), coldry(jlBegin:jl,:), & 339 rrpk_fac0(jlBegin:jl,:,:), rrpk_fac1(jlBegin:jl,:,:), minorfrac(jlBegin:jl,:), & 340 selffac(jlBegin:jl,:),selffrac(jlBegin:jl,:),forfac(jlBegin:jl,:),forfrac(jlBegin:jl,:), & 341 jp(jlBegin:jl,:), rrpk_jt(jlBegin:jl,:,:), (igpt-ngs(ngb(igpt)-1)), indself(jlBegin:jl,:), & 342 indfor(jlBegin:jl,:), indminor(jlBegin:jl,:), & 343 rrpk_taug03(jlBegin:jl,:),fracs(jlBegin:jl,:,ig)) 344 !print *,"Computing" 345 call taumol03_upr(jl,jlBegin,laytrop(1), klev, & 346 rrpk_rat_h2oco2(jlBegin:jl,:,:), colco2(jlBegin:jl,:), colh2o(jlBegin:jl,:), coln2o(jlBegin:jl,:), coldry(jlBegin:jl,:), & 347 rrpk_fac0(jlBegin:jl,:,:), rrpk_fac1(jlBegin:jl,:,:), minorfrac(jlBegin:jl,:), & 348 forfac(jlBegin:jl,:),forfrac(jlBegin:jl,:), & 349 jp(jlBegin:jl,:), rrpk_jt(jlBegin:jl,:,:), (igpt-ngs(ngb(igpt)-1)), & 350 indfor(jlBegin:jl,:), indminor(jlBegin:jl,:), & 351 rrpk_taug03(jlBegin:jl,:),fracs(jlBegin:jl,:,ig)) 352 !print *,"End Computing" 353 END DO 354 CALL system_clock(stop_clock, rate_clock) 355 overall_time=overall_time + (stop_clock-start_clock)/REAL(rate_clock) 356 ENDIF 357 IF(ngb(igpt) == 4) Then 358 !CALL system_clock(start_clock, rate_clock) 359 jl=kproma 360 call taumol04_lwr(jl,laytrop(1), klev, & 361 rrpk_rat_h2oco2(1:jl,:,:), colco2(1:jl,:), colh2o(1:jl,:), coldry(1:jl,:), & 362 rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), & 363 selffac(1:jl,:),selffrac(1:jl,:),forfac(1:jl,:),forfrac(1:jl,:), & 364 jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), indself(1:jl,:), & 365 indfor(1:jl,:), & 366 rrpk_taug04(1:jl,:),fracs(1:jl,:,ig)) 367 call taumol04_upr(jl,laytrop(1), klev, & 368 rrpk_rat_o3co2(1:jl,:,:), colco2(1:jl,:), colo3(1:jl,:), coldry(1:jl,:), & 369 rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), & 370 forfac(1:jl,:),forfrac(1:jl,:), & 371 jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), & 372 indfor(1:jl,:), & 373 rrpk_taug04(1:jl,:),fracs(1:jl,:,ig)) 374 !CALL system_clock(stop_clock, rate_clock) 375 !overall_time=overall_time + (stop_clock-start_clock)/REAL(rate_clock) 376 ENDIF 377 DO jl = 1, kproma 378 ib = ibs(jl, ig) 379 igpt = igs(jl, ig) 380 ! 381 ! Gas concentrations in colxx variables are normalized by 1.e-20_wp in lrtm_coeffs 382 ! CFC gas concentrations (wx) need the same normalization 383 ! Per Eli Mlawer the k values used in gas optics tables have been multiplied by 1e20 384 wx_loc(:,:) = 1.e-20_wp * wx(jl,:,:) 385 IF (ngb(igpt) == 3) THEN 386 taug = rrpk_taug03(jl,:) 387 ELSEIF (ngb(igpt) == 4) THEN 388 taug = rrpk_taug04(jl,:) 389 ELSE 390 CALL gas_optics_lw(klev, igpt, play (jl,:), wx_loc (:,:), coldry (jl,:), laytrop (jl), jp & 391 (jl,:), jt (jl,:), jt1 (jl,:), colh2o (jl,:), colco2 (jl,:), colo3 (jl,:)& 392 , coln2o (jl,:), colco (jl,:), colch4 (jl,:), colo2 (jl,:), colbrd (jl,:), fac00 & 393 (jl,:), fac01 (jl,:), fac10 (jl,:), fac11 (jl,:), rat_h2oco2 (jl,:), rat_h2oco2_1(jl,:), & 394 rat_h2oo3 (jl,:), rat_h2oo3_1 (jl,:), rat_h2on2o (jl,:), rat_h2on2o_1(jl,:), rat_h2och4(jl,:), rat_h2och4_1(& 395 jl,:), rat_n2oco2 (jl,:), rat_n2oco2_1(jl,:), rat_o3co2 (jl,:), rat_o3co2_1 (jl,:), selffac (jl,:), & 396 selffrac (jl,:), indself (jl,:), forfac (jl,:), forfrac (jl,:), indfor (jl,:), minorfrac (& 397 jl,:), scaleminor (jl,:), scaleminorn2(jl,:), indminor (jl,:), fracs (jl,:,ig), taug ) 398 END IF 399 DO jk = 1, klev 400 taut(jl,jk,ig) = taug(jk) + tauaer(jl,jk,ib) 401 END DO 402 END DO ! Loop over columns 403 END DO ! Loop over g point samples - done with gas optical depth calculations 404 PRINT *, "Elapsed time (sec): ", overall_time 405 overall_time=0 406 tautot(1:kproma,:,:) = taut(1:kproma,:,:) + smp_tau(1:kproma,:,:) ! All-sky optical depth. Mask for 0 cloud optical depth? 407 ! 408 ! --- 3.0 Compute radiative transfer. 409 ! -------------------------------- 410 ! 411 ! Initialize fluxes to zero 412 ! 413 uflx(1:kproma,0:klev) = 0.0_wp 414 dflx(1:kproma,0:klev) = 0.0_wp 415 uflxc(1:kproma,0:klev) = 0.0_wp 416 dflxc(1:kproma,0:klev) = 0.0_wp 417 ! 418 ! Planck function in each band at layers and boundaries 419 ! 420 !IBM* ASSERT(NODEPS) 421 DO ig = 1, nbndlw 422 planklay(1:kproma,1:klev,ig) = planckfunction(tlay(1:kproma,1:klev ),ig) 423 planklev(1:kproma,0:klev,ig) = planckfunction(tlev(1:kproma,1:klev+1),ig) 424 plankbnd(1:kproma, ig) = planckfunction(tsfc(1:kproma ),ig) 425 END DO 426 ! 427 ! Precipitable water vapor in each column - this can affect the integration angle secdiff 428 ! 429 pwvcm(1:kproma) = ((amw * sum(wkl(1:kproma,1,1:klev), dim=2)) / (amd * sum(coldry(1:kproma,& 430 1:klev) + wkl(1:kproma,1,1:klev), dim=2))) * (1.e3_wp * psfc(1:kproma)) / (1.e2_wp * grav) 431 ! 432 ! Compute radiative transfer for each set of samples 433 ! 434 DO ig = 1, n_gpts_ts 435 secdiff(1:kproma) = find_secdiff(ibs(1:kproma, ig), pwvcm(1:kproma)) 436 !IBM* ASSERT(NODEPS) 437 DO jl = 1, kproma 438 ib = ibs(jl,ig) 439 layplnk(jl,1:klev) = planklay(jl,1:klev,ib) 440 levplnk(jl,0:klev) = planklev(jl,0:klev,ib) 441 bndplnk(jl) = plankbnd(jl, ib) 442 srfemis(jl) = emis (jl, ib) 443 END DO 444 ! 445 ! All sky fluxes 446 ! 447 CALL lrtm_solver(kproma, kbdim, klev, tautot(:,:,ig), layplnk, levplnk, fracs(:,:,ig), secdiff, bndplnk, srfemis, & 448 zgpfu, zgpfd) 449 uflx(1:kproma,0:klev) = uflx (1:kproma,0:klev) + zgpfu(1:kproma,0:klev) * gpt_scaling 450 dflx(1:kproma,0:klev) = dflx (1:kproma,0:klev) + zgpfd(1:kproma,0:klev) * gpt_scaling 451 ! 452 ! Clear-sky fluxes 453 ! 454 IF (l_do_sep_clear_sky) THEN 455 ! 456 ! Remove clouds and do second RT calculation 457 ! 458 CALL lrtm_solver(kproma, kbdim, klev, taut (:,:,ig), layplnk, levplnk, fracs(:,:,ig), secdiff, bndplnk, & 459 srfemis, zgpcu, zgpcd) 460 uflxc(1:kproma,0:klev) = uflxc(1:kproma,0:klev) + zgpcu(1:kproma,0:klev) * gpt_scaling 461 dflxc(1:kproma,0:klev) = dflxc(1:kproma,0:klev) + zgpcd(1:kproma,0:klev) * gpt_scaling 462 ELSE 463 ! 464 ! Accumulate fluxes by excluding cloudy subcolumns, weighting to account for smaller sample size 465 ! 466 !IBM* ASSERT(NODEPS) 467 DO jk = 0, klev 468 uflxc(1:kproma,jk) = uflxc(1:kproma,jk) & 469 + merge(0._wp, & 470 zgpfu(1:kproma,jk) * clrsky_scaling(1:kproma), & 471 colcldmask(1:kproma,ig)) 472 dflxc(1:kproma,jk) = dflxc(1:kproma,jk) & 473 + merge(0._wp, & 474 zgpfd(1:kproma,jk) * clrsky_scaling(1:kproma), & 475 colcldmask(1:kproma,ig)) 476 END DO 477 END IF 478 END DO ! Loop over samples 479 ! 480 ! --- 3.1 If computing clear-sky fluxes from samples, flag any columns where all samples were cloudy 481 ! 482 ! -------------------------------- 483 IF (.not. l_do_sep_clear_sky) THEN 484 !IBM* ASSERT(NODEPS) 485 DO jl = 1, kproma 486 IF (all(colcldmask(jl,:))) THEN 487 uflxc(jl,0:klev) = rad_undef 488 dflxc(jl,0:klev) = rad_undef 489 END IF 490 END DO 491 END IF 492 END SUBROUTINE lrtm 493 !---------------------------------------------------------------------------- 494 495 elemental FUNCTION planckfunction(temp, band) 496 ! 497 ! Compute the blackbody emission in a given band as a function of temperature 498 ! 499 REAL(KIND=wp), intent(in) :: temp 500 INTEGER, intent(in) :: band 501 REAL(KIND=wp) :: planckfunction 502 INTEGER :: index 503 REAL(KIND=wp) :: fraction 504 index = min(max(1, int(temp - 159._wp)),180) 505 fraction = temp - 159._wp - float(index) 506 planckfunction = totplanck(index, band) + fraction * (totplanck(index+1, band) - totplanck(index, & 507 band)) 508 planckfunction = planckfunction * delwave(band) 509 END FUNCTION planckfunction 510 END MODULE mo_lrtm_driver 511