diff --git a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 index 702da9c..f74bcd4 100644 --- a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 +++ b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 @@ -3694,12 +3694,12 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC class(ty_optical_props_arry), allocatable :: optical_props ! RRTMGP locals - logical :: top_at_1, partial_block, need_aer_optical_props + logical :: top_at_1, need_aer_optical_props logical :: gen_mro, cond_inhomo logical :: rrtmgp_delta_scale, rrtmgp_use_rrtmg_iceflg3_like_forwice integer :: nbnd, ngpt, nmom, icergh integer :: ib, b, nBlocks, colS, colE, ncols_block, & - partial_blockSize, icol, isub, ilay, igpt + icol, isub, ilay, igpt real(wp), allocatable :: t_lev(:) ! (ncol) character(len=ESMF_MAXPATHLEN) :: k_dist_file, cloud_optics_file character(len=ESMF_MAXSTR) :: error_msg @@ -5043,13 +5043,9 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC ! For nstr, must also specify the number of phase function moments (nmom) below. ! ===================================================================================== - ! instantiate optical_props with desired streams - allocate(ty_optical_props_2str::optical_props,__STAT__) ! <-- choose 2-stream SW - - ! initialize spectral discretiz'n and gpt mapping of optical_props - TEST_(optical_props%init(k_dist)) - - ! Used only if nstr (and then must be >= 2) + ! instantiate and initialize optical_props, cloud_props_bnd/gpt, and aer_props + ! inside PROCESS_RRTMGP_BLOCK so each OpenMP thread gets private copies. + ! nmom is only used if nstr (and then must be >= 2) nmom = 2 ! get cloud optical properties (band-only) @@ -5091,31 +5087,6 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC LABEL='RRTMGP_USE_RRTMG_ICEFLG3_LIKE_FORWICE:', & DEFAULT=.TRUE., __RC__) - ! cloud optics file is currently two-stream - ! increment() will handle appropriate stream conversions - allocate(ty_optical_props_2str::cloud_props_bnd_liq,__STAT__) - allocate(ty_optical_props_2str::cloud_props_bnd_ice,__STAT__) - - ! band-only initialization for pre-mcICA cloud optical properties - TEST_(cloud_props_bnd_liq%init(k_dist%get_band_lims_wavenumber())) - TEST_(cloud_props_bnd_ice%init(k_dist%get_band_lims_wavenumber())) - - ! g-point version for McICA sampled cloud optical properties - select type (cloud_props_bnd_liq) - class is (ty_optical_props_2str) - allocate(ty_optical_props_2str::cloud_props_gpt_liq,__STAT__) - class default - TEST_('cloud optical properties (liq) hardwired 2-stream for now') - end select - select type (cloud_props_bnd_ice) - class is (ty_optical_props_2str) - allocate(ty_optical_props_2str::cloud_props_gpt_ice,__STAT__) - class default - TEST_('cloud optical properties (ice) hardwired 2-stream for now') - end select - TEST_(cloud_props_gpt_liq%init(k_dist)) - TEST_(cloud_props_gpt_ice%init(k_dist)) - ! read desired cloud overlap type call MAPL_GetResource( & MAPL, cloud_overlap_type, "RRTMGP_CLOUD_OVERLAP_TYPE_SW:", & @@ -5202,13 +5173,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC ! set up aerosol optical properties need_aer_optical_props = (include_aerosols .and. implements_aerosol_optics) - if (need_aer_optical_props) then - ! aerosol optics system is currently two-stream - ! increment() will handle appropriate stream conversions - allocate(ty_optical_props_2str::aer_props,__STAT__) - ! band-only initialization - TEST_(aer_props%init(k_dist%get_band_lims_wavenumber())) - end if + ! aer_props is allocated and initialized inside PROCESS_RRTMGP_BLOCK !-------------------------------------------------------! ! Loop over blocks of blockSize columns ! @@ -5220,1701 +5185,2360 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC rrtmgp_blockSize, "RRTMGP_SW_BLOCKSIZE:", DEFAULT=4, __RC__) _ASSERT(rrtmgp_blockSize >= 1, 'bad RRTMGP_SW_BLOCKSIZE') - ! for random numbers, for efficiency, reserve the maximum possible - ! subset of columns (rrtmgp_blockSize) since column index is last - allocate(urand(ngpt,LM,rrtmgp_blocksize),__STAT__) - if (gen_mro) then - allocate(urand_aux(ngpt,LM,rrtmgp_blocksize),__STAT__) - if (cond_inhomo) then - allocate(urand_cond (ngpt,LM,rrtmgp_blocksize),__STAT__) - allocate(urand_cond_aux(ngpt,LM,rrtmgp_blocksize),__STAT__) - end if - end if + ! Total number of blocks, including any final partial block. + ! Each block has ncols_block = min(rrtmgp_blockSize, remaining columns), + ! computed at the top of each iteration below. + nBlocks = (ncol + rrtmgp_blockSize - 1) / rrtmgp_blockSize - ! number of FULL blocks by integer division - nBlocks = ncol/rrtmgp_blockSize + ! loop over all blocks + ! RRTMGP is thread-safe (confirmed by RRTMGP developers). + ! All per-block state is private inside PROCESS_RRTMGP_BLOCK. + ! Output arrays are indexed by non-overlapping colS:colE ranges -- no race. + ! NOTE: MAPL timer calls inside the parallel region are guarded by !$OMP CRITICAL. + !$OMP PARALLEL DO SCHEDULE(DYNAMIC) + do b = 1,nBlocks - ! allocate intermediate arrays for FULL blocks - if (nBlocks > 0) then + call PROCESS_RRTMGP_BLOCK( & + b, ncol, rrtmgp_blockSize, LM, ngpt, nbnd, nmom, & + LCLDLM, LCLDMH, include_aerosols, gen_mro, cond_inhomo, & + cloud_overlap_type, IM_World, seeds(2), & + need_aer_optical_props, top_at_1, & + rrtmgp_delta_scale, rrtmgp_use_rrtmg_iceflg3_like_forwice, & + cwp_fac, & + gas_concs, k_dist, cloud_optics, & + p_lay, p_lev, t_lay, QQ3, RR3, dp_wp, dummy_wp, & + CL, dzmid, adl, rdl, Ig1D, Jg1D, band_lims_gpt, & + tsi, mu0, sfc_alb_dir, sfc_alb_dif, taua, ssaa, asya, & + flux_up_clrsky, flux_net_clrsky, & + flux_up_allsky, flux_net_allsky, & + bnd_flux_dn_allsky, bnd_flux_dir_allsky, bnd_flux_net_allsky, & + CLDTS, CLDHS, CLDMS, CLDLS, & + COTTP, COTHP, COTMP, COTLP, & + COTDTP, COTDHP, COTDMP, COTDLP, & + COTNTP, COTNHP, COTNMP, COTNLP, & +#ifdef SOLAR_RADVAL + TAUTP, TAUHP, TAUMP, TAULP, & + COTLDTP, COTLDHP, COTLDMP, COTLDLP, & + COTLNTP, COTLNHP, COTLNMP, COTLNLP, & + COTIDTP, COTIDHP, COTIDMP, COTIDLP, & + COTINTP, COTINHP, COTINMP, COTINLP, & + SSALDTP, SSALDHP, SSALDMP, SSALDLP, & + SSALNTP, SSALNHP, SSALNMP, SSALNLP, & + SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & + SSAINTP, SSAINHP, SSAINMP, SSAINLP, & + ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & + ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & + ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & + ASMINTP, ASMINHP, ASMINMP, ASMINLP, & + CDSDTP, CDSDHP, CDSDMP, CDSDLP, & + CDSNTP, CDSNHP, CDSNMP, CDSNLP, & + CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & + CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & + CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & + CDSINTP, CDSINHP, CDSINMP, CDSINLP, & + SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & + SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & + SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & + SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & + ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & + ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & + ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + FORLDTP, FORLDHP, FORLDMP, FORLDLP, & + FORLNTP, FORLNHP, FORLNMP, FORLNLP, & + FORIDTP, FORIDHP, FORIDMP, FORIDLP, & + FORINTP, FORINHP, FORINMP, FORINLP, & +#endif + MAPL, __RC__) - ! block size UNTIL possible final partial block - ncols_block = rrtmgp_blockSize + end do ! loop over blocks + !$OMP END PARALLEL DO - allocate(toa_flux(ncols_block,ngpt),__STAT__) - allocate(cld_mask(ncols_block,LM,ngpt),__STAT__) - allocate(forwliq(ncols_block,LM,ngpt),__STAT__) - allocate(forwice(ncols_block,LM,ngpt),__STAT__) - if (gen_mro) then - allocate(alpha(ncols_block,LM-1),__STAT__) - if (cond_inhomo) then - allocate(rcorr(ncols_block,LM-1),__STAT__) - allocate(zcw(ncols_block,LM,ngpt),__STAT__) - endif - endif - if (include_aerosols) & - allocate(ClearCounts(4,ncols_block),__STAT__) + call MAPL_TimerOn(MAPL,"--RRTMGP_POST",__RC__) - ! in-cloud cloud optical props - select type (cloud_props_bnd_liq) - class is (ty_optical_props_2str) - TEST_(cloud_props_bnd_liq%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_bnd_ice) - class is (ty_optical_props_2str) - TEST_(cloud_props_bnd_ice%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_gpt_liq) - class is (ty_optical_props_2str) - TEST_(cloud_props_gpt_liq%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_gpt_ice) - class is (ty_optical_props_2str) - TEST_(cloud_props_gpt_ice%alloc_2str(ncols_block,LM)) - end select + ! normalize by incoming solar radiation + allocate(flux_dn_top(ncol), __STAT__) + flux_dn_top(:) = real(max(SLR1D, 1e-7), kind=wp) + do k = 1, LM+1 + flux_up_clrsky (:,k) = flux_up_clrsky (:,k) / flux_dn_top(:) + flux_net_clrsky(:,k) = flux_net_clrsky(:,k) / flux_dn_top(:) + end do + do k = 1, LM+1 + flux_up_allsky (:,k) = flux_up_allsky (:,k) / flux_dn_top(:) + flux_net_allsky(:,k) = flux_net_allsky(:,k) / flux_dn_top(:) + end do + do ib = 1, nbnd + do k = 1, LM+1 + bnd_flux_dn_allsky (:,k,ib) = bnd_flux_dn_allsky (:,k,ib) / flux_dn_top(:) + bnd_flux_net_allsky(:,k,ib) = bnd_flux_net_allsky(:,k,ib) / flux_dn_top(:) + bnd_flux_dir_allsky(:,k,ib) = bnd_flux_dir_allsky(:,k,ib) / flux_dn_top(:) + end do + end do + deallocate(flux_dn_top, __STAT__) - ! aerosol optical props - if (need_aer_optical_props) then - select type (aer_props) - class is (ty_optical_props_2str) - TEST_(aer_props%alloc_2str(ncols_block,LM)) - end select + ! load output arrays + ! clear-sky fluxes + FSCU = real(flux_up_clrsky) + FSC = real(flux_net_clrsky) + ! all-sky fluxes + FSWU = real(flux_up_allsky) + FSW = real(flux_net_allsky) + ! surface net flux per band + do ib = 1, nbnd + FSWBAND(:,ib) = real(bnd_flux_net_allsky(:,LM+1,ib)) + end do + ! surface downwelling direct and diffuse fluxes in bands + if (SOLAR_TO_OBIO .and. include_aerosols) then + do ib = 1, nbnd + DRBAND(:,ib) = real(bnd_flux_dir_allsky(:,LM+1,ib)) + DFBAND(:,ib) = real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) + end do + endif + ! TOA band fluxes + if (include_aerosols) then + if (USE_RRTMG .or. USE_RRTMGP) then + do ib = 1, nbnd + if (band_output(ib)) then + ISRBRGN(ib) % p = real(bnd_flux_dn_allsky(:,1,ib)) + OSRBRGN(ib) % p = real(bnd_flux_dn_allsky(:,1,ib) - bnd_flux_net_allsky(:,1,ib)) + end if + end do end if + end if - ! gas+aer+cld optical properties - select type (optical_props) - class is (ty_optical_props_1scl) - TEST_(optical_props%alloc_1scl(ncols_block,LM)) - class is (ty_optical_props_2str) - TEST_(optical_props%alloc_2str(ncols_block,LM)) - class is (ty_optical_props_nstr) - TEST_(optical_props%alloc_nstr(nmom,ncols_block,LM)) - end select + ! surface direct and diffuse downward in super-bands + ! for *diffuse* downward must subtract direct (downward) from total downward + ! pmn: may later do this using a flux class extension?? + + ! NIR bands (1-9: 820-12850 cm-1, 0.778-12.195 microns) + NIRR = 0.; NIRF = 0. + do ib=1,9 + NIRR = NIRR + real(bnd_flux_dir_allsky(:,LM+1,ib)) + NIRF = NIRF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) + end do + ! PAR bands (11-12: 16000-29000 cm-1, 0.345-0.625 micron) + PARR = 0.; PARF = 0. + do ib=11,12 + PARR = PARR + real(bnd_flux_dir_allsky(:,LM+1,ib)) + PARF = PARF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) + enddo + ! UVR bands (13-14: 29000-50000 cm-1, 0.200-0.345 micron) + UVRR = 0.; UVRF = 0. + do ib=13,14 + UVRR = UVRR + real(bnd_flux_dir_allsky(:,LM+1,ib)) + UVRF = UVRF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) + enddo + ! Transition band (10, 12850-16000 cm-1, 0.625-0.778 micron) + ! split half-and-half to PAR and NIR + NIRR = NIRR + 0.5 * real(bnd_flux_dir_allsky(:,LM+1,10)) + PARR = PARR + 0.5 * real(bnd_flux_dir_allsky(:,LM+1,10)) + NIRF = NIRF + 0.5 * real(bnd_flux_dn_allsky (:,LM+1,10) - bnd_flux_dir_allsky(:,LM+1,10)) + PARF = PARF + 0.5 * real(bnd_flux_dn_allsky (:,LM+1,10) - bnd_flux_dir_allsky(:,LM+1,10)) + ! clean up + deallocate(band_lims_gpt,__STAT__) + deallocate(tsi,mu0,sfc_alb_dir,sfc_alb_dif,__STAT__) + deallocate(dummy_wp,p_lay,t_lay,p_lev,dp_wp,dzmid,__STAT__) + deallocate(flux_up_clrsky,flux_net_clrsky,__STAT__) + deallocate(flux_up_allsky,flux_net_allsky,__STAT__) + deallocate(bnd_flux_dn_allsky,bnd_flux_net_allsky,bnd_flux_dir_allsky,__STAT__) + deallocate(seeds,__STAT__) + if (gen_mro) then + deallocate(adl,__STAT__) + if (cond_inhomo) then + deallocate(rdl,__STAT__) + endif end if + call cloud_optics%finalize() + ! cloud_props_gpt/bnd, aer_props, optical_props are local to PROCESS_RRTMGP_BLOCK + ! and are finalized automatically when that subroutine returns. - ! add final partial block if necessary - partial_block = mod(ncol,rrtmgp_blockSize) /= 0 - if (partial_block) then - partial_blockSize = ncol - nBlocks * rrtmgp_blockSize - nBlocks = nBlocks + 1 - endif + call MAPL_TimerOff(MAPL,"--RRTMGP_POST",__RC__) - ! loop over all blocks - do b = 1,nBlocks + call MAPL_TimerOff(MAPL,"-RRTMGP",__RC__) - ! only the FINAL block can be partial - if (b == nBlocks .and. partial_block) then - ncols_block = partial_blockSize - - if (b > 1) then - ! one or more full blocks already processed - deallocate(toa_flux, __STAT__) - deallocate(cld_mask, __STAT__) - deallocate(forwliq, __STAT__) - deallocate(forwice, __STAT__) - if (gen_mro) then - deallocate(alpha, __STAT__) - if (cond_inhomo) then - deallocate(rcorr,zcw, __STAT__) - endif - endif - if (include_aerosols) & - deallocate(ClearCounts, __STAT__) - endif +#undef TEST_ - allocate(toa_flux(ncols_block,ngpt), __STAT__) - allocate(cld_mask(ncols_block,LM,ngpt), __STAT__) - allocate(forwliq(ncols_block,LM,ngpt), __STAT__) - allocate(forwice(ncols_block,LM,ngpt), __STAT__) - if (gen_mro) then - allocate(alpha(ncols_block,LM-1), __STAT__) - if (cond_inhomo) then - allocate(rcorr(ncols_block,LM-1), __STAT__) - allocate(zcw(ncols_block,LM,ngpt), __STAT__) - endif - endif - if (include_aerosols) & - allocate(ClearCounts(4,ncols_block), __STAT__) + else if (USE_RRTMG) then - ! ty_optical_props routines have an internal deallocation - select type (cloud_props_bnd_liq) - class is (ty_optical_props_2str) - TEST_(cloud_props_bnd_liq%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_bnd_ice) - class is (ty_optical_props_2str) - TEST_(cloud_props_bnd_ice%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_gpt_liq) - class is (ty_optical_props_2str) - TEST_(cloud_props_gpt_liq%alloc_2str(ncols_block,LM)) - end select - select type (cloud_props_gpt_ice) - class is (ty_optical_props_2str) - TEST_(cloud_props_gpt_ice%alloc_2str(ncols_block,LM)) - end select - if (need_aer_optical_props) then - select type (aer_props) - class is (ty_optical_props_2str) - TEST_(aer_props%alloc_2str(ncols_block,LM)) - end select - end if - select type (optical_props) - class is (ty_optical_props_1scl) - TEST_(optical_props%alloc_1scl(ncols_block,LM)) - class is (ty_optical_props_2str) - TEST_(optical_props%alloc_2str(ncols_block,LM)) - class is (ty_optical_props_nstr) - TEST_(optical_props%alloc_nstr(nmom,ncols_block,LM)) - end select + ! regular RRTMG + call MAPL_TimerOn(MAPL,"-RRTMG") - endif ! partial block - - ! prepare block - colS = (b-1) * rrtmgp_blockSize + 1 - colE = colS + ncols_block - 1 - TEST_(gas_concs%get_subset(colS,ncols_block,gas_concs_block)) - - call MAPL_TimerOn(MAPL,"--RRTMGP_GAS_OPTICS",__RC__) - ! gas optics, including source functions - error_msg = k_dist%gas_optics( & - p_lay(colS:colE,:), p_lev(colS:colE,:), t_lay(colS:colE,:), & - gas_concs_block, optical_props, toa_flux) - TEST_(error_msg) - call MAPL_TimerOff(MAPL,"--RRTMGP_GAS_OPTICS",__RC__) - - ! get block of aerosol optical props - if (need_aer_optical_props) then - select type (aer_props) - class is (ty_optical_props_2str) - - ! load un-normalized optical properties from aerosol system - aer_props%tau = real(TAUA(colS:colE,:,:),kind=wp) - aer_props%ssa = real(SSAA(colS:colE,:,:),kind=wp) - aer_props%g = real(ASYA(colS:colE,:,:),kind=wp) - - ! renormalize - where (aer_props%tau > 0._wp .and. aer_props%ssa > 0._wp) - aer_props%g = aer_props%g / aer_props%ssa - aer_props%ssa = aer_props%ssa / aer_props%tau - elsewhere - aer_props%tau = 0._wp - aer_props%ssa = 0._wp - aer_props%g = 0._wp - end where + ! reversed (flipped) vertical dimension arrays and other RRTMG arrays + ! ------------------------------------------------------------------- - ! Because RRTMGP is (currently) compiled at R8, - ! _wp is R8. Apparently with aggressive compiler - ! flags using Intel, it's possible for, say, - ! aer_props%ssa to become slightly greater than one - ! in the above renormalization. So, we add clamps - ! to the values based on the restrictions see in - ! RRTMGP/rte-frontend/mo_optical_props.F90 - ! - ! In testing, the values seen were like 1.00000011905028 - ! so just slightly above one. - - ! tau must be greater than 0.0 - aer_props%tau = max(aer_props%tau, 0._wp) - ! ssa must be between 0.0 and 1.0 - aer_props%ssa = max(min(aer_props%ssa, 1._wp), 0._wp) - ! g must be between -1.0 and 1.0 - aer_props%g = max(min(aer_props%g, 1._wp),-1._wp) - - class default - TEST_('aerosol optical properties hardwired 2-stream for now') - end select - end if + ! interface (between layer) variables + allocate(TLEV (NCOL,LM+1),__STAT__) + allocate(TLEV_R(NCOL,LM+1),__STAT__) + allocate(PLE_R (NCOL,LM+1),__STAT__) + ! cloud physical properties + allocate(FCLD_R(NCOL,LM),__STAT__) + allocate(CLIQWP(NCOL,LM),__STAT__) + allocate(CICEWP(NCOL,LM),__STAT__) + allocate(RELIQ (NCOL,LM),__STAT__) + allocate(REICE (NCOL,LM),__STAT__) + ! aerosol optical properties + allocate(TAUAER(NCOL,LM,NB_RRTMG),__STAT__) + allocate(SSAAER(NCOL,LM,NB_RRTMG),__STAT__) + allocate(ASMAER(NCOL,LM,NB_RRTMG),__STAT__) + ! layer variables + allocate(DPR (NCOL,LM),__STAT__) + allocate(PL_R (NCOL,LM),__STAT__) + allocate(ZL_R (NCOL,LM),__STAT__) + allocate(T_R (NCOL,LM),__STAT__) + allocate(Q_R (NCOL,LM),__STAT__) + allocate(O2_R (NCOL,LM),__STAT__) + allocate(O3_R (NCOL,LM),__STAT__) + allocate(CO2_R (NCOL,LM),__STAT__) + allocate(CH4_R (NCOL,LM),__STAT__) + ! super-layer cloud fractions + allocate(CLEARCOUNTS (NCOL,4),__STAT__) + ! output fluxes + allocate(SWUFLX (NCOL,LM+1),__STAT__) + allocate(SWDFLX (NCOL,LM+1),__STAT__) + allocate(SWUFLXC(NCOL,LM+1),__STAT__) + allocate(SWDFLXC(NCOL,LM+1),__STAT__) + ! un-flipped outputs + allocate(SWUFLXR (NCOL,LM+1),__STAT__) + allocate(SWDFLXR (NCOL,LM+1),__STAT__) + allocate(SWUFLXCR(NCOL,LM+1),__STAT__) + allocate(SWDFLXCR(NCOL,LM+1),__STAT__) - call MAPL_TimerOn(MAPL,"--RRTMGP_CLOUD_OPTICS",__RC__) - - ! Make band in-cloud optical props from cloud_optics and mean in-cloud cloud water paths. - ! These can be scaled later to account for sub-gridscale condensate inhomogeneity. - ! Do phases separately to allow for different forward scattering, etc., per earlier note. - ! liquid ... - error_msg = cloud_optics%cloud_optics( & - real(QQ3(colS:colE,:,2),kind=wp) * dp_wp(colS:colE,:) * cwp_fac, & ! [g/m2] - dummy_wp(colS:colE,:), & - min( max( real(RR3(colS:colE,:,2),kind=wp), & ! [microns] - cloud_optics%get_min_radius_liq()), & - cloud_optics%get_max_radius_liq()), & - dummy_wp(colS:colE,:), & - cloud_props_bnd_liq) - TEST_(error_msg) - ! ice ... - error_msg = cloud_optics%cloud_optics( & - dummy_wp(colS:colE,:), & - real(QQ3(colS:colE,:,1),kind=wp) * dp_wp(colS:colE,:) * cwp_fac, & ! [g/m2] - dummy_wp(colS:colE,:), & - min( max( real(RR3(colS:colE,:,1),kind=wp), & ! [microns] - cloud_optics%get_min_radius_ice()), & - cloud_optics%get_max_radius_ice()), & - cloud_props_bnd_ice) - TEST_(error_msg) - - call MAPL_TimerOff(MAPL,"--RRTMGP_CLOUD_OPTICS",__RC__) - - call MAPL_TimerOn(MAPL,"--RRTMGP_MCICA",__RC__) + ! Set flags related to cloud properties (see RRTMG_SW) + ! ---------------------------------------------------- + call MAPL_GetResource(MAPL,ICEFLGSW,'RRTMG_ICEFLG:',DEFAULT=3,__RC__) + call MAPL_GetResource(MAPL,LIQFLGSW,'RRTMG_LIQFLG:',DEFAULT=1,__RC__) -!!TODO: need to resolve diff between prob of max vs ran and correlation coeff in both paper and code + if (LM > 72) then + call MAPL_GetResource(MAPL,USE_PRECIP_IN_RADIATION,'RRTMGSW_USE_PRECIP_IN_RADIATION:',DEFAULT=.TRUE.,RC=STATUS) + VERIFY_(STATUS) + else + call MAPL_GetResource(MAPL,USE_PRECIP_IN_RADIATION,'RRTMGSW_USE_PRECIP_IN_RADIATION:',DEFAULT=.FALSE.,RC=STATUS) + VERIFY_(STATUS) + endif - ! exponential inter-layer correlations - ! [alpha|rcorr](k) is correlation between layers k and k+1 - ! dzmid(k) is separation between midpoints of layers k and k+1 - if (gen_mro) then - do ilay = 1,LM-1 - ! cloud fraction correlation - alpha(:,ilay) = exp(-abs(dzmid(colS:colE,ilay))/real(adl(colS:colE),kind=wp)) - enddo - if (cond_inhomo) then - do ilay = 1,LM-1 - ! condensate correlation - rcorr(:,ilay) = exp(-abs(dzmid(colS:colE,ilay))/real(rdl(colS:colE),kind=wp)) - enddo - endif - endif + ! Normalize aerosol inputs + ! ------------------------ - ! generate McICA random numbers for block - ! Perhaps later this can be parallelized? - do isub = 1, ncols_block - ! local 1d column index - icol = colS + isub - 1 - ! initialize the Philox PRNG - ! set word1 of key based on GLOBAL location - ! 32-bits can hold all forseeable resolutions - seeds(1) = nint(Jg1D(icol)) * IM_World + nint(Ig1D(icol)) -#ifdef HAVE_MKL - ! instantiate a random number stream for the column - call rng%init(VSL_BRNG_PHILOX4X32X10,seeds) -#else - call rng%init(seeds) -#endif - ! draw the random numbers for the column - urand(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) - if (gen_mro) then - urand_aux(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) - if (cond_inhomo) then - urand_cond (:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) - urand_cond_aux(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) - endif - end if - ! free the rng - call rng%end() - end do + if (num_aero_vars > 0) then + where (TAUA > 0. .and. SSAA > 0. ) + ASYA = ASYA/SSAA + SSAA = SSAA/TAUA + elsewhere + TAUA = 0. + SSAA = 0. + ASYA = 0. + end where + end if - ! cloud sampling to gpoints - select case (cloud_overlap_type) - case ("MAX_RAN_OVERLAP") - error_msg = sampled_mask_max_ran( & - urand(:,:,1:ncols_block), real(CL(colS:colE,:),kind=wp), cld_mask) - TEST_(error_msg) - case ("EXP_RAN_OVERLAP") - ! corr_coeff(ncols_block,LM-1) is an inter-layer correlation coefficient - ! to be provided ... it is not the same as alpha, which is a probability -! error_msg = sampled_mask_exp_ran( & -! urand(:,:,1:ncols_block), real(CL(colS:colE,:),kind=wp), corr_coeff, cld_mask) -! TEST_(error_msg) - TEST_('EXP_RAN_OVERLAP not implemented yet') - case ("GEN_MAX_RAN_OVERLAP") - ! a scheme like Oreopoulos et al. 2012 (doi:10.5194/acp-12-9097-2012) in which both - ! cloud presence and cloud condensate are separately generalized maximum-random: - error_msg = sampled_urand_gen_max_ran(alpha, & - urand(:,:,1:ncols_block),urand_aux(:,:,1:ncols_block)) - TEST_(error_msg) - if (cond_inhomo) then - error_msg = sampled_urand_gen_max_ran(rcorr, & - urand_cond(:,:,1:ncols_block),urand_cond_aux(:,:,1:ncols_block)) - TEST_(error_msg) - end if - do isub = 1,ncols_block - icol = colS + isub - 1 - do ilay = 1,LM - cld_frac = CL(icol,ilay) + ! Flip in vertical, Convert units, and interpolate T, etc. + ! -------------------------------------------------------- + ! RRTMG convention is that vertical indices increase from bot -> top - ! if grid-box clear, no subgrid variability - if (cld_frac <= 0.) then - cld_mask(isub,ilay,:) = .false. - else - ! subgrid-scale cloud mask - cld_mask(isub,ilay,:) = urand(:,ilay,isub) < cld_frac - - ! subgrid-scale condensate - if (cond_inhomo) then - ! level of condensate inhomogeneity based on cloud fraction. - if (cld_frac > 0.99) then - sigma_qcw = 0.5 - elseif (cld_frac > 0.9) then - sigma_qcw = 0.71 - else - sigma_qcw = 1.0 - endif - do igpt = 1,ngpt - if (cld_mask(isub,ilay,igpt)) zcw(isub,ilay,igpt) = & - zcw_lookup(real(urand_cond(igpt,ilay,isub)),sigma_qcw) - end do - end if - end if + call MAPL_TimerOn (MAPL,"--RRTMG_FLIP") - end do - end do + DPR(:,1:LM) = (PLE(:,2:LM+1)-PLE(:,1:LM)) - case default - TEST_('RRTMGP_SW: unknown cloud overlap') - end select + ! cloud water paths converted from g/g to g/m^2 + ! Flip in vertical + do K = 1,LM + LV = LM-K+1 ! LM --> 1 - ! draw McICA optical property samples (band->gpt) - TEST_(draw_samples(cld_mask, cloud_props_bnd_liq, cloud_props_gpt_liq)) - TEST_(draw_samples(cld_mask, cloud_props_bnd_ice, cloud_props_gpt_ice)) + ! Convert content [kg/kg] to path [g/m2] + ! using hydrostatic eqn dp/g ~ rho*dz, + ! so conversion factor is 1000*dp/g ~ 1.02*100*dp. + ! pmn: why not use MAPL_GRAV explicitly? + if (USE_PRECIP_IN_RADIATION) then + ILWT = QQ3(:,LV,KLIQUID)+QQ3(:,LV,KRAIN) + CLIQWP(:,K) = 1.02*100*DPR(:,LV)*ILWT + where (ILWT > 0.0) + RELIQ (:,K) = ( RR3(:,LV,KLIQUID)*QQ3(:,LV,KLIQUID) + & + RR3(:,LV,KRAIN )*QQ3(:,LV,KRAIN ) ) / ILWT + elsewhere + RELIQ (:,K) = 14.0 + end where + ILWT = QQ3(:,LV,KICE)+QQ3(:,LV,KSNOW)+QQ3(:,LV,KGRAUPEL) + CICEWP(:,K) = 1.02*100*DPR(:,LV)*ILWT + WHERE (ILWT > 0.0) + REICE (:,K) = ( RR3(:,LV,KICE )*QQ3(:,LV,KICE ) + & + RR3(:,LV,KSNOW )*QQ3(:,LV,KSNOW ) + & + RR3(:,LV,KGRAUPEL)*QQ3(:,LV,KGRAUPEL) ) / ILWT + elsewhere + REICE (:,K) = 36.0 + end where + else + CLIQWP(:,K) = 1.02*100*DPR(:,LV)*QQ3(:,LV,KLIQUID) + CICEWP(:,K) = 1.02*100*DPR(:,LV)*QQ3(:,LV,KICE) + RELIQ (:,K) = RR3(:,LV,KLIQUID) + REICE (:,K) = RR3(:,LV,KICE ) + endif + enddo - ! Scaling to sub-gridscale water paths: - ! since tau for each phase is linear in the phase's water path - ! and since the scaling zcw applies equally to both phases, the - ! total g-point optical thickness tau will scale with zcw. - if (gen_mro) then - if (cond_inhomo) & - where (cld_mask) cloud_props_gpt_liq%tau = cloud_props_gpt_liq%tau * zcw - where (cld_mask) cloud_props_gpt_ice%tau = cloud_props_gpt_ice%tau * zcw - end if + IF (ICEFLGSW == 0) THEN + WHERE (REICE < 10.) REICE = 10. + WHERE (REICE > 30.) REICE = 30. + ELSE IF (ICEFLGSW == 1) THEN + WHERE (REICE < 13.) REICE = 13. + WHERE (REICE > 130.) REICE = 130. + ELSE IF (ICEFLGSW == 2) THEN + WHERE (REICE < 5.) REICE = 5. + WHERE (REICE > 131.) REICE = 131. + ELSE IF (ICEFLGSW == 3) THEN + WHERE (REICE < 5.) REICE = 5. + WHERE (REICE > 140.) REICE = 140. + ELSE IF (ICEFLGSW == 4) THEN + REICE(:,:) = REICE(:,:)*2. + WHERE (REICE < 1.) REICE = 1. + WHERE (REICE > 200.) REICE = 200. + END IF - call MAPL_TimerOff(MAPL,"--RRTMGP_MCICA",__RC__) + IF (LIQFLGSW == 0) THEN + WHERE (RELIQ < 10.) RELIQ = 10. + WHERE (RELIQ > 30.) RELIQ = 30. + ELSE IF (LIQFLGSW == 1) THEN + WHERE (RELIQ < 2.5) RELIQ = 2.5 + WHERE (RELIQ > 60.) RELIQ = 60. + END IF - ! REFRESH super-layer diagnostics (before delta-scaling TAUs). - ! ** Calculated from subcolumn ensemble, so stochastic ** - ! ------------------------------------------------------- - call MAPL_TimerOn(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - if (include_aerosols) then + ! regular (non-flipped) interface temperatures + TLEV(:,2:LM)=(T(:,1:LM-1)* DPR(:,2:LM) + T(:,2:LM) * DPR(:,1:LM-1)) & + /(DPR(:,1:LM-1) + DPR(:,2:LM)) + TLEV(:,LM+1) = TS(:) + TLEV(:, 1) = TLEV(:,2) - ! super-layer cloud fractions - call clearCounts_threeBand( & - ncols_block, ncols_block, ngpt, LM, LCLDLM, LCLDMH, & - reshape(cld_mask,[LM,ngpt,ncols_block],order=[3,1,2]), & - ClearCounts) - do isub = 1,ncols_block - icol = colS + isub - 1 - CLDTS(icol) = 1. - ClearCounts(1,isub)/float(ngpt) - CLDHS(icol) = 1. - ClearCounts(2,isub)/float(ngpt) - CLDMS(icol) = 1. - ClearCounts(3,isub)/float(ngpt) - CLDLS(icol) = 1. - ClearCounts(4,isub)/float(ngpt) - end do + ! flip in vertical ... - ! in-cloud optical thicknesses in PAR super-band - ! (weighted across and within bands by TOA incident flux) - do isub = 1,ncols_block - icol = colS + isub - 1 + PLE_R (:,1:LM+1) = PLE (:,LM+1:1:-1) / 100. ! hPa + TLEV_R(:,1:LM+1) = TLEV(:,LM+1:1:-1) -#ifdef SOLAR_RADVAL - ! default (no cloud) for TAUx variant - TAUTP(icol) = 0. - TAUHP(icol) = 0. - TAUMP(icol) = 0. - TAULP(icol) = 0. -#endif + PL_R (:,1:LM ) = PL (:,LM:1:-1) / 100. ! hPa + T_R (:,1:LM ) = T (:,LM:1:-1) - ! default (no cloud) for COTx variant - COTTP(icol) = MAPL_UNDEF - COTHP(icol) = MAPL_UNDEF - COTMP(icol) = MAPL_UNDEF - COTLP(icol) = MAPL_UNDEF + ! Specific humidity is converted to Volume Mixing Ratio + Q_R (:,1:LM ) = Q (:,LM:1:-1) / (1.-Q(:,LM:1:-1)) * (MAPL_AIRMW/MAPL_H2OMW) - ! zero denom- and numerator accumulators - COTDTP(icol) = 0.; COTNTP(icol) = 0. - COTDHP(icol) = 0.; COTNHP(icol) = 0. - COTDMP(icol) = 0.; COTNMP(icol) = 0. - COTDLP(icol) = 0.; COTNLP(icol) = 0. -#ifdef SOLAR_RADVAL - COTLDTP(icol) = 0.; COTLNTP(icol) = 0.; COTIDTP(icol) = 0.; COTINTP(icol) = 0. - COTLDHP(icol) = 0.; COTLNHP(icol) = 0.; COTIDHP(icol) = 0.; COTINHP(icol) = 0. - COTLDMP(icol) = 0.; COTLNMP(icol) = 0.; COTIDMP(icol) = 0.; COTINMP(icol) = 0. - COTLDLP(icol) = 0.; COTLNLP(icol) = 0.; COTIDLP(icol) = 0.; COTINLP(icol) = 0. - SSALDTP(icol) = 0.; SSALNTP(icol) = 0.; SSAIDTP(icol) = 0.; SSAINTP(icol) = 0. - SSALDHP(icol) = 0.; SSALNHP(icol) = 0.; SSAIDHP(icol) = 0.; SSAINHP(icol) = 0. - SSALDMP(icol) = 0.; SSALNMP(icol) = 0.; SSAIDMP(icol) = 0.; SSAINMP(icol) = 0. - SSALDLP(icol) = 0.; SSALNLP(icol) = 0.; SSAIDLP(icol) = 0.; SSAINLP(icol) = 0. - ASMLDTP(icol) = 0.; ASMLNTP(icol) = 0.; ASMIDTP(icol) = 0.; ASMINTP(icol) = 0. - ASMLDHP(icol) = 0.; ASMLNHP(icol) = 0.; ASMIDHP(icol) = 0.; ASMINHP(icol) = 0. - ASMLDMP(icol) = 0.; ASMLNMP(icol) = 0.; ASMIDMP(icol) = 0.; ASMINMP(icol) = 0. - ASMLDLP(icol) = 0.; ASMLNLP(icol) = 0.; ASMIDLP(icol) = 0.; ASMINLP(icol) = 0. -#endif + ! Ozone is converted Mass Mixing Ratio to Volume Mixing Ratio + O3_R (:,1:LM ) = O3 (:,LM:1:-1) * (MAPL_AIRMW/MAPL_O3MW) - ! can only be non-zero for potentially cloudy columns - if (any(CL(icol,:) > 0.)) then + ! chemistry and cloud fraction + ! (cloud water paths and effective radii flipped already) + CH4_R (:,1:LM ) = CH4(:,LM:1:-1) + CO2_R (:,1:LM ) = CO2 + O2_R (:,1:LM ) = O2 + FCLD_R(:,1:LM ) = CL (:,LM:1:-1) - ! accumulate over gpts/subcolumns - do ib = 1, nbnd - do igpt = band_lims_gpt(1,ib), band_lims_gpt(2,ib) - - ! band weights for photosynthetically active radiation (PAR) - ! Bands 11-12 (0.345-0.625 um) plus half transition band 10 (0.625-0.778 um) - if (ib >= 11 .and. ib <= 12) then - wgt = 1.0 - else if (ib == 10) then - wgt = 0.5 - else - ! no contribution to PAR - cycle - end if +! Clean up negatives + WHERE (Q_R < 0.) Q_R = 0. + WHERE (O3_R < 0.) O3_R = 0. + WHERE (CH4_R < 0.) CH4_R = 0. + WHERE (CO2_R < 0.) CO2_R = 0. + WHERE (O2_R < 0.) O2_R = 0. + WHERE (FCLD_R < 0.) FCLD_R = 0. - ! TOA flux weighting - ! (note: neither the adjustment of toa_flux to our tsi - ! or for zenith angle are needed yet since this weighting - ! is over gpoint and is normalized for EACH icol) - wgt = wgt * toa_flux(isub,igpt) + ! Adjustment for Earth/Sun distance, from MAPL_SunGetInsolation + ADJES = DIST - ! low pressure layer - sltaulp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt)) - sitaulp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt)) - staulp = sltaulp + sitaulp - if (staulp > 0.) then - COTDLP(icol) = COTDLP(icol) + wgt - COTNLP(icol) = COTNLP(icol) + wgt * staulp - end if -#ifdef SOLAR_RADVAL - sltaussalp = 0.; sltaussaglp = 0. - if (sltaulp > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussalp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt)) - sltaussaglp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%g (isub,LCLDLM:LM,igpt)) - end select - COTLDLP(icol) = COTLDLP(icol) + wgt - COTLNLP(icol) = COTLNLP(icol) + wgt * sltaulp - SSALDLP(icol) = SSALDLP(icol) + wgt * sltaulp - SSALNLP(icol) = SSALNLP(icol) + wgt * sltaussalp - ASMLDLP(icol) = ASMLDLP(icol) + wgt * sltaussalp - ASMLNLP(icol) = ASMLNLP(icol) + wgt * sltaussaglp - end if - sitaussalp = 0.; sitaussaglp = 0. - if (sitaulp > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussalp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt)) - sitaussaglp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%g (isub,LCLDLM:LM,igpt)) - end select - COTIDLP(icol) = COTIDLP(icol) + wgt - COTINLP(icol) = COTINLP(icol) + wgt * sitaulp - SSAIDLP(icol) = SSAIDLP(icol) + wgt * sitaulp - SSAINLP(icol) = SSAINLP(icol) + wgt * sitaussalp - ASMIDLP(icol) = ASMIDLP(icol) + wgt * sitaussalp - ASMINLP(icol) = ASMINLP(icol) + wgt * sitaussaglp - end if -#endif + ! Layer mid-point heights relative to zero at index 1 + ZL_R(:,1) = 0. + do k=2,LM + ! dz = RT/g x dp/p + ! Note: This is correct even though its different from LW. + ! Its because SW uses LE[V]_R 1:LM+1 while LW uses 0:LM. + ZL_R(:,k) = ZL_R(:,k-1) + MAPL_RGAS*TLEV_R(:,k)/MAPL_GRAV*(PL_R(:,k-1)-PL_R(:,k))/PLE_R(:,k) + enddo + + ! aerosols + TAUAER(:,1:LM,:) = TAUA(:,LM:1:-1,:) + SSAAER(:,1:LM,:) = SSAA(:,LM:1:-1,:) + ASMAER(:,1:LM,:) = ASYA(:,LM:1:-1,:) + + call MAPL_TimerOff(MAPL,"--RRTMG_FLIP") + call MAPL_TimerOn (MAPL,"--RRTMG_INIT") + + ! initialize RRTMG SW + call RRTMG_SW_INI + + call MAPL_TimerOff(MAPL,"--RRTMG_INIT") + call MAPL_TimerOn (MAPL,"--RRTMG_RUN") + + ! partition size for columns (profiles) used to improve efficiency + call MAPL_GetResource( MAPL, RPART, 'RRTMGSW_PARTITION_SIZE:', DEFAULT=0, __RC__) + + ! various RRTMG configuration options ... + + IAER = 10 ! Per AER: + ! 0: Turns off aerosols + ! 10: Enables aerosols + + NORMFLX = 1 ! 0: Do not normalize fluxes + ! 1: Normalize fluxes + + DYOFYR = DOY ! Day of year + + call MAPL_GetResource( MAPL, ISOLVAR ,'ISOLVAR:', DEFAULT=0, __RC__) + + ! ISOLVAR: + ! Flag for solar variability method + ! -1 = (when scon .eq. 0.0): No solar variability + ! and no solar cycle (Kurucz solar irradiance + ! of 1368.22 Wm-2 only); + ! (when scon .ne. 0.0): Kurucz solar irradiance + ! scaled to scon and solar variability defined + ! (optional) by setting non-zero scale factors + ! for each band in bndsolvar + ! 0 = (when SCON .eq. 0.0): No solar variability + ! and no solar cycle (NRLSSI2 solar constant of + ! 1360.85 Wm-2 for the 100-50000 cm-1 spectral + ! range only), with facular and sunspot effects + ! fixed to the mean of Solar Cycles 13-24; + ! (when SCON .ne. 0.0): No solar variability + ! and no solar cycle (NRLSSI2 solar constant of + ! 1360.85 Wm-2 for the 100-50000 cm-1 spectral + ! range only), is scaled to SCON + ! 1 = Solar variability (using NRLSSI2 solar + ! model) with solar cycle contribution + ! determined by fraction of solar cycle + ! with facular and sunspot variations + ! fixed to their mean variations over the + ! average of Solar Cycles 13-24; + ! two amplitude scale factors allow + ! facular and sunspot adjustments from + ! mean solar cycle as defined by indsolvar + ! 2 = Solar variability (using NRLSSI2 solar + ! model) over solar cycle determined by + ! direct specification of Mg (facular) + ! and SB (sunspot) indices provided + ! in indsolvar (scon = 0.0 only) + ! 3 = (when scon .eq. 0.0): No solar variability + ! and no solar cycle (NRLSSI2 solar irradiance + ! of 1360.85 Wm-2 only); + ! (when scon .ne. 0.0): NRLSSI2 solar irradiance + ! scaled to scon and solar variability defined + ! (optional) by setting non-zero scale factors + ! for each band in bndsolvar + + if (ISOLVAR == 1) then + if (MAPL_AM_I_ROOT()) then + write (*,*) "ERROR in RRTMG_SW" + write (*,*) "ISOLVAR==1 is currently unsupported as we have no" + write (*,*) "way of correctly setting solcycfrac." + end if + _FAIL('RRTMG SW: ISOLVAR==1 currently unsupported') + end if - ! mid pressure layer - sltaump = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt)) - sitaump = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt)) - staump = sltaump + sitaump - if (staump > 0.) then - COTDMP(icol) = COTDMP(icol) + wgt - COTNMP(icol) = COTNMP(icol) + wgt * staump - end if -#ifdef SOLAR_RADVAL - sltaussamp = 0.; sltaussagmp = 0. - if (sltaump > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussamp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt)) - sltaussagmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%g (isub,LCLDMH:LCLDLM-1,igpt)) - end select - COTLDMP(icol) = COTLDMP(icol) + wgt - COTLNMP(icol) = COTLNMP(icol) + wgt * sltaump - SSALDMP(icol) = SSALDMP(icol) + wgt * sltaump - SSALNMP(icol) = SSALNMP(icol) + wgt * sltaussamp - ASMLDMP(icol) = ASMLDMP(icol) + wgt * sltaussamp - ASMLNMP(icol) = ASMLNMP(icol) + wgt * sltaussagmp - end if - sitaussamp = 0.; sitaussagmp = 0. - if (sitaump > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussamp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt)) - sitaussagmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%g (isub,LCLDMH:LCLDLM-1,igpt)) - end select - COTIDMP(icol) = COTIDMP(icol) + wgt - COTINMP(icol) = COTINMP(icol) + wgt * sitaump - SSAIDMP(icol) = SSAIDMP(icol) + wgt * sitaump - SSAINMP(icol) = SSAINMP(icol) + wgt * sitaussamp - ASMIDMP(icol) = ASMIDMP(icol) + wgt * sitaussamp - ASMINMP(icol) = ASMINMP(icol) + wgt * sitaussagmp - end if -#endif + ! INDSOLVAR = Facular and sunspot amplitude + ! scale factors (isolvar=1), or + ! Mg and SB indices (isolvar=2) + ! Dimensions: (2) - ! high pressure layer - sltauhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt)) - sitauhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt)) - stauhp = sltauhp + sitauhp - if (stauhp > 0.) then - COTDHP(icol) = COTDHP(icol) + wgt - COTNHP(icol) = COTNHP(icol) + wgt * stauhp - end if -#ifdef SOLAR_RADVAL - sltaussahp = 0.; sltaussaghp = 0. - if (sltauhp > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussahp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt)) - sltaussaghp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%g (isub,1:LCLDMH-1,igpt)) - end select - COTLDHP(icol) = COTLDHP(icol) + wgt - COTLNHP(icol) = COTLNHP(icol) + wgt * sltauhp - SSALDHP(icol) = SSALDHP(icol) + wgt * sltauhp - SSALNHP(icol) = SSALNHP(icol) + wgt * sltaussahp - ASMLDHP(icol) = ASMLDHP(icol) + wgt * sltaussahp - ASMLNHP(icol) = ASMLNHP(icol) + wgt * sltaussaghp - end if - sitaussahp = 0.; sitaussaghp = 0. - if (sitauhp > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussahp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt)) - sitaussaghp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%g (isub,1:LCLDMH-1,igpt)) - end select - COTIDHP(icol) = COTIDHP(icol) + wgt - COTINHP(icol) = COTINHP(icol) + wgt * sitauhp - SSAIDHP(icol) = SSAIDHP(icol) + wgt * sitauhp - SSAINHP(icol) = SSAINHP(icol) + wgt * sitaussahp - ASMIDHP(icol) = ASMIDHP(icol) + wgt * sitaussahp - ASMINHP(icol) = ASMINHP(icol) + wgt * sitaussaghp - end if -#endif + if (ISOLVAR == 2) then - ! whole subcolumn - sltautp = sltaulp + sltaump + sltauhp - sitautp = sitaulp + sitaump + sitauhp - stautp = staulp + staump + stauhp - if (stautp > 0.) then - COTDTP(icol) = COTDTP(icol) + wgt - COTNTP(icol) = COTNTP(icol) + wgt * stautp - end if -#ifdef SOLAR_RADVAL - sltaussatp = sltaussalp + sltaussamp + sltaussahp - sltaussagtp = sltaussaglp + sltaussagmp + sltaussaghp - if (sltautp > 0.) then - COTLDTP(icol) = COTLDTP(icol) + wgt - COTLNTP(icol) = COTLNTP(icol) + wgt * sltautp - SSALDTP(icol) = SSALDTP(icol) + wgt * sltautp - SSALNTP(icol) = SSALNTP(icol) + wgt * sltaussatp - ASMLDTP(icol) = ASMLDTP(icol) + wgt * sltaussatp - ASMLNTP(icol) = ASMLNTP(icol) + wgt * sltaussagtp - end if - sitaussatp = sitaussalp + sitaussamp + sitaussahp - sitaussagtp = sitaussaglp + sitaussagmp + sitaussaghp - if (sitautp > 0.) then - COTIDTP(icol) = COTIDTP(icol) + wgt - COTINTP(icol) = COTINTP(icol) + wgt * sitautp - SSAIDTP(icol) = SSAIDTP(icol) + wgt * sitautp - SSAINTP(icol) = SSAINTP(icol) + wgt * sitaussatp - ASMIDTP(icol) = ASMIDTP(icol) + wgt * sitaussatp - ASMINTP(icol) = ASMINTP(icol) + wgt * sitaussagtp - end if -#endif + ! Solar indices from our file + INDSOLVAR(1) = MG + INDSOLVAR(2) = SB - end do ! igpt - end do ! ib + else - ! normalize - ! Note: TAUx defaults zero, COTx defaults MAPL_UNDEF - if (COTDTP(icol) > 0. .and. COTNTP(icol) > 0.) then - COTTP(icol) = COTNTP(icol) / COTDTP(icol) -#ifdef SOLAR_RADVAL - TAUTP(icol) = COTTP(icol) -#endif - end if + call MAPL_GetResource( MAPL, INDSOLVAR(1) ,'INDSOLVAR_1:', DEFAULT=1.0, __RC__) + call MAPL_GetResource( MAPL, INDSOLVAR(2) ,'INDSOLVAR_2:', DEFAULT=1.0, __RC__) - if (COTDHP(icol) > 0. .and. COTNHP(icol) > 0.) then - COTHP(icol) = COTNHP(icol) / COTDHP(icol) -#ifdef SOLAR_RADVAL - TAUHP(icol) = COTHP(icol) -#endif - end if + end if - if (COTDMP(icol) > 0. .and. COTNMP(icol) > 0.) then - COTMP(icol) = COTNMP(icol) / COTDMP(icol) -#ifdef SOLAR_RADVAL - TAUMP(icol) = COTMP(icol) -#endif - end if - if (COTDLP(icol) > 0. .and. COTNLP(icol) > 0.) then - COTLP(icol) = COTNLP(icol) / COTDLP(icol) -#ifdef SOLAR_RADVAL - TAULP(icol) = COTLP(icol) -#endif - end if + BNDSOLVAR = 1.0 ! Solar variability scale factors + ! for each shortwave band + ! Dimensions: (nbndsw=14) - end if ! potentially cloudy column - end do ! isub - end if ! include_aerosols - call MAPL_TimerOff(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) + ! SOLCYCFRAC: Fraction of averaged 11-year solar cycle (0-1) + ! at current time (isolvar=1) + ! 0.0 represents the first day of year 1 + ! 1.0 represents the last day of year 11 - ! delta-scaling of cloud optical properties (accounts for forward scattering) - call MAPL_TimerOn(MAPL,"--RRTMGP_DELTA_SCALE",__RC__) - forwliq = 0.; forwice = 0. ! default for no delta-scaling - if (rrtmgp_delta_scale) then + ! MAT: Note while we don't currently use SOLCYCFRAC, we set it to something + ! to avoid an optional variable on GPUs - ! default delta-scaling for liquid - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - forwliq = cloud_props_gpt_liq%g ** 2 - end select - TEST_(cloud_props_gpt_liq%delta_scale(forwliq)) - - if (rrtmgp_use_rrtmg_iceflg3_like_forwice) then - ! non-default delta-scaling for ice (as in RRTMG iceflag==3) - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - radice_lwr = cloud_optics%get_min_radius_ice() - radice_upr = cloud_optics%get_max_radius_ice() - do isub = 1,ncols_block - icol = colS + isub - 1 - do ilay = 1,LM - ! only if at least potentially cloudy ... - if (CL(icol,ilay) > 0.) then - - ! prepare for radice interpolation ... - ! first get radice consistent with RRTMGP ice cloud optics - radice = min(max(real(RR3(icol,ilay,1),kind=wp),radice_lwr),radice_upr) - ! now force into RRTMG's iceflag==3 reice binning range [5,140]um. - radice = min(max(radice,5._wp),140._wp) - ! RRTMG has 46 reice bins with 5um->radidx==1, 140um->radidx==46, - ! but radidx is forced to [1,45] so LIN2_ARG1 interpolation works. - radfac = (radice - 2._wp) / 3._wp - radidx = min(max(int(radfac),1),45) - rfint = radfac - real(radidx,kind=wp) - - do ib = 1,nbnd - ! interpolate fdelta in radice for band ib - fdelta = LIN2_ARG1(fdlice3_rrtmgp,radidx,ib,rfint) - - ! forwice calc for each g-point - do igpt = band_lims_gpt(1,ib),band_lims_gpt(2,ib) - if (cloud_props_gpt_ice%tau(isub,ilay,igpt) > 0.) then - forwice(isub,ilay,igpt) = min( & - fdelta + 0.5_wp / cloud_props_gpt_ice%ssa(isub,ilay,igpt), & - cloud_props_gpt_ice%g(isub,ilay,igpt)) - endif - enddo ! g-points - enddo ! bands - - endif ! potentially cloudy - enddo ! layers - enddo ! columns - end select - TEST_(cloud_props_gpt_ice%delta_scale(forwice)) - else - ! default delta-scaling for ice - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - forwice = cloud_props_gpt_ice%g ** 2 - end select - TEST_(cloud_props_gpt_ice%delta_scale(forwice)) - endif - endif - call MAPL_TimerOff(MAPL,"--RRTMGP_DELTA_SCALE",__RC__) + call MAPL_GetResource( MAPL, SOLCYCFRAC ,'SOLCYCFRAC:', DEFAULT=1.0, __RC__) + + ! call RRTMG SW + ! ------------- + + call RRTMG_SW (MAPL, & + RPART, NCOL, LM, & + SC, ADJES, ZT, ISOLVAR, & + PL_R, PLE_R, T_R, & + Q_R, O3_R, CO2_R, CH4_R, O2_R, & + ICEFLGSW, LIQFLGSW, & + FCLD_R, CICEWP, CLIQWP, REICE, RELIQ, & + DYOFYR, ZL_R, ALAT, & + IAER, TAUAER, SSAAER, ASMAER, & + ALBVR, ALBVF, ALBNR, ALBNF, & + LM-LCLDLM+1, LM-LCLDMH+1, NORMFLX, & + CLEARCOUNTS, SWUFLX, SWDFLX, SWUFLXC, SWDFLXC, & + NIRR, NIRF, PARR, PARF, UVRR, UVRF, FSWBAND, & + + COTDTP, COTDHP, COTDMP, COTDLP, & + COTNTP, COTNHP, COTNMP, COTNLP, & #ifdef SOLAR_RADVAL - ! REFRESH super-layer diagnostics (after delta-scaling TAUs). - ! ** Calculated from subcolumn ensemble, so stochastic ** - ! ------------------------------------------------------- - call MAPL_TimerOn(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - if (include_aerosols) then + CDSDTP, CDSDHP, CDSDMP, CDSDLP, & + CDSNTP, CDSNHP, CDSNMP, CDSNLP, & - ! in-cloud optical thicknesses in PAR super-band - ! (weighted across and within bands by TOA incident flux) - do isub = 1,ncols_block - icol = colS + isub - 1 + COTLDTP, COTLDHP, COTLDMP, COTLDLP, & + COTLNTP, COTLNHP, COTLNMP, COTLNLP, & + CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & + CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & + COTIDTP, COTIDHP, COTIDMP, COTIDLP, & + COTINTP, COTINHP, COTINMP, COTINLP, & + CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & + CDSINTP, CDSINHP, CDSINMP, CDSINLP, & - ! zero denom- and numerator accumulators - CDSDTP(icol) = 0.; CDSNTP(icol) = 0. - CDSDHP(icol) = 0.; CDSNHP(icol) = 0. - CDSDMP(icol) = 0.; CDSNMP(icol) = 0. - CDSDLP(icol) = 0.; CDSNLP(icol) = 0. - - CDSLDTP(icol) = 0.; CDSLNTP(icol) = 0.; CDSIDTP(icol) = 0.; CDSINTP(icol) = 0. - CDSLDHP(icol) = 0.; CDSLNHP(icol) = 0.; CDSIDHP(icol) = 0.; CDSINHP(icol) = 0. - CDSLDMP(icol) = 0.; CDSLNMP(icol) = 0.; CDSIDMP(icol) = 0.; CDSINMP(icol) = 0. - CDSLDLP(icol) = 0.; CDSLNLP(icol) = 0.; CDSIDLP(icol) = 0.; CDSINLP(icol) = 0. - - SDSLDTP(icol) = 0.; SDSLNTP(icol) = 0.; SDSIDTP(icol) = 0.; SDSINTP(icol) = 0. - SDSLDHP(icol) = 0.; SDSLNHP(icol) = 0.; SDSIDHP(icol) = 0.; SDSINHP(icol) = 0. - SDSLDMP(icol) = 0.; SDSLNMP(icol) = 0.; SDSIDMP(icol) = 0.; SDSINMP(icol) = 0. - SDSLDLP(icol) = 0.; SDSLNLP(icol) = 0.; SDSIDLP(icol) = 0.; SDSINLP(icol) = 0. - - ADSLDTP(icol) = 0.; ADSLNTP(icol) = 0.; ADSIDTP(icol) = 0.; ADSINTP(icol) = 0. - ADSLDHP(icol) = 0.; ADSLNHP(icol) = 0.; ADSIDHP(icol) = 0.; ADSINHP(icol) = 0. - ADSLDMP(icol) = 0.; ADSLNMP(icol) = 0.; ADSIDMP(icol) = 0.; ADSINMP(icol) = 0. - ADSLDLP(icol) = 0.; ADSLNLP(icol) = 0.; ADSIDLP(icol) = 0.; ADSINLP(icol) = 0. - - FORLDTP(icol) = 0.; FORLNTP(icol) = 0.; FORIDTP(icol) = 0.; FORINTP(icol) = 0. - FORLDHP(icol) = 0.; FORLNHP(icol) = 0.; FORIDHP(icol) = 0.; FORINHP(icol) = 0. - FORLDMP(icol) = 0.; FORLNMP(icol) = 0.; FORIDMP(icol) = 0.; FORINMP(icol) = 0. - FORLDLP(icol) = 0.; FORLNLP(icol) = 0.; FORIDLP(icol) = 0.; FORINLP(icol) = 0. + SSALDTP, SSALDHP, SSALDMP, SSALDLP, & + SSALNTP, SSALNHP, SSALNMP, SSALNLP, & + SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & + SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & + SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & + SSAINTP, SSAINHP, SSAINMP, SSAINLP, & + SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & + SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + + ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & + ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & + ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & + ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & + ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & + ASMINTP, ASMINHP, ASMINMP, ASMINLP, & + ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & + ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + + FORLDTP, FORLDHP, FORLDMP, FORLDLP, & + FORLNTP, FORLNHP, FORLNMP, FORLNLP, & + FORIDTP, FORIDHP, FORIDMP, FORIDLP, & + FORINTP, FORINHP, FORINMP, FORINLP, & +#endif + + SOLAR_TO_OBIO .and. include_aerosols, DRBAND, DFBAND, & + BAND_OUTPUT .and. include_aerosols, ISRBRGN, OSRBRGN, & + BNDSOLVAR, INDSOLVAR, SOLCYCFRAC, & + __RC__) + + call MAPL_TimerOff(MAPL,"--RRTMG_RUN") + call MAPL_TimerOn (MAPL,"--RRTMG_FLIP") + + ! unflip the outputs in the vertical + ! ---------------------------------- + + SWUFLXR (:,1:LM+1) = SWUFLX (:,LM+1:1:-1) + SWDFLXR (:,1:LM+1) = SWDFLX (:,LM+1:1:-1) + SWUFLXCR(:,1:LM+1) = SWUFLXC(:,LM+1:1:-1) + SWDFLXCR(:,1:LM+1) = SWDFLXC(:,LM+1:1:-1) + + call MAPL_TimerOff(MAPL,"--RRTMG_FLIP") - ! can only be non-zero for potentially cloudy columns - if (any(CL(icol,:) > 0.)) then + ! required outputs + ! ---------------- - ! accumulate over gpts/subcolumns - do ib = 1, nbnd - do igpt = band_lims_gpt(1,ib), band_lims_gpt(2,ib) + ! convert super-layer clearCounts to cloud fractions + if (include_aerosols) then + CLDTS(:) = 1. - CLEARCOUNTS(:,1)/float(NGPTSW) + CLDHS(:) = 1. - CLEARCOUNTS(:,2)/float(NGPTSW) + CLDMS(:) = 1. - CLEARCOUNTS(:,3)/float(NGPTSW) + CLDLS(:) = 1. - CLEARCOUNTS(:,4)/float(NGPTSW) + end if - ! band weights for photosynthetically active radiation (PAR) - ! Bands 11-12 (0.345-0.625 um) plus half transition band 10 (0.625-0.778 um) - if (ib >= 11 .and. ib <= 12) then - wgt = 1.0 - else if (ib == 10) then - wgt = 0.5 - else - ! no contribution to PAR - cycle - end if + ! undef versions of cloud optical thicknesses + ! We cannot use a merge here because some compilers (e.g. Intel) + ! will will evaluate the first argument first and if both are + ! zero, it will return NaNs. + where (COTNTP > 0. .and. COTDTP > 0.) + COTTP = COTNTP/COTDTP + elsewhere + COTTP = MAPL_UNDEF + end where - ! TOA flux weighting - ! (note: neither the adjustment of toa_flux to our tsi - ! or for zenith angle are needed yet since this weighting - ! is over gpoint and is normalized for EACH icol) - wgt = wgt * toa_flux(isub,igpt) + where (COTNHP > 0. .and. COTDHP > 0.) + COTHP = COTNHP/COTDHP + elsewhere + COTHP = MAPL_UNDEF + end where - ! low pressure layer - sltaulp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt)) - sitaulp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt)) - staulp = sltaulp + sitaulp - if (staulp > 0.) then - CDSNLP(icol) = CDSNLP(icol) + wgt * staulp - CDSDLP(icol) = CDSDLP(icol) + wgt - end if - sltaussalp = 0.; sltaussaglp = 0.; sltaussaflp = 0. - if (sltaulp > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussalp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt)) - sltaussaglp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%g (isub,LCLDLM:LM,igpt)) - sltaussaflp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & - forwliq(isub,LCLDLM:LM,igpt)) - end select - CDSLDLP(icol) = CDSLDLP(icol) + wgt - CDSLNLP(icol) = CDSLNLP(icol) + wgt * sltaulp - SDSLDLP(icol) = SDSLDLP(icol) + wgt * sltaulp - SDSLNLP(icol) = SDSLNLP(icol) + wgt * sltaussalp - ADSLDLP(icol) = ADSLDLP(icol) + wgt * sltaussalp - ADSLNLP(icol) = ADSLNLP(icol) + wgt * sltaussaglp - FORLDLP(icol) = FORLDLP(icol) + wgt * sltaussalp - FORLNLP(icol) = FORLNLP(icol) + wgt * sltaussaflp - end if - sitaussalp = 0.; sitaussaglp = 0.; sitaussaflp = 0. - if (sitaulp > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussalp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt)) - sitaussaglp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%g (isub,LCLDLM:LM,igpt)) - sitaussaflp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & - forwice(isub,LCLDLM:LM,igpt)) - end select - CDSIDLP(icol) = CDSIDLP(icol) + wgt - CDSINLP(icol) = CDSINLP(icol) + wgt * sitaulp - SDSIDLP(icol) = SDSIDLP(icol) + wgt * sitaulp - SDSINLP(icol) = SDSINLP(icol) + wgt * sitaussalp - ADSIDLP(icol) = ADSIDLP(icol) + wgt * sitaussalp - ADSINLP(icol) = ADSINLP(icol) + wgt * sitaussaglp - FORIDLP(icol) = FORIDLP(icol) + wgt * sitaussalp - FORINLP(icol) = FORINLP(icol) + wgt * sitaussaflp - end if + where (COTNMP > 0. .and. COTDMP > 0.) + COTMP = COTNMP/COTDMP + elsewhere + COTMP = MAPL_UNDEF + end where - ! mid pressure layer - sltaump = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt)) - sitaump = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt)) - staump = sltaump + sitaump - if (staump > 0.) then - CDSNMP(icol) = CDSNMP(icol) + wgt * staump - CDSDMP(icol) = CDSDMP(icol) + wgt - end if - sltaussamp = 0.; sltaussagmp = 0.; sltaussafmp = 0. - if (sltaump > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussamp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt)) - sltaussagmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%g (isub,LCLDMH:LCLDLM-1,igpt)) - sltaussafmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - forwliq(isub,LCLDMH:LCLDLM-1,igpt)) - end select - CDSLDMP(icol) = CDSLDMP(icol) + wgt - CDSLNMP(icol) = CDSLNMP(icol) + wgt * sltaump - SDSLDMP(icol) = SDSLDMP(icol) + wgt * sltaump - SDSLNMP(icol) = SDSLNMP(icol) + wgt * sltaussamp - ADSLDMP(icol) = ADSLDMP(icol) + wgt * sltaussamp - ADSLNMP(icol) = ADSLNMP(icol) + wgt * sltaussagmp - FORLDMP(icol) = FORLDMP(icol) + wgt * sltaussamp - FORLNMP(icol) = FORLNMP(icol) + wgt * sltaussafmp - end if - sitaussamp = 0.; sitaussagmp = 0.; sitaussafmp = 0. - if (sitaump > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussamp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt)) - sitaussagmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%g (isub,LCLDMH:LCLDLM-1,igpt)) - sitaussafmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & - forwice(isub,LCLDMH:LCLDLM-1,igpt)) - end select - CDSIDMP(icol) = CDSIDMP(icol) + wgt - CDSINMP(icol) = CDSINMP(icol) + wgt * sitaump - SDSIDMP(icol) = SDSIDMP(icol) + wgt * sitaump - SDSINMP(icol) = SDSINMP(icol) + wgt * sitaussamp - ADSIDMP(icol) = ADSIDMP(icol) + wgt * sitaussamp - ADSINMP(icol) = ADSINMP(icol) + wgt * sitaussagmp - FORIDMP(icol) = FORIDMP(icol) + wgt * sitaussamp - FORINMP(icol) = FORINMP(icol) + wgt * sitaussafmp - end if + where (COTNLP > 0. .and. COTDLP > 0.) + COTLP = COTNLP/COTDLP + elsewhere + COTLP = MAPL_UNDEF + end where - ! high pressure layer - sltauhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt)) - sitauhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt)) - stauhp = sltauhp + sitauhp - if (stauhp > 0.) then - CDSNHP(icol) = CDSNHP(icol) + wgt * stauhp - CDSDHP(icol) = CDSDHP(icol) + wgt - end if - sltaussahp = 0.; sltaussaghp = 0.; sltaussafhp = 0. - if (sltauhp > 0.) then - select type(cloud_props_gpt_liq) - class is (ty_optical_props_2str) - sltaussahp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt)) - sltaussaghp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%g (isub,1:LCLDMH-1,igpt)) - sltaussafhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & - forwliq(isub,1:LCLDMH-1,igpt)) - end select - CDSLDHP(icol) = CDSLDHP(icol) + wgt - CDSLNHP(icol) = CDSLNHP(icol) + wgt * sltauhp - SDSLDHP(icol) = SDSLDHP(icol) + wgt * sltauhp - SDSLNHP(icol) = SDSLNHP(icol) + wgt * sltaussahp - ADSLDHP(icol) = ADSLDHP(icol) + wgt * sltaussahp - ADSLNHP(icol) = ADSLNHP(icol) + wgt * sltaussaghp - FORLDHP(icol) = FORLDHP(icol) + wgt * sltaussahp - FORLNHP(icol) = FORLNHP(icol) + wgt * sltaussafhp - end if - sitaussahp = 0.; sitaussaghp = 0.; sitaussafhp = 0. - if (sitauhp > 0.) then - select type(cloud_props_gpt_ice) - class is (ty_optical_props_2str) - sitaussahp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt)) - sitaussaghp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%g (isub,1:LCLDMH-1,igpt)) - sitaussafhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & - cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & - forwice(isub,1:LCLDMH-1,igpt)) - end select - CDSIDHP(icol) = CDSIDHP(icol) + wgt - CDSINHP(icol) = CDSINHP(icol) + wgt * sitauhp - SDSIDHP(icol) = SDSIDHP(icol) + wgt * sitauhp - SDSINHP(icol) = SDSINHP(icol) + wgt * sitaussahp - ADSIDHP(icol) = ADSIDHP(icol) + wgt * sitaussahp - ADSINHP(icol) = ADSINHP(icol) + wgt * sitaussaghp - FORIDHP(icol) = FORIDHP(icol) + wgt * sitaussahp - FORINHP(icol) = FORINHP(icol) + wgt * sitaussafhp - end if +#ifdef SOLAR_RADVAL + ! zero versions of cloud optical thicknesses + ! We can use merge() here because we cannot divide by zero + TAUTP = merge(COTTP, 0., COTNTP > 0. .and. COTDTP > 0.) + TAUHP = merge(COTHP, 0., COTNHP > 0. .and. COTDHP > 0.) + TAUMP = merge(COTMP, 0., COTNMP > 0. .and. COTDMP > 0.) + TAULP = merge(COTLP, 0., COTNLP > 0. .and. COTDLP > 0.) +#endif - ! whole subcolumn - sltautp = sltaulp + sltaump + sltauhp - sitautp = sitaulp + sitaump + sitauhp - stautp = staulp + staump + stauhp - if (stautp > 0.) then - CDSNTP(icol) = CDSNTP(icol) + wgt * stautp - CDSDTP(icol) = CDSDTP(icol) + wgt - end if - sltaussatp = sltaussalp + sltaussamp + sltaussahp - sltaussagtp = sltaussaglp + sltaussagmp + sltaussaghp - sltaussaftp = sltaussaflp + sltaussafmp + sltaussafhp - if (sltautp > 0.) then - CDSLDTP(icol) = CDSLDTP(icol) + wgt - CDSLNTP(icol) = CDSLNTP(icol) + wgt * sltautp - SDSLDTP(icol) = SDSLDTP(icol) + wgt * sltautp - SDSLNTP(icol) = SDSLNTP(icol) + wgt * sltaussatp - ADSLDTP(icol) = ADSLDTP(icol) + wgt * sltaussatp - ADSLNTP(icol) = ADSLNTP(icol) + wgt * sltaussagtp - FORLDTP(icol) = FORLDTP(icol) + wgt * sltaussatp - FORLNTP(icol) = FORLNTP(icol) + wgt * sltaussaftp - end if - sitaussatp = sitaussalp + sitaussamp + sitaussahp - sitaussagtp = sitaussaglp + sitaussagmp + sitaussaghp - sitaussaftp = sitaussaflp + sitaussafmp + sitaussafhp - if (sitautp > 0.) then - CDSIDTP(icol) = CDSIDTP(icol) + wgt - CDSINTP(icol) = CDSINTP(icol) + wgt * sitautp - SDSIDTP(icol) = SDSIDTP(icol) + wgt * sitautp - SDSINTP(icol) = SDSINTP(icol) + wgt * sitaussatp - ADSIDTP(icol) = ADSIDTP(icol) + wgt * sitaussatp - ADSINTP(icol) = ADSINTP(icol) + wgt * sitaussagtp - FORIDTP(icol) = FORIDTP(icol) + wgt * sitaussatp - FORINTP(icol) = FORINTP(icol) + wgt * sitaussaftp - end if + ! fluxes + FSW = SWDFLXR - SWUFLXR + FSC = SWDFLXCR - SWUFLXCR + FSWU = SWUFLXR + FSCU = SWUFLXCR - end do ! igpt - end do ! ib + ! Deallocate the working inputs + !------------------------------ + deallocate(TLEV ,__STAT__) + deallocate(TLEV_R,__STAT__) + deallocate(PLE_R ,__STAT__) + deallocate(FCLD_R,__STAT__) + deallocate(CLIQWP,__STAT__) + deallocate(CICEWP,__STAT__) + deallocate(RELIQ ,__STAT__) + deallocate(REICE ,__STAT__) - end if ! potentially cloudy column - end do ! isub - end if ! include_aerosols - call MAPL_TimerOff(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) -#endif + deallocate(TAUAER,__STAT__) + deallocate(SSAAER,__STAT__) + deallocate(ASMAER,__STAT__) + deallocate(DPR ,__STAT__) + deallocate(PL_R ,__STAT__) + deallocate(ZL_R ,__STAT__) + deallocate(T_R ,__STAT__) + deallocate(Q_R ,__STAT__) + deallocate(O2_R ,__STAT__) + deallocate(O3_R ,__STAT__) + deallocate(CO2_R ,__STAT__) + deallocate(CH4_R ,__STAT__) - call MAPL_TimerOn(MAPL,"--RRTMGP_RT",__RC__) + deallocate(CLEARCOUNTS ,__STAT__) - ! scale to our tsi - ! (both toa_flux and tsi are NORMAL to solar beam, [W/m2]) - toa_flux = toa_flux * spread(tsi(colS:colE)/sum(toa_flux,dim=2), 2, ngpt) + deallocate(SWUFLX ,__STAT__) + deallocate(SWDFLX ,__STAT__) + deallocate(SWUFLXC,__STAT__) + deallocate(SWDFLXC,__STAT__) - ! add in aerosol optical properties if requested and available - if (need_aer_optical_props) then - TEST_(aer_props%increment(optical_props)) - end if + deallocate(SWUFLXR ,__STAT__) + deallocate(SWDFLXR ,__STAT__) + deallocate(SWUFLXCR,__STAT__) + deallocate(SWDFLXCR,__STAT__) - ! clear-sky radiative transfer - fluxes_clrsky%flux_up => flux_up_clrsky(colS:colE,:) - fluxes_clrsky%flux_net => flux_net_clrsky(colS:colE,:) - error_msg = rte_sw( & - optical_props, top_at_1, mu0(colS:colE), toa_flux, & - sfc_alb_dir(:,colS:colE), sfc_alb_dif(:,colS:colE), & - fluxes_clrsky) - TEST_(error_msg) - - ! add in cloud optical properties - ! add ice first since its optical depths are usually smaller - TEST_(cloud_props_gpt_ice%increment(optical_props)) - TEST_(cloud_props_gpt_liq%increment(optical_props)) - - ! all-sky radiative transfer - fluxes_allsky%flux_up => flux_up_allsky(colS:colE,:) - fluxes_allsky%flux_net => flux_net_allsky(colS:colE,:) - fluxes_allsky%bnd_flux_dn => bnd_flux_dn_allsky(colS:colE,:,:) - fluxes_allsky%bnd_flux_dn_dir => bnd_flux_dir_allsky(colS:colE,:,:) - fluxes_allsky%bnd_flux_net => bnd_flux_net_allsky(colS:colE,:,:) - error_msg = rte_sw( & - optical_props, top_at_1, mu0(colS:colE), toa_flux, & - sfc_alb_dir(:,colS:colE), sfc_alb_dif(:,colS:colE), & - fluxes_allsky) - TEST_(error_msg) - - call MAPL_TimerOff(MAPL,"--RRTMGP_RT",__RC__) + call MAPL_TimerOff(MAPL,"-RRTMG") - end do ! loop over blocks + else - call MAPL_TimerOn(MAPL,"--RRTMGP_POST",__RC__) + _FAIL('unknown SW radiation scheme!') - ! normalize by incoming solar radiation - allocate(flux_dn_top(ncol), __STAT__) - flux_dn_top(:) = real(max(SLR1D, 1e-7), kind=wp) - do k = 1, LM+1 - flux_up_clrsky (:,k) = flux_up_clrsky (:,k) / flux_dn_top(:) - flux_net_clrsky(:,k) = flux_net_clrsky(:,k) / flux_dn_top(:) - end do - do k = 1, LM+1 - flux_up_allsky (:,k) = flux_up_allsky (:,k) / flux_dn_top(:) - flux_net_allsky(:,k) = flux_net_allsky(:,k) / flux_dn_top(:) - end do - do ib = 1, nbnd - do k = 1, LM+1 - bnd_flux_dn_allsky (:,k,ib) = bnd_flux_dn_allsky (:,k,ib) / flux_dn_top(:) - bnd_flux_net_allsky(:,k,ib) = bnd_flux_net_allsky(:,k,ib) / flux_dn_top(:) - bnd_flux_dir_allsky(:,k,ib) = bnd_flux_dir_allsky(:,k,ib) / flux_dn_top(:) - end do - end do - deallocate(flux_dn_top, __STAT__) + end if SCHEME - ! load output arrays - ! clear-sky fluxes - FSCU = real(flux_up_clrsky) - FSC = real(flux_net_clrsky) - ! all-sky fluxes - FSWU = real(flux_up_allsky) - FSW = real(flux_net_allsky) - ! surface net flux per band - do ib = 1, nbnd - FSWBAND(:,ib) = real(bnd_flux_net_allsky(:,LM+1,ib)) - end do - ! surface downwelling direct and diffuse fluxes in bands - if (SOLAR_TO_OBIO .and. include_aerosols) then - do ib = 1, nbnd - DRBAND(:,ib) = real(bnd_flux_dir_allsky(:,LM+1,ib)) - DFBAND(:,ib) = real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) - end do - endif - ! TOA band fluxes - if (include_aerosols) then - if (USE_RRTMG .or. USE_RRTMGP) then - do ib = 1, nbnd - if (band_output(ib)) then - ISRBRGN(ib) % p = real(bnd_flux_dn_allsky(:,1,ib)) - OSRBRGN(ib) % p = real(bnd_flux_dn_allsky(:,1,ib) - bnd_flux_net_allsky(:,1,ib)) - end if - end do - end if - end if + ! Deallocate the working inputs + !------------------------------ - ! surface direct and diffuse downward in super-bands - ! for *diffuse* downward must subtract direct (downward) from total downward - ! pmn: may later do this using a flux class extension?? + deallocate (PL, RH, PLhPa) + deallocate (QQ3, RR3) + deallocate (ILWT) + deallocate (O3) + deallocate (TAUA, SSAA, ASYA) - ! NIR bands (1-9: 820-12850 cm-1, 0.778-12.195 microns) - NIRR = 0.; NIRF = 0. - do ib=1,9 - NIRR = NIRR + real(bnd_flux_dir_allsky(:,LM+1,ib)) - NIRF = NIRF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) - end do - ! PAR bands (11-12: 16000-29000 cm-1, 0.345-0.625 micron) - PARR = 0.; PARF = 0. - do ib=11,12 - PARR = PARR + real(bnd_flux_dir_allsky(:,LM+1,ib)) - PARF = PARF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) - enddo - ! UVR bands (13-14: 29000-50000 cm-1, 0.200-0.345 micron) - UVRR = 0.; UVRF = 0. - do ib=13,14 - UVRR = UVRR + real(bnd_flux_dir_allsky(:,LM+1,ib)) - UVRF = UVRF + real(bnd_flux_dn_allsky (:,LM+1,ib) - bnd_flux_dir_allsky(:,LM+1,ib)) - enddo - ! Transition band (10, 12850-16000 cm-1, 0.625-0.778 micron) - ! split half-and-half to PAR and NIR - NIRR = NIRR + 0.5 * real(bnd_flux_dir_allsky(:,LM+1,10)) - PARR = PARR + 0.5 * real(bnd_flux_dir_allsky(:,LM+1,10)) - NIRF = NIRF + 0.5 * real(bnd_flux_dn_allsky (:,LM+1,10) - bnd_flux_dir_allsky(:,LM+1,10)) - PARF = PARF + 0.5 * real(bnd_flux_dn_allsky (:,LM+1,10) - bnd_flux_dir_allsky(:,LM+1,10)) + ! Complete load balancing by retrieving work done remotely + !--------------------------------------------------------- - ! clean up - deallocate(band_lims_gpt,__STAT__) - deallocate(tsi,mu0,sfc_alb_dir,sfc_alb_dif,toa_flux,__STAT__) - deallocate(dummy_wp,p_lay,t_lay,p_lev,dp_wp,dzmid,__STAT__) - deallocate(flux_up_clrsky,flux_net_clrsky,__STAT__) - deallocate(flux_up_allsky,flux_net_allsky,__STAT__) - deallocate(bnd_flux_dn_allsky,bnd_flux_net_allsky,bnd_flux_dir_allsky,__STAT__) - deallocate(seeds,urand,cld_mask,__STAT__) - deallocate(forwliq,forwice,__STAT__) - if (gen_mro) then - deallocate(adl,alpha,urand_aux,__STAT__) - if (cond_inhomo) then - deallocate(rdl,rcorr,urand_cond,urand_cond_aux,zcw,__STAT__) - endif - end if - if (include_aerosols) deallocate(ClearCounts,__STAT__) - call cloud_optics%finalize() - call cloud_props_gpt_liq%finalize() - call cloud_props_gpt_ice%finalize() - call cloud_props_bnd_liq%finalize() - call cloud_props_bnd_ice%finalize() - if (need_aer_optical_props) call aer_props%finalize() - call optical_props%finalize() + call MAPL_TimerOn(MAPL,"-BALANCE") - call MAPL_TimerOff(MAPL,"--RRTMGP_POST",__RC__) + call MAPL_TimerOn(MAPL,"--RETRIEVE") + if (LoadBalance) then + if (size(BufOut) > 0) call MAPL_BalanceWork(BufOut, NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + if (size(BufInOut) > 0) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + end if + call MAPL_TimerOff(MAPL,"--RETRIEVE") - call MAPL_TimerOff(MAPL,"-RRTMGP",__RC__) + ! Unpack the results. Fills "masked" (night) locations with default value from internal state + !-------------------------------------------------------------------------------------------- + ! resulting internals are then contiguous versions + ! Note: InOut variables do not fill unmasked locations with a default, + ! since the unmasked locations may contain potentially useful aged data. -#undef TEST_ + i1InOut = 1; i1Out = 1 + INT_VARS_3: do k=1,NumInt + if (SlicesInt(k) == 0) cycle - else if (USE_RRTMG) then + if (IntInOut(k)) then + pi1 => i1InOut + else + pi1 => i1Out + call MAPL_VarSpecGet(InternalSpec(k),DEFAULT=def,__RC__) + endif - ! regular RRTMG - call MAPL_TimerOn(MAPL,"-RRTMG") + if (ugDim(k) > 0) then + select case(rgDim(k)) + case(MAPL_DIMSHORZVERT) + call ESMFL_StateGetPointerToData(INTERNAL,ptr4,NamesInt(k),__RC__) + if (IntInOut(k)) then + do j=1,ugDim(k) + call UnPackIt(BufInOut(pi1+(j-1)*size(ptr4,3)*NumMax),ptr4(:,:,:,j), & + daytime,NumMax,HorzDims,size(ptr4,3)) + end do + else + do j=1,ugDim(k) + call UnPackIt(BufOut (pi1+(j-1)*size(ptr4,3)*NumMax),ptr4(:,:,:,j), & + daytime,NumMax,HorzDims,size(ptr4,3),def) + end do + endif + case(MAPL_DIMSHORZONLY) + call ESMFL_StateGetPointerToData(INTERNAL,ptr3,NamesInt(k),__RC__) + if (IntInOut(k)) then + call UnPackIt(BufInOut(pi1),ptr3,daytime,NumMax,HorzDims,ugDim(k)) + else + call UnPackIt(BufOut (pi1),ptr3,daytime,NumMax,HorzDims,ugDim(k),def) + endif + end select + else + select case(rgDim(k)) + case(MAPL_DIMSHORZVERT) + call ESMFL_StateGetPointerToData(INTERNAL,ptr3,NamesInt(k),__RC__) + if (IntInOut(k)) then + call UnPackIt(BufInOut(pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3)) + else + call UnPackIt(BufOut (pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3),def) + endif + case(MAPL_DIMSHORZONLY) + call ESMFL_StateGetPointerToData(INTERNAL,ptr2,NamesInt(k),__RC__) + if (IntInOut(k)) then + call UnPackIt(BufInOut(pi1),ptr2,daytime,NumMax,HorzDims,1) + else + call UnPackIt(BufOut (pi1),ptr2,daytime,NumMax,HorzDims,1,def) + endif + end select + end if + pi1 = pi1 + NumMax*SlicesInt(k) - ! reversed (flipped) vertical dimension arrays and other RRTMG arrays - ! ------------------------------------------------------------------- + enddo INT_VARS_3 - ! interface (between layer) variables - allocate(TLEV (NCOL,LM+1),__STAT__) - allocate(TLEV_R(NCOL,LM+1),__STAT__) - allocate(PLE_R (NCOL,LM+1),__STAT__) - ! cloud physical properties - allocate(FCLD_R(NCOL,LM),__STAT__) - allocate(CLIQWP(NCOL,LM),__STAT__) - allocate(CICEWP(NCOL,LM),__STAT__) - allocate(RELIQ (NCOL,LM),__STAT__) - allocate(REICE (NCOL,LM),__STAT__) - ! aerosol optical properties - allocate(TAUAER(NCOL,LM,NB_RRTMG),__STAT__) - allocate(SSAAER(NCOL,LM,NB_RRTMG),__STAT__) - allocate(ASMAER(NCOL,LM,NB_RRTMG),__STAT__) - ! layer variables - allocate(DPR (NCOL,LM),__STAT__) - allocate(PL_R (NCOL,LM),__STAT__) - allocate(ZL_R (NCOL,LM),__STAT__) - allocate(T_R (NCOL,LM),__STAT__) - allocate(Q_R (NCOL,LM),__STAT__) - allocate(O2_R (NCOL,LM),__STAT__) - allocate(O3_R (NCOL,LM),__STAT__) - allocate(CO2_R (NCOL,LM),__STAT__) - allocate(CH4_R (NCOL,LM),__STAT__) - ! super-layer cloud fractions - allocate(CLEARCOUNTS (NCOL,4),__STAT__) - ! output fluxes - allocate(SWUFLX (NCOL,LM+1),__STAT__) - allocate(SWDFLX (NCOL,LM+1),__STAT__) - allocate(SWUFLXC(NCOL,LM+1),__STAT__) - allocate(SWDFLXC(NCOL,LM+1),__STAT__) - ! un-flipped outputs - allocate(SWUFLXR (NCOL,LM+1),__STAT__) - allocate(SWDFLXR (NCOL,LM+1),__STAT__) - allocate(SWUFLXCR(NCOL,LM+1),__STAT__) - allocate(SWDFLXCR(NCOL,LM+1),__STAT__) + ! clean up + deallocate(SlicesInp,NamesInp,__STAT__) + deallocate(SlicesInt,NamesInt,__STAT__) + deallocate(IntInOut,rgDim,ugDim,__STAT__) + deallocate(BufInp,BufInOut,BufOut,__STAT__) + call MAPL_TimerOn(MAPL,"--DESTROY") + if (LoadBalance) call MAPL_BalanceDestroy(Handle=SolarBalanceHandle, __RC__) + call MAPL_TimerOff(MAPL,"--DESTROY") - ! Set flags related to cloud properties (see RRTMG_SW) - ! ---------------------------------------------------- - call MAPL_GetResource(MAPL,ICEFLGSW,'RRTMG_ICEFLG:',DEFAULT=3,__RC__) - call MAPL_GetResource(MAPL,LIQFLGSW,'RRTMG_LIQFLG:',DEFAULT=1,__RC__) + call MAPL_TimerOff(MAPL,"-BALANCE") - if (LM > 72) then - call MAPL_GetResource(MAPL,USE_PRECIP_IN_RADIATION,'RRTMGSW_USE_PRECIP_IN_RADIATION:',DEFAULT=.TRUE.,RC=STATUS) - VERIFY_(STATUS) - else - call MAPL_GetResource(MAPL,USE_PRECIP_IN_RADIATION,'RRTMGSW_USE_PRECIP_IN_RADIATION:',DEFAULT=.FALSE.,RC=STATUS) - VERIFY_(STATUS) - endif + RETURN_(ESMF_SUCCESS) - ! Normalize aerosol inputs - ! ------------------------ + end subroutine SORADCORE - if (num_aero_vars > 0) then - where (TAUA > 0. .and. SSAA > 0. ) - ASYA = ASYA/SSAA - SSAA = SSAA/TAUA - elsewhere - TAUA = 0. - SSAA = 0. - ASYA = 0. - end where - end if + ! --------------------------------------------------------------------------- + ! Compute gas optical properties for one block of columns. + ! gas_concs_block and error_msg are local (thread-private in future OMP use). + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_gas_optics(colS, colE, ncols_block, LM, & + gas_concs, k_dist, p_lay, p_lev, t_lay, & + optical_props, toa_flux, MAPL, RC) - ! Flip in vertical, Convert units, and interpolate T, etc. - ! -------------------------------------------------------- - ! RRTMG convention is that vertical indices increase from bot -> top + use mo_gas_concentrations, only: ty_gas_concs + use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_2str + use mo_rte_kind, only: wp + + integer, intent(in) :: colS, colE, ncols_block, LM + type(ty_gas_concs), intent(in) :: gas_concs + type(ty_gas_optics_rrtmgp), intent(in) :: k_dist + real(wp), intent(in) :: p_lay(:,:), p_lev(:,:), t_lay(:,:) + class(ty_optical_props_arry), intent(inout) :: optical_props + real(wp), intent(out) :: toa_flux(:,:) + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! locals -- will be thread-private under future !$OMP PARALLEL DO + type(ty_gas_concs) :: gas_concs_block + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS + + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_GAS_OPTICS",__RC__) + + TEST_(gas_concs%get_subset(colS, ncols_block, gas_concs_block)) + + ! gas optics, including source functions + error_msg = k_dist%gas_optics( & + p_lay(colS:colE,:), p_lev(colS:colE,:), t_lay(colS:colE,:), & + gas_concs_block, optical_props, toa_flux) + TEST_(error_msg) - call MAPL_TimerOn (MAPL,"--RRTMG_FLIP") + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_GAS_OPTICS",__RC__) - DPR(:,1:LM) = (PLE(:,2:LM+1)-PLE(:,1:LM)) + RETURN_(ESMF_SUCCESS) - ! cloud water paths converted from g/g to g/m^2 - ! Flip in vertical - do K = 1,LM - LV = LM-K+1 ! LM --> 1 + end subroutine compute_gas_optics +#undef TEST_ - ! Convert content [kg/kg] to path [g/m2] - ! using hydrostatic eqn dp/g ~ rho*dz, - ! so conversion factor is 1000*dp/g ~ 1.02*100*dp. - ! pmn: why not use MAPL_GRAV explicitly? - if (USE_PRECIP_IN_RADIATION) then - ILWT = QQ3(:,LV,KLIQUID)+QQ3(:,LV,KRAIN) - CLIQWP(:,K) = 1.02*100*DPR(:,LV)*ILWT - where (ILWT > 0.0) - RELIQ (:,K) = ( RR3(:,LV,KLIQUID)*QQ3(:,LV,KLIQUID) + & - RR3(:,LV,KRAIN )*QQ3(:,LV,KRAIN ) ) / ILWT - elsewhere - RELIQ (:,K) = 14.0 - end where - ILWT = QQ3(:,LV,KICE)+QQ3(:,LV,KSNOW)+QQ3(:,LV,KGRAUPEL) - CICEWP(:,K) = 1.02*100*DPR(:,LV)*ILWT - WHERE (ILWT > 0.0) - REICE (:,K) = ( RR3(:,LV,KICE )*QQ3(:,LV,KICE ) + & - RR3(:,LV,KSNOW )*QQ3(:,LV,KSNOW ) + & - RR3(:,LV,KGRAUPEL)*QQ3(:,LV,KGRAUPEL) ) / ILWT - elsewhere - REICE (:,K) = 36.0 - end where - else - CLIQWP(:,K) = 1.02*100*DPR(:,LV)*QQ3(:,LV,KLIQUID) - CICEWP(:,K) = 1.02*100*DPR(:,LV)*QQ3(:,LV,KICE) - RELIQ (:,K) = RR3(:,LV,KLIQUID) - REICE (:,K) = RR3(:,LV,KICE ) - endif - enddo + ! --------------------------------------------------------------------------- + ! Load and normalize aerosol optical properties for one block of columns. + ! No thread-private locals beyond scalars; aer_props is intent(inout). + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_aer_optics(colS, colE, need_aer_optical_props, & + taua, ssaa, asya, aer_props, RC) - IF (ICEFLGSW == 0) THEN - WHERE (REICE < 10.) REICE = 10. - WHERE (REICE > 30.) REICE = 30. - ELSE IF (ICEFLGSW == 1) THEN - WHERE (REICE < 13.) REICE = 13. - WHERE (REICE > 130.) REICE = 130. - ELSE IF (ICEFLGSW == 2) THEN - WHERE (REICE < 5.) REICE = 5. - WHERE (REICE > 131.) REICE = 131. - ELSE IF (ICEFLGSW == 3) THEN - WHERE (REICE < 5.) REICE = 5. - WHERE (REICE > 140.) REICE = 140. - ELSE IF (ICEFLGSW == 4) THEN - REICE(:,:) = REICE(:,:)*2. - WHERE (REICE < 1.) REICE = 1. - WHERE (REICE > 200.) REICE = 200. - END IF + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_2str + use mo_rte_kind, only: wp - IF (LIQFLGSW == 0) THEN - WHERE (RELIQ < 10.) RELIQ = 10. - WHERE (RELIQ > 30.) RELIQ = 30. - ELSE IF (LIQFLGSW == 1) THEN - WHERE (RELIQ < 2.5) RELIQ = 2.5 - WHERE (RELIQ > 60.) RELIQ = 60. - END IF + integer, intent(in) :: colS, colE + logical, intent(in) :: need_aer_optical_props + real, dimension(:,:,:), intent(in) :: taua, ssaa, asya + class(ty_optical_props_arry), intent(inout) :: aer_props + integer, optional, intent(out) :: RC - ! regular (non-flipped) interface temperatures - TLEV(:,2:LM)=(T(:,1:LM-1)* DPR(:,2:LM) + T(:,2:LM) * DPR(:,1:LM-1)) & - /(DPR(:,1:LM-1) + DPR(:,2:LM)) - TLEV(:,LM+1) = TS(:) - TLEV(:, 1) = TLEV(:,2) + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS - ! flip in vertical ... + if (.not. need_aer_optical_props) then + RETURN_(ESMF_SUCCESS) + end if - PLE_R (:,1:LM+1) = PLE (:,LM+1:1:-1) / 100. ! hPa - TLEV_R(:,1:LM+1) = TLEV(:,LM+1:1:-1) + select type (aer_props) + class is (ty_optical_props_2str) - PL_R (:,1:LM ) = PL (:,LM:1:-1) / 100. ! hPa - T_R (:,1:LM ) = T (:,LM:1:-1) + ! load un-normalized optical properties from aerosol system + aer_props%tau = real(taua(colS:colE,:,:),kind=wp) + aer_props%ssa = real(ssaa(colS:colE,:,:),kind=wp) + aer_props%g = real(asya(colS:colE,:,:),kind=wp) + + ! renormalize + where (aer_props%tau > 0._wp .and. aer_props%ssa > 0._wp) + aer_props%g = aer_props%g / aer_props%ssa + aer_props%ssa = aer_props%ssa / aer_props%tau + elsewhere + aer_props%tau = 0._wp + aer_props%ssa = 0._wp + aer_props%g = 0._wp + end where + + ! Because RRTMGP is (currently) compiled at R8, + ! _wp is R8. Apparently with aggressive compiler + ! flags using Intel, it's possible for, say, + ! aer_props%ssa to become slightly greater than one + ! in the above renormalization. So, we add clamps + ! to the values based on the restrictions see in + ! RRTMGP/rte-frontend/mo_optical_props.F90 + ! + ! In testing, the values seen were like 1.00000011905028 + ! so just slightly above one. + + ! tau must be greater than 0.0 + aer_props%tau = max(aer_props%tau, 0._wp) + ! ssa must be between 0.0 and 1.0 + aer_props%ssa = max(min(aer_props%ssa, 1._wp), 0._wp) + ! g must be between -1.0 and 1.0 + aer_props%g = max(min(aer_props%g, 1._wp),-1._wp) - ! Specific humidity is converted to Volume Mixing Ratio - Q_R (:,1:LM ) = Q (:,LM:1:-1) / (1.-Q(:,LM:1:-1)) * (MAPL_AIRMW/MAPL_H2OMW) + class default + TEST_('aerosol optical properties hardwired 2-stream for now') + end select - ! Ozone is converted Mass Mixing Ratio to Volume Mixing Ratio - O3_R (:,1:LM ) = O3 (:,LM:1:-1) * (MAPL_AIRMW/MAPL_O3MW) + RETURN_(ESMF_SUCCESS) - ! chemistry and cloud fraction - ! (cloud water paths and effective radii flipped already) - CH4_R (:,1:LM ) = CH4(:,LM:1:-1) - CO2_R (:,1:LM ) = CO2 - O2_R (:,1:LM ) = O2 - FCLD_R(:,1:LM ) = CL (:,LM:1:-1) + end subroutine compute_aer_optics +#undef TEST_ -! Clean up negatives - WHERE (Q_R < 0.) Q_R = 0. - WHERE (O3_R < 0.) O3_R = 0. - WHERE (CH4_R < 0.) CH4_R = 0. - WHERE (CO2_R < 0.) CO2_R = 0. - WHERE (O2_R < 0.) O2_R = 0. - WHERE (FCLD_R < 0.) FCLD_R = 0. + ! --------------------------------------------------------------------------- + ! Compute cloud optical properties and McICA sampling for one block. + ! All per-block temporaries (alpha, rcorr, zcw, urand*, cld_mask, rng, + ! seeds) are local -- thread-private in future !$OMP PARALLEL DO use. + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_cloud_optics_mcica( & + colS, colE, ncols_block, LM, ngpt, & + gen_mro, cond_inhomo, cloud_overlap_type, cwp_fac, IM_World, & + seeds_time_key, & + QQ3, RR3, dp_wp, dummy_wp, CL, dzmid, adl, rdl, Ig1D, Jg1D, & + cloud_optics, & + cloud_props_bnd_liq, cloud_props_bnd_ice, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + cld_mask, & + MAPL, RC) - ! Adjustment for Earth/Sun distance, from MAPL_SunGetInsolation - ADJES = DIST + use mo_cloud_sampling, only: draw_samples, & + sampled_mask_max_ran, & + sampled_urand_gen_max_ran + use mo_cloud_optics_rrtmgp, only: ty_cloud_optics_rrtmgp + use mo_optical_props, only: ty_optical_props_arry + use mo_rte_kind, only: wp + use cloud_condensate_inhomogeneity, only: zcw_lookup +#ifdef HAVE_MKL + use MKL_VSL_TYPE + use mo_rng_mklvsl_plus, only: ty_rng_mklvsl_plus +#else + use mo_rng_mt19937, only: ty_rng_mt +#endif - ! Layer mid-point heights relative to zero at index 1 - ZL_R(:,1) = 0. - do k=2,LM - ! dz = RT/g x dp/p - ! Note: This is correct even though its different from LW. - ! Its because SW uses LE[V]_R 1:LM+1 while LW uses 0:LM. - ZL_R(:,k) = ZL_R(:,k-1) + MAPL_RGAS*TLEV_R(:,k)/MAPL_GRAV*(PL_R(:,k-1)-PL_R(:,k))/PLE_R(:,k) - enddo + integer, intent(in) :: colS, colE, ncols_block, LM, ngpt + logical, intent(in) :: gen_mro, cond_inhomo + character(len=*), intent(in) :: cloud_overlap_type + real(wp), intent(in) :: cwp_fac + integer, intent(in) :: IM_World + integer, intent(in) :: seeds_time_key + real, dimension(:,:,:), intent(in) :: QQ3, RR3 + real(wp), dimension(:,:), intent(in) :: dp_wp, dummy_wp + real, dimension(:,:), intent(in) :: CL + real(wp), dimension(:,:), intent(in) :: dzmid + real, dimension(:), intent(in) :: adl, rdl + real, dimension(:), intent(in) :: Ig1D, Jg1D + type(ty_cloud_optics_rrtmgp), intent(inout) :: cloud_optics + class(ty_optical_props_arry), intent(inout) :: cloud_props_bnd_liq, cloud_props_bnd_ice + class(ty_optical_props_arry), intent(inout) :: cloud_props_gpt_liq, cloud_props_gpt_ice + logical, allocatable, intent(out) :: cld_mask(:,:,:) + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! locals -- all thread-private under future !$OMP PARALLEL DO + real(wp), allocatable :: alpha(:,:), rcorr(:,:) + real(wp), allocatable :: zcw(:,:,:) + real(wp), allocatable :: urand(:,:,:), urand_aux(:,:,:) + real(wp), allocatable :: urand_cond(:,:,:), urand_cond_aux(:,:,:) + integer, allocatable :: seeds(:) +#ifdef HAVE_MKL + type(ty_rng_mklvsl_plus) :: rng +#else + type(ty_rng_mt) :: rng +#endif + integer :: isub, icol, ilay, igpt + real :: cld_frac, sigma_qcw + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS + + ! allocate locals sized to this block + allocate(urand(ngpt,LM,ncols_block)) + allocate(cld_mask(ncols_block,LM,ngpt)) + if (gen_mro) then + allocate(urand_aux(ngpt,LM,ncols_block)) + allocate(alpha(ncols_block,LM-1)) + if (cond_inhomo) then + allocate(urand_cond (ngpt,LM,ncols_block)) + allocate(urand_cond_aux(ngpt,LM,ncols_block)) + allocate(rcorr(ncols_block,LM-1)) + allocate(zcw(ncols_block,LM,ngpt)) + end if + end if + allocate(seeds(3)) ! 2-word key plus word1 of counter + seeds = 0 + seeds(2) = seeds_time_key ! time part of key (computed once outside block loop) + ! for SW start at counter=65,536 + seeds(3) = 65536 - ! aerosols - TAUAER(:,1:LM,:) = TAUA(:,LM:1:-1,:) - SSAAER(:,1:LM,:) = SSAA(:,LM:1:-1,:) - ASMAER(:,1:LM,:) = ASYA(:,LM:1:-1,:) + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_CLOUD_OPTICS",__RC__) + + ! Make band in-cloud optical props from cloud_optics and mean in-cloud cloud water paths. + ! These can be scaled later to account for sub-gridscale condensate inhomogeneity. + ! Do phases separately to allow for different forward scattering, etc., per earlier note. + ! liquid ... + error_msg = cloud_optics%cloud_optics( & + real(QQ3(colS:colE,:,2),kind=wp) * dp_wp(colS:colE,:) * cwp_fac, & ! [g/m2] + dummy_wp(colS:colE,:), & + min( max( real(RR3(colS:colE,:,2),kind=wp), & ! [microns] + cloud_optics%get_min_radius_liq()), & + cloud_optics%get_max_radius_liq()), & + dummy_wp(colS:colE,:), & + cloud_props_bnd_liq) + TEST_(error_msg) + ! ice ... + error_msg = cloud_optics%cloud_optics( & + dummy_wp(colS:colE,:), & + real(QQ3(colS:colE,:,1),kind=wp) * dp_wp(colS:colE,:) * cwp_fac, & ! [g/m2] + dummy_wp(colS:colE,:), & + min( max( real(RR3(colS:colE,:,1),kind=wp), & ! [microns] + cloud_optics%get_min_radius_ice()), & + cloud_optics%get_max_radius_ice()), & + cloud_props_bnd_ice) + TEST_(error_msg) - call MAPL_TimerOff(MAPL,"--RRTMG_FLIP") - call MAPL_TimerOn (MAPL,"--RRTMG_INIT") + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_CLOUD_OPTICS",__RC__) - ! initialize RRTMG SW - call RRTMG_SW_INI + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_MCICA",__RC__) - call MAPL_TimerOff(MAPL,"--RRTMG_INIT") - call MAPL_TimerOn (MAPL,"--RRTMG_RUN") +!!TODO: need to resolve diff between prob of max vs ran and correlation coeff in both paper and code - ! partition size for columns (profiles) used to improve efficiency - call MAPL_GetResource( MAPL, RPART, 'RRTMGSW_PARTITION_SIZE:', DEFAULT=0, __RC__) + ! exponential inter-layer correlations + ! [alpha|rcorr](k) is correlation between layers k and k+1 + ! dzmid(k) is separation between midpoints of layers k and k+1 + if (gen_mro) then + do ilay = 1,LM-1 + ! cloud fraction correlation + alpha(:,ilay) = exp(-abs(dzmid(colS:colE,ilay))/real(adl(colS:colE),kind=wp)) + enddo + if (cond_inhomo) then + do ilay = 1,LM-1 + ! condensate correlation + rcorr(:,ilay) = exp(-abs(dzmid(colS:colE,ilay))/real(rdl(colS:colE),kind=wp)) + enddo + endif + endif - ! various RRTMG configuration options ... + ! generate McICA random numbers for block + do isub = 1, ncols_block + ! local 1d column index + icol = colS + isub - 1 + ! initialize the Philox PRNG + ! set word1 of key based on GLOBAL location + ! 32-bits can hold all forseeable resolutions + seeds(1) = nint(Jg1D(icol)) * IM_World + nint(Ig1D(icol)) +#ifdef HAVE_MKL + ! instantiate a random number stream for the column + call rng%init(VSL_BRNG_PHILOX4X32X10,seeds) +#else + call rng%init(seeds) +#endif + ! draw the random numbers for the column + urand(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) + if (gen_mro) then + urand_aux(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) + if (cond_inhomo) then + urand_cond (:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) + urand_cond_aux(:,:,isub) = reshape(rng%get_random(ngpt*LM),(/ngpt,LM/)) + endif + end if + ! free the rng + call rng%end() + end do - IAER = 10 ! Per AER: - ! 0: Turns off aerosols - ! 10: Enables aerosols + ! cloud sampling to gpoints + select case (cloud_overlap_type) + case ("MAX_RAN_OVERLAP") + error_msg = sampled_mask_max_ran( & + urand, real(CL(colS:colE,:),kind=wp), cld_mask) + TEST_(error_msg) + case ("EXP_RAN_OVERLAP") + ! corr_coeff(ncols_block,LM-1) is an inter-layer correlation coefficient + ! to be provided ... it is not the same as alpha, which is a probability +! error_msg = sampled_mask_exp_ran( & +! urand, real(CL(colS:colE,:),kind=wp), corr_coeff, cld_mask) +! TEST_(error_msg) + TEST_('EXP_RAN_OVERLAP not implemented yet') + case ("GEN_MAX_RAN_OVERLAP") + ! a scheme like Oreopoulos et al. 2012 (doi:10.5194/acp-12-9097-2012) in which both + ! cloud presence and cloud condensate are separately generalized maximum-random: + error_msg = sampled_urand_gen_max_ran(alpha, urand, urand_aux) + TEST_(error_msg) + if (cond_inhomo) then + error_msg = sampled_urand_gen_max_ran(rcorr, urand_cond, urand_cond_aux) + TEST_(error_msg) + end if + do isub = 1,ncols_block + icol = colS + isub - 1 + do ilay = 1,LM + cld_frac = CL(icol,ilay) - NORMFLX = 1 ! 0: Do not normalize fluxes - ! 1: Normalize fluxes + ! if grid-box clear, no subgrid variability + if (cld_frac <= 0.) then + cld_mask(isub,ilay,:) = .false. + else + ! subgrid-scale cloud mask + cld_mask(isub,ilay,:) = urand(:,ilay,isub) < cld_frac + + ! subgrid-scale condensate + if (cond_inhomo) then + ! level of condensate inhomogeneity based on cloud fraction. + if (cld_frac > 0.99) then + sigma_qcw = 0.5 + elseif (cld_frac > 0.9) then + sigma_qcw = 0.71 + else + sigma_qcw = 1.0 + endif + do igpt = 1,ngpt + if (cld_mask(isub,ilay,igpt)) zcw(isub,ilay,igpt) = & + zcw_lookup(real(urand_cond(igpt,ilay,isub)),sigma_qcw) + end do + end if + end if - DYOFYR = DOY ! Day of year + end do + end do - call MAPL_GetResource( MAPL, ISOLVAR ,'ISOLVAR:', DEFAULT=0, __RC__) + case default + TEST_('RRTMGP_SW: unknown cloud overlap') + end select - ! ISOLVAR: - ! Flag for solar variability method - ! -1 = (when scon .eq. 0.0): No solar variability - ! and no solar cycle (Kurucz solar irradiance - ! of 1368.22 Wm-2 only); - ! (when scon .ne. 0.0): Kurucz solar irradiance - ! scaled to scon and solar variability defined - ! (optional) by setting non-zero scale factors - ! for each band in bndsolvar - ! 0 = (when SCON .eq. 0.0): No solar variability - ! and no solar cycle (NRLSSI2 solar constant of - ! 1360.85 Wm-2 for the 100-50000 cm-1 spectral - ! range only), with facular and sunspot effects - ! fixed to the mean of Solar Cycles 13-24; - ! (when SCON .ne. 0.0): No solar variability - ! and no solar cycle (NRLSSI2 solar constant of - ! 1360.85 Wm-2 for the 100-50000 cm-1 spectral - ! range only), is scaled to SCON - ! 1 = Solar variability (using NRLSSI2 solar - ! model) with solar cycle contribution - ! determined by fraction of solar cycle - ! with facular and sunspot variations - ! fixed to their mean variations over the - ! average of Solar Cycles 13-24; - ! two amplitude scale factors allow - ! facular and sunspot adjustments from - ! mean solar cycle as defined by indsolvar - ! 2 = Solar variability (using NRLSSI2 solar - ! model) over solar cycle determined by - ! direct specification of Mg (facular) - ! and SB (sunspot) indices provided - ! in indsolvar (scon = 0.0 only) - ! 3 = (when scon .eq. 0.0): No solar variability - ! and no solar cycle (NRLSSI2 solar irradiance - ! of 1360.85 Wm-2 only); - ! (when scon .ne. 0.0): NRLSSI2 solar irradiance - ! scaled to scon and solar variability defined - ! (optional) by setting non-zero scale factors - ! for each band in bndsolvar + ! draw McICA optical property samples (band->gpt) + TEST_(draw_samples(cld_mask, cloud_props_bnd_liq, cloud_props_gpt_liq)) + TEST_(draw_samples(cld_mask, cloud_props_bnd_ice, cloud_props_gpt_ice)) - if (ISOLVAR == 1) then - if (MAPL_AM_I_ROOT()) then - write (*,*) "ERROR in RRTMG_SW" - write (*,*) "ISOLVAR==1 is currently unsupported as we have no" - write (*,*) "way of correctly setting solcycfrac." - end if - _FAIL('RRTMG SW: ISOLVAR==1 currently unsupported') + ! Scaling to sub-gridscale water paths: + ! since tau for each phase is linear in the phase's water path + ! and since the scaling zcw applies equally to both phases, the + ! total g-point optical thickness tau will scale with zcw. + if (gen_mro) then + if (cond_inhomo) & + where (cld_mask) cloud_props_gpt_liq%tau = cloud_props_gpt_liq%tau * zcw + where (cld_mask) cloud_props_gpt_ice%tau = cloud_props_gpt_ice%tau * zcw end if - ! INDSOLVAR = Facular and sunspot amplitude - ! scale factors (isolvar=1), or - ! Mg and SB indices (isolvar=2) - ! Dimensions: (2) + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_MCICA",__RC__) - if (ISOLVAR == 2) then + RETURN_(ESMF_SUCCESS) - ! Solar indices from our file - INDSOLVAR(1) = MG - INDSOLVAR(2) = SB + end subroutine compute_cloud_optics_mcica +#undef TEST_ - else + ! --------------------------------------------------------------------------- + ! Super-layer cloud fraction and in-cloud optical depth diagnostics, + ! computed from the McICA subcolumn ensemble BEFORE delta-scaling. + ! All scalar temporaries are local (thread-private under future OMP). + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_sprlyr_diags_predelta( & + colS, ncols_block, LM, ngpt, nbnd, LCLDLM, LCLDMH, & + include_aerosols, & + cld_mask, ClearCounts, CL, toa_flux, band_lims_gpt, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + CLDTS, CLDHS, CLDMS, CLDLS, & + COTTP, COTHP, COTMP, COTLP, & + COTDTP, COTDHP, COTDMP, COTDLP, & + COTNTP, COTNHP, COTNMP, COTNLP, & +#ifdef SOLAR_RADVAL + TAUTP, TAUHP, TAUMP, TAULP, & + COTLDTP, COTLDHP, COTLDMP, COTLDLP, & + COTLNTP, COTLNHP, COTLNMP, COTLNLP, & + COTIDTP, COTIDHP, COTIDMP, COTIDLP, & + COTINTP, COTINHP, COTINMP, COTINLP, & + SSALDTP, SSALDHP, SSALDMP, SSALDLP, & + SSALNTP, SSALNHP, SSALNMP, SSALNLP, & + SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & + SSAINTP, SSAINHP, SSAINMP, SSAINLP, & + ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & + ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & + ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & + ASMINTP, ASMINHP, ASMINMP, ASMINLP, & +#endif + MAPL, RC) + + use mo_optical_props, only: ty_optical_props_arry + + integer, intent(in) :: colS, ncols_block, LM, ngpt, nbnd + integer, intent(in) :: LCLDLM, LCLDMH + logical, intent(in) :: include_aerosols + logical, intent(in) :: cld_mask(:,:,:) + integer, intent(inout) :: ClearCounts(:,:) + real, intent(in) :: CL(:,:) + real(wp), intent(in) :: toa_flux(:,:) + integer, intent(in) :: band_lims_gpt(:,:) + class(ty_optical_props_arry), intent(in) :: cloud_props_gpt_liq, cloud_props_gpt_ice + real, intent(inout) :: CLDTS(:), CLDHS(:), CLDMS(:), CLDLS(:) + real, intent(inout) :: COTTP(:), COTHP(:), COTMP(:), COTLP(:) + real, intent(inout) :: COTDTP(:), COTDHP(:), COTDMP(:), COTDLP(:) + real, intent(inout) :: COTNTP(:), COTNHP(:), COTNMP(:), COTNLP(:) +#ifdef SOLAR_RADVAL + real, intent(inout) :: TAUTP(:), TAUHP(:), TAUMP(:), TAULP(:) + real, intent(inout) :: COTLDTP(:), COTLDHP(:), COTLDMP(:), COTLDLP(:) + real, intent(inout) :: COTLNTP(:), COTLNHP(:), COTLNMP(:), COTLNLP(:) + real, intent(inout) :: COTIDTP(:), COTIDHP(:), COTIDMP(:), COTIDLP(:) + real, intent(inout) :: COTINTP(:), COTINHP(:), COTINMP(:), COTINLP(:) + real, intent(inout) :: SSALDTP(:), SSALDHP(:), SSALDMP(:), SSALDLP(:) + real, intent(inout) :: SSALNTP(:), SSALNHP(:), SSALNMP(:), SSALNLP(:) + real, intent(inout) :: SSAIDTP(:), SSAIDHP(:), SSAIDMP(:), SSAIDLP(:) + real, intent(inout) :: SSAINTP(:), SSAINHP(:), SSAINMP(:), SSAINLP(:) + real, intent(inout) :: ASMLDTP(:), ASMLDHP(:), ASMLDMP(:), ASMLDLP(:) + real, intent(inout) :: ASMLNTP(:), ASMLNHP(:), ASMLNMP(:), ASMLNLP(:) + real, intent(inout) :: ASMIDTP(:), ASMIDHP(:), ASMIDMP(:), ASMIDLP(:) + real, intent(inout) :: ASMINTP(:), ASMINHP(:), ASMINMP(:), ASMINLP(:) +#endif + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! locals -- all thread-private under future !$OMP PARALLEL DO + integer :: isub, icol, ib, igpt + real :: wgt + real :: stautp, stauhp, staump, staulp + real :: sltautp, sltauhp, sltaump, sltaulp + real :: sitautp, sitauhp, sitaump, sitaulp +#ifdef SOLAR_RADVAL + real :: sltaussatp, sltaussahp, sltaussamp, sltaussalp + real :: sitaussatp, sitaussahp, sitaussamp, sitaussalp + real :: sltaussagtp, sltaussaghp, sltaussagmp, sltaussaglp + real :: sitaussagtp, sitaussaghp, sitaussagmp, sitaussaglp +#endif + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS - call MAPL_GetResource( MAPL, INDSOLVAR(1) ,'INDSOLVAR_1:', DEFAULT=1.0, __RC__) - call MAPL_GetResource( MAPL, INDSOLVAR(2) ,'INDSOLVAR_2:', DEFAULT=1.0, __RC__) + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - end if + if (include_aerosols) then + + ! super-layer cloud fractions + call clearCounts_threeBand( & + ncols_block, ncols_block, ngpt, LM, LCLDLM, LCLDMH, & + reshape(cld_mask,[LM,ngpt,ncols_block],order=[3,1,2]), & + ClearCounts) + do isub = 1,ncols_block + icol = colS + isub - 1 + CLDTS(icol) = 1. - ClearCounts(1,isub)/float(ngpt) + CLDHS(icol) = 1. - ClearCounts(2,isub)/float(ngpt) + CLDMS(icol) = 1. - ClearCounts(3,isub)/float(ngpt) + CLDLS(icol) = 1. - ClearCounts(4,isub)/float(ngpt) + end do + + ! in-cloud optical thicknesses in PAR super-band + ! (weighted across and within bands by TOA incident flux) + do isub = 1,ncols_block + icol = colS + isub - 1 + +#ifdef SOLAR_RADVAL + ! default (no cloud) for TAUx variant + TAUTP(icol) = 0. + TAUHP(icol) = 0. + TAUMP(icol) = 0. + TAULP(icol) = 0. +#endif + + ! default (no cloud) for COTx variant + COTTP(icol) = MAPL_UNDEF + COTHP(icol) = MAPL_UNDEF + COTMP(icol) = MAPL_UNDEF + COTLP(icol) = MAPL_UNDEF + + ! zero denom- and numerator accumulators + COTDTP(icol) = 0.; COTNTP(icol) = 0. + COTDHP(icol) = 0.; COTNHP(icol) = 0. + COTDMP(icol) = 0.; COTNMP(icol) = 0. + COTDLP(icol) = 0.; COTNLP(icol) = 0. +#ifdef SOLAR_RADVAL + COTLDTP(icol) = 0.; COTLNTP(icol) = 0.; COTIDTP(icol) = 0.; COTINTP(icol) = 0. + COTLDHP(icol) = 0.; COTLNHP(icol) = 0.; COTIDHP(icol) = 0.; COTINHP(icol) = 0. + COTLDMP(icol) = 0.; COTLNMP(icol) = 0.; COTIDMP(icol) = 0.; COTINMP(icol) = 0. + COTLDLP(icol) = 0.; COTLNLP(icol) = 0.; COTIDLP(icol) = 0.; COTINLP(icol) = 0. + SSALDTP(icol) = 0.; SSALNTP(icol) = 0.; SSAIDTP(icol) = 0.; SSAINTP(icol) = 0. + SSALDHP(icol) = 0.; SSALNHP(icol) = 0.; SSAIDHP(icol) = 0.; SSAINHP(icol) = 0. + SSALDMP(icol) = 0.; SSALNMP(icol) = 0.; SSAIDMP(icol) = 0.; SSAINMP(icol) = 0. + SSALDLP(icol) = 0.; SSALNLP(icol) = 0.; SSAIDLP(icol) = 0.; SSAINLP(icol) = 0. + ASMLDTP(icol) = 0.; ASMLNTP(icol) = 0.; ASMIDTP(icol) = 0.; ASMINTP(icol) = 0. + ASMLDHP(icol) = 0.; ASMLNHP(icol) = 0.; ASMIDHP(icol) = 0.; ASMINHP(icol) = 0. + ASMLDMP(icol) = 0.; ASMLNMP(icol) = 0.; ASMIDMP(icol) = 0.; ASMINMP(icol) = 0. + ASMLDLP(icol) = 0.; ASMLNLP(icol) = 0.; ASMIDLP(icol) = 0.; ASMINLP(icol) = 0. +#endif + ! can only be non-zero for potentially cloudy columns + if (any(CL(icol,:) > 0.)) then + + ! accumulate over gpts/subcolumns + do ib = 1, nbnd + do igpt = band_lims_gpt(1,ib), band_lims_gpt(2,ib) + + ! band weights for photosynthetically active radiation (PAR) + ! Bands 11-12 (0.345-0.625 um) plus half transition band 10 (0.625-0.778 um) + if (ib >= 11 .and. ib <= 12) then + wgt = 1.0 + else if (ib == 10) then + wgt = 0.5 + else + ! no contribution to PAR + cycle + end if + + ! TOA flux weighting + ! (note: neither the adjustment of toa_flux to our tsi + ! or for zenith angle are needed yet since this weighting + ! is over gpoint and is normalized for EACH icol) + wgt = wgt * toa_flux(isub,igpt) + + ! low pressure layer + sltaulp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt)) + sitaulp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt)) + staulp = sltaulp + sitaulp + if (staulp > 0.) then + COTDLP(icol) = COTDLP(icol) + wgt + COTNLP(icol) = COTNLP(icol) + wgt * staulp + end if +#ifdef SOLAR_RADVAL + sltaussalp = 0.; sltaussaglp = 0. + if (sltaulp > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussalp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt)) + sltaussaglp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%g (isub,LCLDLM:LM,igpt)) + end select + COTLDLP(icol) = COTLDLP(icol) + wgt + COTLNLP(icol) = COTLNLP(icol) + wgt * sltaulp + SSALDLP(icol) = SSALDLP(icol) + wgt * sltaulp + SSALNLP(icol) = SSALNLP(icol) + wgt * sltaussalp + ASMLDLP(icol) = ASMLDLP(icol) + wgt * sltaussalp + ASMLNLP(icol) = ASMLNLP(icol) + wgt * sltaussaglp + end if + sitaussalp = 0.; sitaussaglp = 0. + if (sitaulp > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussalp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt)) + sitaussaglp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%g (isub,LCLDLM:LM,igpt)) + end select + COTIDLP(icol) = COTIDLP(icol) + wgt + COTINLP(icol) = COTINLP(icol) + wgt * sitaulp + SSAIDLP(icol) = SSAIDLP(icol) + wgt * sitaulp + SSAINLP(icol) = SSAINLP(icol) + wgt * sitaussalp + ASMIDLP(icol) = ASMIDLP(icol) + wgt * sitaussalp + ASMINLP(icol) = ASMINLP(icol) + wgt * sitaussaglp + end if +#endif - BNDSOLVAR = 1.0 ! Solar variability scale factors - ! for each shortwave band - ! Dimensions: (nbndsw=14) + ! mid pressure layer + sltaump = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt)) + sitaump = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt)) + staump = sltaump + sitaump + if (staump > 0.) then + COTDMP(icol) = COTDMP(icol) + wgt + COTNMP(icol) = COTNMP(icol) + wgt * staump + end if +#ifdef SOLAR_RADVAL + sltaussamp = 0.; sltaussagmp = 0. + if (sltaump > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussamp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt)) + sltaussagmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%g (isub,LCLDMH:LCLDLM-1,igpt)) + end select + COTLDMP(icol) = COTLDMP(icol) + wgt + COTLNMP(icol) = COTLNMP(icol) + wgt * sltaump + SSALDMP(icol) = SSALDMP(icol) + wgt * sltaump + SSALNMP(icol) = SSALNMP(icol) + wgt * sltaussamp + ASMLDMP(icol) = ASMLDMP(icol) + wgt * sltaussamp + ASMLNMP(icol) = ASMLNMP(icol) + wgt * sltaussagmp + end if + sitaussamp = 0.; sitaussagmp = 0. + if (sitaump > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussamp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt)) + sitaussagmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%g (isub,LCLDMH:LCLDLM-1,igpt)) + end select + COTIDMP(icol) = COTIDMP(icol) + wgt + COTINMP(icol) = COTINMP(icol) + wgt * sitaump + SSAIDMP(icol) = SSAIDMP(icol) + wgt * sitaump + SSAINMP(icol) = SSAINMP(icol) + wgt * sitaussamp + ASMIDMP(icol) = ASMIDMP(icol) + wgt * sitaussamp + ASMINMP(icol) = ASMINMP(icol) + wgt * sitaussagmp + end if +#endif - ! SOLCYCFRAC: Fraction of averaged 11-year solar cycle (0-1) - ! at current time (isolvar=1) - ! 0.0 represents the first day of year 1 - ! 1.0 represents the last day of year 11 + ! high pressure layer + sltauhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt)) + sitauhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt)) + stauhp = sltauhp + sitauhp + if (stauhp > 0.) then + COTDHP(icol) = COTDHP(icol) + wgt + COTNHP(icol) = COTNHP(icol) + wgt * stauhp + end if +#ifdef SOLAR_RADVAL + sltaussahp = 0.; sltaussaghp = 0. + if (sltauhp > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussahp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt)) + sltaussaghp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%g (isub,1:LCLDMH-1,igpt)) + end select + COTLDHP(icol) = COTLDHP(icol) + wgt + COTLNHP(icol) = COTLNHP(icol) + wgt * sltauhp + SSALDHP(icol) = SSALDHP(icol) + wgt * sltauhp + SSALNHP(icol) = SSALNHP(icol) + wgt * sltaussahp + ASMLDHP(icol) = ASMLDHP(icol) + wgt * sltaussahp + ASMLNHP(icol) = ASMLNHP(icol) + wgt * sltaussaghp + end if + sitaussahp = 0.; sitaussaghp = 0. + if (sitauhp > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussahp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt)) + sitaussaghp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%g (isub,1:LCLDMH-1,igpt)) + end select + COTIDHP(icol) = COTIDHP(icol) + wgt + COTINHP(icol) = COTINHP(icol) + wgt * sitauhp + SSAIDHP(icol) = SSAIDHP(icol) + wgt * sitauhp + SSAINHP(icol) = SSAINHP(icol) + wgt * sitaussahp + ASMIDHP(icol) = ASMIDHP(icol) + wgt * sitaussahp + ASMINHP(icol) = ASMINHP(icol) + wgt * sitaussaghp + end if +#endif - ! MAT: Note while we don't currently use SOLCYCFRAC, we set it to something - ! to avoid an optional variable on GPUs + ! whole subcolumn + sltautp = sltaulp + sltaump + sltauhp + sitautp = sitaulp + sitaump + sitauhp + stautp = staulp + staump + stauhp + if (stautp > 0.) then + COTDTP(icol) = COTDTP(icol) + wgt + COTNTP(icol) = COTNTP(icol) + wgt * stautp + end if +#ifdef SOLAR_RADVAL + sltaussatp = sltaussalp + sltaussamp + sltaussahp + sltaussagtp = sltaussaglp + sltaussagmp + sltaussaghp + if (sltautp > 0.) then + COTLDTP(icol) = COTLDTP(icol) + wgt + COTLNTP(icol) = COTLNTP(icol) + wgt * sltautp + SSALDTP(icol) = SSALDTP(icol) + wgt * sltautp + SSALNTP(icol) = SSALNTP(icol) + wgt * sltaussatp + ASMLDTP(icol) = ASMLDTP(icol) + wgt * sltaussatp + ASMLNTP(icol) = ASMLNTP(icol) + wgt * sltaussagtp + end if + sitaussatp = sitaussalp + sitaussamp + sitaussahp + sitaussagtp = sitaussaglp + sitaussagmp + sitaussaghp + if (sitautp > 0.) then + COTIDTP(icol) = COTIDTP(icol) + wgt + COTINTP(icol) = COTINTP(icol) + wgt * sitautp + SSAIDTP(icol) = SSAIDTP(icol) + wgt * sitautp + SSAINTP(icol) = SSAINTP(icol) + wgt * sitaussatp + ASMIDTP(icol) = ASMIDTP(icol) + wgt * sitaussatp + ASMINTP(icol) = ASMINTP(icol) + wgt * sitaussagtp + end if +#endif - call MAPL_GetResource( MAPL, SOLCYCFRAC ,'SOLCYCFRAC:', DEFAULT=1.0, __RC__) + end do ! igpt + end do ! ib - ! call RRTMG SW - ! ------------- + ! normalize + ! Note: TAUx defaults zero, COTx defaults MAPL_UNDEF + if (COTDTP(icol) > 0. .and. COTNTP(icol) > 0.) then + COTTP(icol) = COTNTP(icol) / COTDTP(icol) +#ifdef SOLAR_RADVAL + TAUTP(icol) = COTTP(icol) +#endif + end if - call RRTMG_SW (MAPL, & - RPART, NCOL, LM, & - SC, ADJES, ZT, ISOLVAR, & - PL_R, PLE_R, T_R, & - Q_R, O3_R, CO2_R, CH4_R, O2_R, & - ICEFLGSW, LIQFLGSW, & - FCLD_R, CICEWP, CLIQWP, REICE, RELIQ, & - DYOFYR, ZL_R, ALAT, & - IAER, TAUAER, SSAAER, ASMAER, & - ALBVR, ALBVF, ALBNR, ALBNF, & - LM-LCLDLM+1, LM-LCLDMH+1, NORMFLX, & - CLEARCOUNTS, SWUFLX, SWDFLX, SWUFLXC, SWDFLXC, & - NIRR, NIRF, PARR, PARF, UVRR, UVRF, FSWBAND, & + if (COTDHP(icol) > 0. .and. COTNHP(icol) > 0.) then + COTHP(icol) = COTNHP(icol) / COTDHP(icol) +#ifdef SOLAR_RADVAL + TAUHP(icol) = COTHP(icol) +#endif + end if - COTDTP, COTDHP, COTDMP, COTDLP, & - COTNTP, COTNHP, COTNMP, COTNLP, & + if (COTDMP(icol) > 0. .and. COTNMP(icol) > 0.) then + COTMP(icol) = COTNMP(icol) / COTDMP(icol) +#ifdef SOLAR_RADVAL + TAUMP(icol) = COTMP(icol) +#endif + end if + if (COTDLP(icol) > 0. .and. COTNLP(icol) > 0.) then + COTLP(icol) = COTNLP(icol) / COTDLP(icol) #ifdef SOLAR_RADVAL - CDSDTP, CDSDHP, CDSDMP, CDSDLP, & - CDSNTP, CDSNHP, CDSNMP, CDSNLP, & + TAULP(icol) = COTLP(icol) +#endif + end if - COTLDTP, COTLDHP, COTLDMP, COTLDLP, & - COTLNTP, COTLNHP, COTLNMP, COTLNLP, & - CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & - CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & - COTIDTP, COTIDHP, COTIDMP, COTIDLP, & - COTINTP, COTINHP, COTINMP, COTINLP, & - CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & - CDSINTP, CDSINHP, CDSINMP, CDSINLP, & + end if ! potentially cloudy column + end do ! isub + end if ! include_aerosols - SSALDTP, SSALDHP, SSALDMP, SSALDLP, & - SSALNTP, SSALNHP, SSALNMP, SSALNLP, & - SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & - SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & - SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & - SSAINTP, SSAINHP, SSAINMP, SSAINLP, & - SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & - SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & - ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & - ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & - ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & - ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & - ASMINTP, ASMINHP, ASMINMP, ASMINLP, & - ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & - ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + RETURN_(ESMF_SUCCESS) - FORLDTP, FORLDHP, FORLDMP, FORLDLP, & - FORLNTP, FORLNHP, FORLNMP, FORLNLP, & - FORIDTP, FORIDHP, FORIDMP, FORIDLP, & - FORINTP, FORINHP, FORINMP, FORINLP, & -#endif + end subroutine compute_sprlyr_diags_predelta +#undef TEST_ - SOLAR_TO_OBIO .and. include_aerosols, DRBAND, DFBAND, & - BAND_OUTPUT .and. include_aerosols, ISRBRGN, OSRBRGN, & - BNDSOLVAR, INDSOLVAR, SOLCYCFRAC, & - __RC__) + ! --------------------------------------------------------------------------- + ! Delta-scale cloud optical properties (liquid and ice) for one block. + ! forwliq and forwice are computed here and returned for use in Step 2f. + ! All scalar temporaries are local (thread-private under future OMP). + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_delta_scale( & + colS, ncols_block, LM, ngpt, nbnd, & + rrtmgp_delta_scale, rrtmgp_use_rrtmg_iceflg3_like_forwice, & + CL, RR3, band_lims_gpt, & + cloud_optics, cloud_props_gpt_liq, cloud_props_gpt_ice, & + forwliq, forwice, & + MAPL, RC) + + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_2str + use mo_cloud_optics_rrtmgp, only: ty_cloud_optics_rrtmgp + use mo_rte_kind, only: wp + + integer, intent(in) :: colS, ncols_block, LM, ngpt, nbnd + logical, intent(in) :: rrtmgp_delta_scale + logical, intent(in) :: rrtmgp_use_rrtmg_iceflg3_like_forwice + real, dimension(:,:), intent(in) :: CL + real, dimension(:,:,:), intent(in) :: RR3 + integer, dimension(:,:), intent(in) :: band_lims_gpt + type(ty_cloud_optics_rrtmgp), intent(inout) :: cloud_optics + class(ty_optical_props_arry), intent(inout) :: cloud_props_gpt_liq, cloud_props_gpt_ice + real(wp), dimension(:,:,:), intent(out) :: forwliq, forwice + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! locals -- all thread-private under future !$OMP PARALLEL DO + integer :: isub, icol, ilay, ib, igpt, radidx + real(wp) :: radice_lwr, radice_upr, radice, radfac, rfint, fdelta + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS + + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_DELTA_SCALE",__RC__) + + forwliq = 0.; forwice = 0. ! default for no delta-scaling + if (rrtmgp_delta_scale) then + + ! default delta-scaling for liquid + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + forwliq = cloud_props_gpt_liq%g ** 2 + end select + TEST_(cloud_props_gpt_liq%delta_scale(forwliq)) - call MAPL_TimerOff(MAPL,"--RRTMG_RUN") - call MAPL_TimerOn (MAPL,"--RRTMG_FLIP") + if (rrtmgp_use_rrtmg_iceflg3_like_forwice) then + ! non-default delta-scaling for ice (as in RRTMG iceflag==3) + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + radice_lwr = cloud_optics%get_min_radius_ice() + radice_upr = cloud_optics%get_max_radius_ice() + do isub = 1,ncols_block + icol = colS + isub - 1 + do ilay = 1,LM + ! only if at least potentially cloudy ... + if (CL(icol,ilay) > 0.) then + + ! prepare for radice interpolation ... + ! first get radice consistent with RRTMGP ice cloud optics + radice = min(max(real(RR3(icol,ilay,1),kind=wp),radice_lwr),radice_upr) + ! now force into RRTMG's iceflag==3 reice binning range [5,140]um. + radice = min(max(radice,5._wp),140._wp) + ! RRTMG has 46 reice bins with 5um->radidx==1, 140um->radidx==46, + ! but radidx is forced to [1,45] so LIN2_ARG1 interpolation works. + radfac = (radice - 2._wp) / 3._wp + radidx = min(max(int(radfac),1),45) + rfint = radfac - real(radidx,kind=wp) + + do ib = 1,nbnd + ! interpolate fdelta in radice for band ib + fdelta = LIN2_ARG1(fdlice3_rrtmgp,radidx,ib,rfint) + + ! forwice calc for each g-point + do igpt = band_lims_gpt(1,ib),band_lims_gpt(2,ib) + if (cloud_props_gpt_ice%tau(isub,ilay,igpt) > 0.) then + forwice(isub,ilay,igpt) = min( & + fdelta + 0.5_wp / cloud_props_gpt_ice%ssa(isub,ilay,igpt), & + cloud_props_gpt_ice%g(isub,ilay,igpt)) + endif + enddo ! g-points + enddo ! bands + + endif ! potentially cloudy + enddo ! layers + enddo ! columns + end select + TEST_(cloud_props_gpt_ice%delta_scale(forwice)) + else + ! default delta-scaling for ice + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + forwice = cloud_props_gpt_ice%g ** 2 + end select + TEST_(cloud_props_gpt_ice%delta_scale(forwice)) + endif + endif - ! unflip the outputs in the vertical - ! ---------------------------------- + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_DELTA_SCALE",__RC__) - SWUFLXR (:,1:LM+1) = SWUFLX (:,LM+1:1:-1) - SWDFLXR (:,1:LM+1) = SWDFLX (:,LM+1:1:-1) - SWUFLXCR(:,1:LM+1) = SWUFLXC(:,LM+1:1:-1) - SWDFLXCR(:,1:LM+1) = SWDFLXC(:,LM+1:1:-1) + RETURN_(ESMF_SUCCESS) - call MAPL_TimerOff(MAPL,"--RRTMG_FLIP") + end subroutine compute_delta_scale +#undef TEST_ - ! required outputs - ! ---------------- + ! --------------------------------------------------------------------------- + ! Super-layer cloud optical depth diagnostics AFTER delta-scaling. + ! Only compiled/called under #ifdef SOLAR_RADVAL. + ! forwliq and forwice (from compute_delta_scale) are intent(in) here. + ! All scalar temporaries are local (thread-private under future OMP). + ! --------------------------------------------------------------------------- +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_sprlyr_diags_postdelta( & + colS, ncols_block, LM, ngpt, nbnd, LCLDLM, LCLDMH, & + include_aerosols, & + CL, toa_flux, band_lims_gpt, & + forwliq, forwice, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + CDSDTP, CDSDHP, CDSDMP, CDSDLP, & + CDSNTP, CDSNHP, CDSNMP, CDSNLP, & + CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & + CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & + CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & + CDSINTP, CDSINHP, CDSINMP, CDSINLP, & + SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & + SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & + SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & + SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & + ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & + ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & + ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + FORLDTP, FORLDHP, FORLDMP, FORLDLP, & + FORLNTP, FORLNHP, FORLNMP, FORLNLP, & + FORIDTP, FORIDHP, FORIDMP, FORIDLP, & + FORINTP, FORINHP, FORINMP, FORINLP, & + MAPL, RC) + + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_2str + use mo_rte_kind, only: wp + + integer, intent(in) :: colS, ncols_block, LM, ngpt, nbnd + integer, intent(in) :: LCLDLM, LCLDMH + logical, intent(in) :: include_aerosols + real, intent(in) :: CL(:,:) + real(wp), intent(in) :: toa_flux(:,:) + integer, intent(in) :: band_lims_gpt(:,:) + real(wp), intent(in) :: forwliq(:,:,:), forwice(:,:,:) + class(ty_optical_props_arry), intent(in) :: cloud_props_gpt_liq, cloud_props_gpt_ice + real, intent(inout) :: CDSDTP(:), CDSDHP(:), CDSDMP(:), CDSDLP(:) + real, intent(inout) :: CDSNTP(:), CDSNHP(:), CDSNMP(:), CDSNLP(:) + real, intent(inout) :: CDSLDTP(:), CDSLDHP(:), CDSLDMP(:), CDSLDLP(:) + real, intent(inout) :: CDSLNTP(:), CDSLNHP(:), CDSLNMP(:), CDSLNLP(:) + real, intent(inout) :: CDSIDTP(:), CDSIDHP(:), CDSIDMP(:), CDSIDLP(:) + real, intent(inout) :: CDSINTP(:), CDSINHP(:), CDSINMP(:), CDSINLP(:) + real, intent(inout) :: SDSLDTP(:), SDSLDHP(:), SDSLDMP(:), SDSLDLP(:) + real, intent(inout) :: SDSLNTP(:), SDSLNHP(:), SDSLNMP(:), SDSLNLP(:) + real, intent(inout) :: SDSIDTP(:), SDSIDHP(:), SDSIDMP(:), SDSIDLP(:) + real, intent(inout) :: SDSINTP(:), SDSINHP(:), SDSINMP(:), SDSINLP(:) + real, intent(inout) :: ADSLDTP(:), ADSLDHP(:), ADSLDMP(:), ADSLDLP(:) + real, intent(inout) :: ADSLNTP(:), ADSLNHP(:), ADSLNMP(:), ADSLNLP(:) + real, intent(inout) :: ADSIDTP(:), ADSIDHP(:), ADSIDMP(:), ADSIDLP(:) + real, intent(inout) :: ADSINTP(:), ADSINHP(:), ADSINMP(:), ADSINLP(:) + real, intent(inout) :: FORLDTP(:), FORLDHP(:), FORLDMP(:), FORLDLP(:) + real, intent(inout) :: FORLNTP(:), FORLNHP(:), FORLNMP(:), FORLNLP(:) + real, intent(inout) :: FORIDTP(:), FORIDHP(:), FORIDMP(:), FORIDLP(:) + real, intent(inout) :: FORINTP(:), FORINHP(:), FORINMP(:), FORINLP(:) + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! locals -- all thread-private under future !$OMP PARALLEL DO + integer :: isub, icol, ib, igpt + real :: wgt + real :: stautp, stauhp, staump, staulp + real :: sltautp, sltauhp, sltaump, sltaulp + real :: sitautp, sitauhp, sitaump, sitaulp + real :: sltaussatp, sltaussahp, sltaussamp, sltaussalp + real :: sitaussatp, sitaussahp, sitaussamp, sitaussalp + real :: sltaussagtp, sltaussaghp, sltaussagmp, sltaussaglp + real :: sitaussagtp, sitaussaghp, sitaussagmp, sitaussaglp + real :: sltaussaftp, sltaussafhp, sltaussafmp, sltaussaflp + real :: sitaussaftp, sitaussafhp, sitaussafmp, sitaussaflp + character(len=ESMF_MAXSTR) :: error_msg + integer :: STATUS + + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - ! convert super-layer clearCounts to cloud fractions if (include_aerosols) then - CLDTS(:) = 1. - CLEARCOUNTS(:,1)/float(NGPTSW) - CLDHS(:) = 1. - CLEARCOUNTS(:,2)/float(NGPTSW) - CLDMS(:) = 1. - CLEARCOUNTS(:,3)/float(NGPTSW) - CLDLS(:) = 1. - CLEARCOUNTS(:,4)/float(NGPTSW) - end if - ! undef versions of cloud optical thicknesses - ! We cannot use a merge here because some compilers (e.g. Intel) - ! will will evaluate the first argument first and if both are - ! zero, it will return NaNs. - where (COTNTP > 0. .and. COTDTP > 0.) - COTTP = COTNTP/COTDTP - elsewhere - COTTP = MAPL_UNDEF - end where + ! in-cloud optical thicknesses in PAR super-band + ! (weighted across and within bands by TOA incident flux) + do isub = 1,ncols_block + icol = colS + isub - 1 - where (COTNHP > 0. .and. COTDHP > 0.) - COTHP = COTNHP/COTDHP - elsewhere - COTHP = MAPL_UNDEF - end where + ! zero denom- and numerator accumulators + CDSDTP(icol) = 0.; CDSNTP(icol) = 0. + CDSDHP(icol) = 0.; CDSNHP(icol) = 0. + CDSDMP(icol) = 0.; CDSNMP(icol) = 0. + CDSDLP(icol) = 0.; CDSNLP(icol) = 0. + + CDSLDTP(icol) = 0.; CDSLNTP(icol) = 0.; CDSIDTP(icol) = 0.; CDSINTP(icol) = 0. + CDSLDHP(icol) = 0.; CDSLNHP(icol) = 0.; CDSIDHP(icol) = 0.; CDSINHP(icol) = 0. + CDSLDMP(icol) = 0.; CDSLNMP(icol) = 0.; CDSIDMP(icol) = 0.; CDSINMP(icol) = 0. + CDSLDLP(icol) = 0.; CDSLNLP(icol) = 0.; CDSIDLP(icol) = 0.; CDSINLP(icol) = 0. + + SDSLDTP(icol) = 0.; SDSLNTP(icol) = 0.; SDSIDTP(icol) = 0.; SDSINTP(icol) = 0. + SDSLDHP(icol) = 0.; SDSLNHP(icol) = 0.; SDSIDHP(icol) = 0.; SDSINHP(icol) = 0. + SDSLDMP(icol) = 0.; SDSLNMP(icol) = 0.; SDSIDMP(icol) = 0.; SDSINMP(icol) = 0. + SDSLDLP(icol) = 0.; SDSLNLP(icol) = 0.; SDSIDLP(icol) = 0.; SDSINLP(icol) = 0. + + ADSLDTP(icol) = 0.; ADSLNTP(icol) = 0.; ADSIDTP(icol) = 0.; ADSINTP(icol) = 0. + ADSLDHP(icol) = 0.; ADSLNHP(icol) = 0.; ADSIDHP(icol) = 0.; ADSINHP(icol) = 0. + ADSLDMP(icol) = 0.; ADSLNMP(icol) = 0.; ADSIDMP(icol) = 0.; ADSINMP(icol) = 0. + ADSLDLP(icol) = 0.; ADSLNLP(icol) = 0.; ADSIDLP(icol) = 0.; ADSINLP(icol) = 0. + + FORLDTP(icol) = 0.; FORLNTP(icol) = 0.; FORIDTP(icol) = 0.; FORINTP(icol) = 0. + FORLDHP(icol) = 0.; FORLNHP(icol) = 0.; FORIDHP(icol) = 0.; FORINHP(icol) = 0. + FORLDMP(icol) = 0.; FORLNMP(icol) = 0.; FORIDMP(icol) = 0.; FORINMP(icol) = 0. + FORLDLP(icol) = 0.; FORLNLP(icol) = 0.; FORIDLP(icol) = 0.; FORINLP(icol) = 0. + + ! can only be non-zero for potentially cloudy columns + if (any(CL(icol,:) > 0.)) then + + ! accumulate over gpts/subcolumns + do ib = 1, nbnd + do igpt = band_lims_gpt(1,ib), band_lims_gpt(2,ib) + + ! band weights for photosynthetically active radiation (PAR) + ! Bands 11-12 (0.345-0.625 um) plus half transition band 10 (0.625-0.778 um) + if (ib >= 11 .and. ib <= 12) then + wgt = 1.0 + else if (ib == 10) then + wgt = 0.5 + else + ! no contribution to PAR + cycle + end if - where (COTNMP > 0. .and. COTDMP > 0.) - COTMP = COTNMP/COTDMP - elsewhere - COTMP = MAPL_UNDEF - end where + wgt = wgt * toa_flux(isub,igpt) - where (COTNLP > 0. .and. COTDLP > 0.) - COTLP = COTNLP/COTDLP - elsewhere - COTLP = MAPL_UNDEF - end where + ! low pressure layer + sltaulp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt)) + sitaulp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt)) + staulp = sltaulp + sitaulp + if (staulp > 0.) then + CDSNLP(icol) = CDSNLP(icol) + wgt * staulp + CDSDLP(icol) = CDSDLP(icol) + wgt + end if + sltaussalp = 0.; sltaussaglp = 0.; sltaussaflp = 0. + if (sltaulp > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussalp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt)) + sltaussaglp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%g (isub,LCLDLM:LM,igpt)) + sltaussaflp = sum(cloud_props_gpt_liq%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDLM:LM,igpt) * & + forwliq(isub,LCLDLM:LM,igpt)) + end select + CDSLDLP(icol) = CDSLDLP(icol) + wgt + CDSLNLP(icol) = CDSLNLP(icol) + wgt * sltaulp + SDSLDLP(icol) = SDSLDLP(icol) + wgt * sltaulp + SDSLNLP(icol) = SDSLNLP(icol) + wgt * sltaussalp + ADSLDLP(icol) = ADSLDLP(icol) + wgt * sltaussalp + ADSLNLP(icol) = ADSLNLP(icol) + wgt * sltaussaglp + FORLDLP(icol) = FORLDLP(icol) + wgt * sltaussalp + FORLNLP(icol) = FORLNLP(icol) + wgt * sltaussaflp + end if + sitaussalp = 0.; sitaussaglp = 0.; sitaussaflp = 0. + if (sitaulp > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussalp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt)) + sitaussaglp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%g (isub,LCLDLM:LM,igpt)) + sitaussaflp = sum(cloud_props_gpt_ice%tau(isub,LCLDLM:LM,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDLM:LM,igpt) * & + forwice(isub,LCLDLM:LM,igpt)) + end select + CDSIDLP(icol) = CDSIDLP(icol) + wgt + CDSINLP(icol) = CDSINLP(icol) + wgt * sitaulp + SDSIDLP(icol) = SDSIDLP(icol) + wgt * sitaulp + SDSINLP(icol) = SDSINLP(icol) + wgt * sitaussalp + ADSIDLP(icol) = ADSIDLP(icol) + wgt * sitaussalp + ADSINLP(icol) = ADSINLP(icol) + wgt * sitaussaglp + FORIDLP(icol) = FORIDLP(icol) + wgt * sitaussalp + FORINLP(icol) = FORINLP(icol) + wgt * sitaussaflp + end if -#ifdef SOLAR_RADVAL - ! zero versions of cloud optical thicknesses - ! We can use merge() here because we cannot divide by zero - TAUTP = merge(COTTP, 0., COTNTP > 0. .and. COTDTP > 0.) - TAUHP = merge(COTHP, 0., COTNHP > 0. .and. COTDHP > 0.) - TAUMP = merge(COTMP, 0., COTNMP > 0. .and. COTDMP > 0.) - TAULP = merge(COTLP, 0., COTNLP > 0. .and. COTDLP > 0.) -#endif + ! mid pressure layer + sltaump = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt)) + sitaump = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt)) + staump = sltaump + sitaump + if (staump > 0.) then + CDSNMP(icol) = CDSNMP(icol) + wgt * staump + CDSDMP(icol) = CDSDMP(icol) + wgt + end if + sltaussamp = 0.; sltaussagmp = 0.; sltaussafmp = 0. + if (sltaump > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussamp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt)) + sltaussagmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%g (isub,LCLDMH:LCLDLM-1,igpt)) + sltaussafmp = sum(cloud_props_gpt_liq%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + forwliq(isub,LCLDMH:LCLDLM-1,igpt)) + end select + CDSLDMP(icol) = CDSLDMP(icol) + wgt + CDSLNMP(icol) = CDSLNMP(icol) + wgt * sltaump + SDSLDMP(icol) = SDSLDMP(icol) + wgt * sltaump + SDSLNMP(icol) = SDSLNMP(icol) + wgt * sltaussamp + ADSLDMP(icol) = ADSLDMP(icol) + wgt * sltaussamp + ADSLNMP(icol) = ADSLNMP(icol) + wgt * sltaussagmp + FORLDMP(icol) = FORLDMP(icol) + wgt * sltaussamp + FORLNMP(icol) = FORLNMP(icol) + wgt * sltaussafmp + end if + sitaussamp = 0.; sitaussagmp = 0.; sitaussafmp = 0. + if (sitaump > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussamp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt)) + sitaussagmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%g (isub,LCLDMH:LCLDLM-1,igpt)) + sitaussafmp = sum(cloud_props_gpt_ice%tau(isub,LCLDMH:LCLDLM-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,LCLDMH:LCLDLM-1,igpt) * & + forwice(isub,LCLDMH:LCLDLM-1,igpt)) + end select + CDSIDMP(icol) = CDSIDMP(icol) + wgt + CDSINMP(icol) = CDSINMP(icol) + wgt * sitaump + SDSIDMP(icol) = SDSIDMP(icol) + wgt * sitaump + SDSINMP(icol) = SDSINMP(icol) + wgt * sitaussamp + ADSIDMP(icol) = ADSIDMP(icol) + wgt * sitaussamp + ADSINMP(icol) = ADSINMP(icol) + wgt * sitaussagmp + FORIDMP(icol) = FORIDMP(icol) + wgt * sitaussamp + FORINMP(icol) = FORINMP(icol) + wgt * sitaussafmp + end if - ! fluxes - FSW = SWDFLXR - SWUFLXR - FSC = SWDFLXCR - SWUFLXCR - FSWU = SWUFLXR - FSCU = SWUFLXCR + ! high pressure layer + sltauhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt)) + sitauhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt)) + stauhp = sltauhp + sitauhp + if (stauhp > 0.) then + CDSNHP(icol) = CDSNHP(icol) + wgt * stauhp + CDSDHP(icol) = CDSDHP(icol) + wgt + end if + sltaussahp = 0.; sltaussaghp = 0.; sltaussafhp = 0. + if (sltauhp > 0.) then + select type(cloud_props_gpt_liq) + class is (ty_optical_props_2str) + sltaussahp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt)) + sltaussaghp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%g (isub,1:LCLDMH-1,igpt)) + sltaussafhp = sum(cloud_props_gpt_liq%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_liq%ssa(isub,1:LCLDMH-1,igpt) * & + forwliq(isub,1:LCLDMH-1,igpt)) + end select + CDSLDHP(icol) = CDSLDHP(icol) + wgt + CDSLNHP(icol) = CDSLNHP(icol) + wgt * sltauhp + SDSLDHP(icol) = SDSLDHP(icol) + wgt * sltauhp + SDSLNHP(icol) = SDSLNHP(icol) + wgt * sltaussahp + ADSLDHP(icol) = ADSLDHP(icol) + wgt * sltaussahp + ADSLNHP(icol) = ADSLNHP(icol) + wgt * sltaussaghp + FORLDHP(icol) = FORLDHP(icol) + wgt * sltaussahp + FORLNHP(icol) = FORLNHP(icol) + wgt * sltaussafhp + end if + sitaussahp = 0.; sitaussaghp = 0.; sitaussafhp = 0. + if (sitauhp > 0.) then + select type(cloud_props_gpt_ice) + class is (ty_optical_props_2str) + sitaussahp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt)) + sitaussaghp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%g (isub,1:LCLDMH-1,igpt)) + sitaussafhp = sum(cloud_props_gpt_ice%tau(isub,1:LCLDMH-1,igpt) * & + cloud_props_gpt_ice%ssa(isub,1:LCLDMH-1,igpt) * & + forwice(isub,1:LCLDMH-1,igpt)) + end select + CDSIDHP(icol) = CDSIDHP(icol) + wgt + CDSINHP(icol) = CDSINHP(icol) + wgt * sitauhp + SDSIDHP(icol) = SDSIDHP(icol) + wgt * sitauhp + SDSINHP(icol) = SDSINHP(icol) + wgt * sitaussahp + ADSIDHP(icol) = ADSIDHP(icol) + wgt * sitaussahp + ADSINHP(icol) = ADSINHP(icol) + wgt * sitaussaghp + FORIDHP(icol) = FORIDHP(icol) + wgt * sitaussahp + FORINHP(icol) = FORINHP(icol) + wgt * sitaussafhp + end if - ! Deallocate the working inputs - !------------------------------ - deallocate(TLEV ,__STAT__) - deallocate(TLEV_R,__STAT__) - deallocate(PLE_R ,__STAT__) - deallocate(FCLD_R,__STAT__) - deallocate(CLIQWP,__STAT__) - deallocate(CICEWP,__STAT__) - deallocate(RELIQ ,__STAT__) - deallocate(REICE ,__STAT__) + ! whole subcolumn + sltautp = sltaulp + sltaump + sltauhp + sitautp = sitaulp + sitaump + sitauhp + stautp = staulp + staump + stauhp + if (stautp > 0.) then + CDSNTP(icol) = CDSNTP(icol) + wgt * stautp + CDSDTP(icol) = CDSDTP(icol) + wgt + end if + sltaussatp = sltaussalp + sltaussamp + sltaussahp + sltaussagtp = sltaussaglp + sltaussagmp + sltaussaghp + sltaussaftp = sltaussaflp + sltaussafmp + sltaussafhp + if (sltautp > 0.) then + CDSLDTP(icol) = CDSLDTP(icol) + wgt + CDSLNTP(icol) = CDSLNTP(icol) + wgt * sltautp + SDSLDTP(icol) = SDSLDTP(icol) + wgt * sltautp + SDSLNTP(icol) = SDSLNTP(icol) + wgt * sltaussatp + ADSLDTP(icol) = ADSLDTP(icol) + wgt * sltaussatp + ADSLNTP(icol) = ADSLNTP(icol) + wgt * sltaussagtp + FORLDTP(icol) = FORLDTP(icol) + wgt * sltaussatp + FORLNTP(icol) = FORLNTP(icol) + wgt * sltaussaftp + end if + sitaussatp = sitaussalp + sitaussamp + sitaussahp + sitaussagtp = sitaussaglp + sitaussagmp + sitaussaghp + sitaussaftp = sitaussaflp + sitaussafmp + sitaussafhp + if (sitautp > 0.) then + CDSIDTP(icol) = CDSIDTP(icol) + wgt + CDSINTP(icol) = CDSINTP(icol) + wgt * sitautp + SDSIDTP(icol) = SDSIDTP(icol) + wgt * sitautp + SDSINTP(icol) = SDSINTP(icol) + wgt * sitaussatp + ADSIDTP(icol) = ADSIDTP(icol) + wgt * sitaussatp + ADSINTP(icol) = ADSINTP(icol) + wgt * sitaussagtp + FORIDTP(icol) = FORIDTP(icol) + wgt * sitaussatp + FORINTP(icol) = FORINTP(icol) + wgt * sitaussaftp + end if - deallocate(TAUAER,__STAT__) - deallocate(SSAAER,__STAT__) - deallocate(ASMAER,__STAT__) - deallocate(DPR ,__STAT__) - deallocate(PL_R ,__STAT__) - deallocate(ZL_R ,__STAT__) - deallocate(T_R ,__STAT__) - deallocate(Q_R ,__STAT__) - deallocate(O2_R ,__STAT__) - deallocate(O3_R ,__STAT__) - deallocate(CO2_R ,__STAT__) - deallocate(CH4_R ,__STAT__) + end do ! igpt + end do ! ib - deallocate(CLEARCOUNTS ,__STAT__) + end if ! potentially cloudy column + end do ! isub + end if ! include_aerosols - deallocate(SWUFLX ,__STAT__) - deallocate(SWDFLX ,__STAT__) - deallocate(SWUFLXC,__STAT__) - deallocate(SWDFLXC,__STAT__) + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_SPRLYR_DIAGS",__RC__) - deallocate(SWUFLXR ,__STAT__) - deallocate(SWDFLXR ,__STAT__) - deallocate(SWUFLXCR,__STAT__) - deallocate(SWDFLXCR,__STAT__) + RETURN_(ESMF_SUCCESS) - call MAPL_TimerOff(MAPL,"-RRTMG") + end subroutine compute_sprlyr_diags_postdelta +#undef TEST_ - else +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine compute_rte_sw( & + colS, colE, ngpt, & + tsi, toa_flux, optical_props, & + top_at_1, mu0, sfc_alb_dir, sfc_alb_dif, & + fluxes_clrsky, flux_up_clrsky, flux_net_clrsky, & + fluxes_allsky, flux_up_allsky, flux_net_allsky, & + bnd_flux_dn_allsky, bnd_flux_dir_allsky, bnd_flux_net_allsky, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + MAPL, RC) + + use mo_optical_props, only: ty_optical_props_arry + use mo_rte_kind, only: wp + use mo_rte_sw, only: rte_sw + use mo_fluxes_byband, only: ty_fluxes_byband + + integer, intent(in) :: colS, colE, ngpt + real(wp), intent(in) :: tsi(:) + real(wp), intent(inout) :: toa_flux(:,:) + class(ty_optical_props_arry), intent(inout) :: optical_props + logical, intent(in) :: top_at_1 + real(wp), intent(in) :: mu0(:) + real(wp), intent(in) :: sfc_alb_dir(:,:) + real(wp), intent(in) :: sfc_alb_dif(:,:) + type(ty_fluxes_byband), intent(inout) :: fluxes_clrsky + real(wp), target, intent(inout) :: flux_up_clrsky(:,:) + real(wp), target, intent(inout) :: flux_net_clrsky(:,:) + type(ty_fluxes_byband), intent(inout) :: fluxes_allsky + real(wp), target, intent(inout) :: flux_up_allsky(:,:) + real(wp), target, intent(inout) :: flux_net_allsky(:,:) + real(wp), target, intent(inout) :: bnd_flux_dn_allsky(:,:,:) + real(wp), target, intent(inout) :: bnd_flux_dir_allsky(:,:,:) + real(wp), target, intent(inout) :: bnd_flux_net_allsky(:,:,:) + class(ty_optical_props_arry), intent(inout) :: cloud_props_gpt_liq + class(ty_optical_props_arry), intent(inout) :: cloud_props_gpt_ice + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + character(len=512) :: error_msg + integer :: STATUS + + ! MAPL_TimerOn disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOn(MAPL,"--RRTMGP_RT",__RC__) + + ! scale to our tsi + ! (both toa_flux and tsi are NORMAL to solar beam, [W/m2]) + toa_flux = toa_flux * spread(tsi(colS:colE)/sum(toa_flux,dim=2), 2, ngpt) + + ! clear-sky radiative transfer + fluxes_clrsky%flux_up => flux_up_clrsky + fluxes_clrsky%flux_net => flux_net_clrsky + error_msg = rte_sw( & + optical_props, top_at_1, mu0(colS:colE), toa_flux, & + sfc_alb_dir(:,colS:colE), sfc_alb_dif(:,colS:colE), & + fluxes_clrsky) + TEST_(error_msg) - _FAIL('unknown SW radiation scheme!') + ! add in cloud optical properties + ! add ice first since its optical depths are usually smaller + TEST_(cloud_props_gpt_ice%increment(optical_props)) + TEST_(cloud_props_gpt_liq%increment(optical_props)) + + ! all-sky radiative transfer + fluxes_allsky%flux_up => flux_up_allsky + fluxes_allsky%flux_net => flux_net_allsky + fluxes_allsky%bnd_flux_dn => bnd_flux_dn_allsky + fluxes_allsky%bnd_flux_dn_dir => bnd_flux_dir_allsky + fluxes_allsky%bnd_flux_net => bnd_flux_net_allsky + error_msg = rte_sw( & + optical_props, top_at_1, mu0(colS:colE), toa_flux, & + sfc_alb_dir(:,colS:colE), sfc_alb_dif(:,colS:colE), & + fluxes_allsky) + TEST_(error_msg) - end if SCHEME + ! MAPL_TimerOff disabled inside OMP parallel region (not thread-safe) + ! call MAPL_TimerOff(MAPL,"--RRTMGP_RT",__RC__) - ! Deallocate the working inputs - !------------------------------ + RETURN + end subroutine compute_rte_sw +#undef TEST_ - deallocate (PL, RH, PLhPa) - deallocate (QQ3, RR3) - deallocate (ILWT) - deallocate (O3) - deallocate (TAUA, SSAA, ASYA) +#define TEST_(A) error_msg = A; if (trim(error_msg)/="") then; _FAIL("RRTMGP Error: "//trim(error_msg)); endif + subroutine PROCESS_RRTMGP_BLOCK( & + b, ncol, rrtmgp_blockSize, LM, ngpt, nbnd, nmom, & + LCLDLM, LCLDMH, include_aerosols, gen_mro, cond_inhomo, & + cloud_overlap_type, IM_World, seeds_time_key, & + need_aer_optical_props, top_at_1, & + rrtmgp_delta_scale, rrtmgp_use_rrtmg_iceflg3_like_forwice, & + cwp_fac_arg, & + gas_concs, k_dist, cloud_optics, & + p_lay, p_lev, t_lay, QQ3, RR3, dp_wp, dummy_wp, & + CL, dzmid, adl, rdl, Ig1D, Jg1D, band_lims_gpt, & + tsi, mu0, sfc_alb_dir, sfc_alb_dif, taua, ssaa, asya, & + flux_up_clrsky, flux_net_clrsky, & + flux_up_allsky, flux_net_allsky, & + bnd_flux_dn_allsky, bnd_flux_dir_allsky, bnd_flux_net_allsky, & + CLDTS, CLDHS, CLDMS, CLDLS, & + COTTP, COTHP, COTMP, COTLP, & + COTDTP, COTDHP, COTDMP, COTDLP, & + COTNTP, COTNHP, COTNMP, COTNLP, & +#ifdef SOLAR_RADVAL + TAUTP, TAUHP, TAUMP, TAULP, & + COTLDTP, COTLDHP, COTLDMP, COTLDLP, & + COTLNTP, COTLNHP, COTLNMP, COTLNLP, & + COTIDTP, COTIDHP, COTIDMP, COTIDLP, & + COTINTP, COTINHP, COTINMP, COTINLP, & + SSALDTP, SSALDHP, SSALDMP, SSALDLP, & + SSALNTP, SSALNHP, SSALNMP, SSALNLP, & + SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & + SSAINTP, SSAINHP, SSAINMP, SSAINLP, & + ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & + ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & + ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & + ASMINTP, ASMINHP, ASMINMP, ASMINLP, & + CDSDTP, CDSDHP, CDSDMP, CDSDLP, & + CDSNTP, CDSNHP, CDSNMP, CDSNLP, & + CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & + CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & + CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & + CDSINTP, CDSINHP, CDSINMP, CDSINLP, & + SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & + SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & + SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & + SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & + ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & + ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & + ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + FORLDTP, FORLDHP, FORLDMP, FORLDLP, & + FORLNTP, FORLNHP, FORLNMP, FORLNLP, & + FORIDTP, FORIDHP, FORIDMP, FORIDLP, & + FORINTP, FORINHP, FORINMP, FORINLP, & +#endif + MAPL, RC) + + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_1scl, & + ty_optical_props_2str, ty_optical_props_nstr + use mo_rte_kind, only: wp + use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp + use mo_gas_concentrations, only: ty_gas_concs + use mo_cloud_optics_rrtmgp, only: ty_cloud_optics_rrtmgp + use mo_fluxes_byband, only: ty_fluxes_byband + + integer, intent(in) :: b, ncol, rrtmgp_blockSize + integer, intent(in) :: LM, ngpt, nbnd, nmom + integer, intent(in) :: LCLDLM, LCLDMH + logical, intent(in) :: include_aerosols + logical, intent(in) :: gen_mro, cond_inhomo + character(len=*), intent(in) :: cloud_overlap_type + integer, intent(in) :: IM_World + integer, intent(in) :: seeds_time_key + logical, intent(in) :: need_aer_optical_props + logical, intent(in) :: top_at_1 + logical, intent(in) :: rrtmgp_delta_scale + logical, intent(in) :: rrtmgp_use_rrtmg_iceflg3_like_forwice + real(wp), intent(in) :: cwp_fac_arg + type(ty_gas_concs), intent(in) :: gas_concs + type(ty_gas_optics_rrtmgp), intent(in) :: k_dist + type(ty_cloud_optics_rrtmgp), intent(inout) :: cloud_optics + real(wp), intent(in) :: p_lay(:,:), p_lev(:,:), t_lay(:,:) + real, intent(in) :: QQ3(:,:,:), RR3(:,:,:) + real(wp), intent(in) :: dp_wp(:,:), dummy_wp(:,:) + real, intent(in) :: CL(:,:) + real(wp), intent(in) :: dzmid(:,:) + real, intent(in) :: adl(:), rdl(:) + real, intent(in) :: Ig1D(:), Jg1D(:) + integer, intent(in) :: band_lims_gpt(:,:) + real(wp), intent(in) :: tsi(:), mu0(:) + real(wp), intent(in) :: sfc_alb_dir(:,:), sfc_alb_dif(:,:) + real, intent(in) :: taua(:,:,:), ssaa(:,:,:), asya(:,:,:) + real(wp), target, intent(inout) :: flux_up_clrsky(:,:), flux_net_clrsky(:,:) + real(wp), target, intent(inout) :: flux_up_allsky(:,:), flux_net_allsky(:,:) + real(wp), target, intent(inout) :: bnd_flux_dn_allsky(:,:,:) + real(wp), target, intent(inout) :: bnd_flux_dir_allsky(:,:,:) + real(wp), target, intent(inout) :: bnd_flux_net_allsky(:,:,:) + real, intent(inout) :: CLDTS(:), CLDHS(:), CLDMS(:), CLDLS(:) + real, intent(inout) :: COTTP(:), COTHP(:), COTMP(:), COTLP(:) + real, intent(inout) :: COTDTP(:), COTDHP(:), COTDMP(:), COTDLP(:) + real, intent(inout) :: COTNTP(:), COTNHP(:), COTNMP(:), COTNLP(:) +#ifdef SOLAR_RADVAL + real, intent(inout) :: TAUTP(:), TAUHP(:), TAUMP(:), TAULP(:) + real, intent(inout) :: COTLDTP(:), COTLDHP(:), COTLDMP(:), COTLDLP(:) + real, intent(inout) :: COTLNTP(:), COTLNHP(:), COTLNMP(:), COTLNLP(:) + real, intent(inout) :: COTIDTP(:), COTIDHP(:), COTIDMP(:), COTIDLP(:) + real, intent(inout) :: COTINTP(:), COTINHP(:), COTINMP(:), COTINLP(:) + real, intent(inout) :: SSALDTP(:), SSALDHP(:), SSALDMP(:), SSALDLP(:) + real, intent(inout) :: SSALNTP(:), SSALNHP(:), SSALNMP(:), SSALNLP(:) + real, intent(inout) :: SSAIDTP(:), SSAIDHP(:), SSAIDMP(:), SSAIDLP(:) + real, intent(inout) :: SSAINTP(:), SSAINHP(:), SSAINMP(:), SSAINLP(:) + real, intent(inout) :: ASMLDTP(:), ASMLDHP(:), ASMLDMP(:), ASMLDLP(:) + real, intent(inout) :: ASMLNTP(:), ASMLNHP(:), ASMLNMP(:), ASMLNLP(:) + real, intent(inout) :: ASMIDTP(:), ASMIDHP(:), ASMIDMP(:), ASMIDLP(:) + real, intent(inout) :: ASMINTP(:), ASMINHP(:), ASMINMP(:), ASMINLP(:) + real, intent(inout) :: CDSDTP(:), CDSDHP(:), CDSDMP(:), CDSDLP(:) + real, intent(inout) :: CDSNTP(:), CDSNHP(:), CDSNMP(:), CDSNLP(:) + real, intent(inout) :: CDSLDTP(:), CDSLDHP(:), CDSLDMP(:), CDSLDLP(:) + real, intent(inout) :: CDSLNTP(:), CDSLNHP(:), CDSLNMP(:), CDSLNLP(:) + real, intent(inout) :: CDSIDTP(:), CDSIDHP(:), CDSIDMP(:), CDSIDLP(:) + real, intent(inout) :: CDSINTP(:), CDSINHP(:), CDSINMP(:), CDSINLP(:) + real, intent(inout) :: SDSLDTP(:), SDSLDHP(:), SDSLDMP(:), SDSLDLP(:) + real, intent(inout) :: SDSLNTP(:), SDSLNHP(:), SDSLNMP(:), SDSLNLP(:) + real, intent(inout) :: SDSIDTP(:), SDSIDHP(:), SDSIDMP(:), SDSIDLP(:) + real, intent(inout) :: SDSINTP(:), SDSINHP(:), SDSINMP(:), SDSINLP(:) + real, intent(inout) :: ADSLDTP(:), ADSLDHP(:), ADSLDMP(:), ADSLDLP(:) + real, intent(inout) :: ADSLNTP(:), ADSLNHP(:), ADSLNMP(:), ADSLNLP(:) + real, intent(inout) :: ADSIDTP(:), ADSIDHP(:), ADSIDMP(:), ADSIDLP(:) + real, intent(inout) :: ADSINTP(:), ADSINHP(:), ADSINMP(:), ADSINLP(:) + real, intent(inout) :: FORLDTP(:), FORLDHP(:), FORLDMP(:), FORLDLP(:) + real, intent(inout) :: FORLNTP(:), FORLNHP(:), FORLNMP(:), FORLNLP(:) + real, intent(inout) :: FORIDTP(:), FORIDHP(:), FORIDMP(:), FORIDLP(:) + real, intent(inout) :: FORINTP(:), FORINHP(:), FORINMP(:), FORINLP(:) +#endif + type(MAPL_MetaComp), intent(inout) :: MAPL + integer, optional, intent(out) :: RC + + ! per-block private local variables + integer :: colS, colE, ncols_block + integer, allocatable :: ClearCounts(:,:) + logical, allocatable :: cld_mask(:,:,:) + real(wp), allocatable :: toa_flux(:,:) + real(wp), allocatable :: forwliq(:,:,:), forwice(:,:,:) + class(ty_optical_props_arry), allocatable :: optical_props + class(ty_optical_props_arry), allocatable :: cloud_props_bnd_liq, cloud_props_bnd_ice + class(ty_optical_props_arry), allocatable :: cloud_props_gpt_liq, cloud_props_gpt_ice + class(ty_optical_props_arry), allocatable :: aer_props + type(ty_fluxes_byband) :: fluxes_clrsky, fluxes_allsky + + character(len=512) :: error_msg + integer :: STATUS + + ! compute block column range; the final block may be partial + ncols_block = min(rrtmgp_blockSize, ncol - (b-1)*rrtmgp_blockSize) + colS = (b-1) * rrtmgp_blockSize + 1 + colE = colS + ncols_block - 1 + + ! allocate per-block arrays + allocate(toa_flux(ncols_block,ngpt), __STAT__) + allocate(forwliq(ncols_block,LM,ngpt), __STAT__) + allocate(forwice(ncols_block,LM,ngpt), __STAT__) + if (include_aerosols) & + allocate(ClearCounts(4,ncols_block), __STAT__) - ! Complete load balancing by retrieving work done remotely - !--------------------------------------------------------- + ! instantiate optical_props with desired streams + allocate(ty_optical_props_2str::optical_props,__STAT__) ! <-- choose 2-stream SW + TEST_(optical_props%init(k_dist)) - call MAPL_TimerOn(MAPL,"-BALANCE") + ! cloud optics: band-space and g-point allocatables (thread-private) + allocate(ty_optical_props_2str::cloud_props_bnd_liq,__STAT__) + allocate(ty_optical_props_2str::cloud_props_bnd_ice,__STAT__) + TEST_(cloud_props_bnd_liq%init(k_dist%get_band_lims_wavenumber())) + TEST_(cloud_props_bnd_ice%init(k_dist%get_band_lims_wavenumber())) - call MAPL_TimerOn(MAPL,"--RETRIEVE") - if (LoadBalance) then - if (size(BufOut) > 0) call MAPL_BalanceWork(BufOut, NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) - if (size(BufInOut) > 0) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + select type (cloud_props_bnd_liq) + class is (ty_optical_props_2str) + allocate(ty_optical_props_2str::cloud_props_gpt_liq,__STAT__) + class default + TEST_('cloud optical properties (liq) hardwired 2-stream for now') + end select + select type (cloud_props_bnd_ice) + class is (ty_optical_props_2str) + allocate(ty_optical_props_2str::cloud_props_gpt_ice,__STAT__) + class default + TEST_('cloud optical properties (ice) hardwired 2-stream for now') + end select + TEST_(cloud_props_gpt_liq%init(k_dist)) + TEST_(cloud_props_gpt_ice%init(k_dist)) + + ! aerosol optical properties (thread-private) + if (need_aer_optical_props) then + allocate(ty_optical_props_2str::aer_props,__STAT__) + TEST_(aer_props%init(k_dist%get_band_lims_wavenumber())) end if - call MAPL_TimerOff(MAPL,"--RETRIEVE") - ! Unpack the results. Fills "masked" (night) locations with default value from internal state - !-------------------------------------------------------------------------------------------- - ! resulting internals are then contiguous versions - ! Note: InOut variables do not fill unmasked locations with a default, - ! since the unmasked locations may contain potentially useful aged data. + ! ty_optical_props routines have an internal deallocation + select type (cloud_props_bnd_liq) + class is (ty_optical_props_2str) + TEST_(cloud_props_bnd_liq%alloc_2str(ncols_block,LM)) + end select + select type (cloud_props_bnd_ice) + class is (ty_optical_props_2str) + TEST_(cloud_props_bnd_ice%alloc_2str(ncols_block,LM)) + end select + select type (cloud_props_gpt_liq) + class is (ty_optical_props_2str) + TEST_(cloud_props_gpt_liq%alloc_2str(ncols_block,LM)) + end select + select type (cloud_props_gpt_ice) + class is (ty_optical_props_2str) + TEST_(cloud_props_gpt_ice%alloc_2str(ncols_block,LM)) + end select + if (need_aer_optical_props) then + select type (aer_props) + class is (ty_optical_props_2str) + TEST_(aer_props%alloc_2str(ncols_block,LM)) + end select + end if + select type (optical_props) + class is (ty_optical_props_1scl) + TEST_(optical_props%alloc_1scl(ncols_block,LM)) + class is (ty_optical_props_2str) + TEST_(optical_props%alloc_2str(ncols_block,LM)) + class is (ty_optical_props_nstr) + TEST_(optical_props%alloc_nstr(nmom,ncols_block,LM)) + end select - i1InOut = 1; i1Out = 1 - INT_VARS_3: do k=1,NumInt - if (SlicesInt(k) == 0) cycle + call compute_gas_optics(colS, colE, ncols_block, LM, & + gas_concs, k_dist, p_lay, p_lev, t_lay, & + optical_props, toa_flux, MAPL, __RC__) - if (IntInOut(k)) then - pi1 => i1InOut - else - pi1 => i1Out - call MAPL_VarSpecGet(InternalSpec(k),DEFAULT=def,__RC__) - endif + if (need_aer_optical_props) then + call compute_aer_optics(colS, colE, need_aer_optical_props, & + taua, ssaa, asya, aer_props, __RC__) + end if - if (ugDim(k) > 0) then - select case(rgDim(k)) - case(MAPL_DIMSHORZVERT) - call ESMFL_StateGetPointerToData(INTERNAL,ptr4,NamesInt(k),__RC__) - if (IntInOut(k)) then - do j=1,ugDim(k) - call UnPackIt(BufInOut(pi1+(j-1)*size(ptr4,3)*NumMax),ptr4(:,:,:,j), & - daytime,NumMax,HorzDims,size(ptr4,3)) - end do - else - do j=1,ugDim(k) - call UnPackIt(BufOut (pi1+(j-1)*size(ptr4,3)*NumMax),ptr4(:,:,:,j), & - daytime,NumMax,HorzDims,size(ptr4,3),def) - end do - endif - case(MAPL_DIMSHORZONLY) - call ESMFL_StateGetPointerToData(INTERNAL,ptr3,NamesInt(k),__RC__) - if (IntInOut(k)) then - call UnPackIt(BufInOut(pi1),ptr3,daytime,NumMax,HorzDims,ugDim(k)) - else - call UnPackIt(BufOut (pi1),ptr3,daytime,NumMax,HorzDims,ugDim(k),def) - endif - end select - else - select case(rgDim(k)) - case(MAPL_DIMSHORZVERT) - call ESMFL_StateGetPointerToData(INTERNAL,ptr3,NamesInt(k),__RC__) - if (IntInOut(k)) then - call UnPackIt(BufInOut(pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3)) - else - call UnPackIt(BufOut (pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3),def) - endif - case(MAPL_DIMSHORZONLY) - call ESMFL_StateGetPointerToData(INTERNAL,ptr2,NamesInt(k),__RC__) - if (IntInOut(k)) then - call UnPackIt(BufInOut(pi1),ptr2,daytime,NumMax,HorzDims,1) - else - call UnPackIt(BufOut (pi1),ptr2,daytime,NumMax,HorzDims,1,def) - endif - end select - end if - pi1 = pi1 + NumMax*SlicesInt(k) + call compute_cloud_optics_mcica( & + colS, colE, ncols_block, LM, ngpt, & + gen_mro, cond_inhomo, cloud_overlap_type, cwp_fac_arg, IM_World, & + seeds_time_key, & + QQ3, RR3, dp_wp, dummy_wp, CL, dzmid, adl, rdl, Ig1D, Jg1D, & + cloud_optics, & + cloud_props_bnd_liq, cloud_props_bnd_ice, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + cld_mask, & + MAPL, __RC__) + + ! REFRESH super-layer diagnostics (before delta-scaling TAUs). + ! ** Calculated from subcolumn ensemble, so stochastic ** + ! ------------------------------------------------------- + call compute_sprlyr_diags_predelta( & + colS, ncols_block, LM, ngpt, nbnd, LCLDLM, LCLDMH, & + include_aerosols, & + cld_mask, ClearCounts, CL, toa_flux, band_lims_gpt, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + CLDTS, CLDHS, CLDMS, CLDLS, & + COTTP, COTHP, COTMP, COTLP, & + COTDTP, COTDHP, COTDMP, COTDLP, & + COTNTP, COTNHP, COTNMP, COTNLP, & +#ifdef SOLAR_RADVAL + TAUTP, TAUHP, TAUMP, TAULP, & + COTLDTP, COTLDHP, COTLDMP, COTLDLP, & + COTLNTP, COTLNHP, COTLNMP, COTLNLP, & + COTIDTP, COTIDHP, COTIDMP, COTIDLP, & + COTINTP, COTINHP, COTINMP, COTINLP, & + SSALDTP, SSALDHP, SSALDMP, SSALDLP, & + SSALNTP, SSALNHP, SSALNMP, SSALNLP, & + SSAIDTP, SSAIDHP, SSAIDMP, SSAIDLP, & + SSAINTP, SSAINHP, SSAINMP, SSAINLP, & + ASMLDTP, ASMLDHP, ASMLDMP, ASMLDLP, & + ASMLNTP, ASMLNHP, ASMLNMP, ASMLNLP, & + ASMIDTP, ASMIDHP, ASMIDMP, ASMIDLP, & + ASMINTP, ASMINHP, ASMINMP, ASMINLP, & +#endif + MAPL, __RC__) - enddo INT_VARS_3 + ! delta-scaling of cloud optical properties (accounts for forward scattering) + call compute_delta_scale( & + colS, ncols_block, LM, ngpt, nbnd, & + rrtmgp_delta_scale, rrtmgp_use_rrtmg_iceflg3_like_forwice, & + CL, RR3, band_lims_gpt, & + cloud_optics, cloud_props_gpt_liq, cloud_props_gpt_ice, & + forwliq, forwice, & + MAPL, __RC__) - ! clean up - deallocate(SlicesInp,NamesInp,__STAT__) - deallocate(SlicesInt,NamesInt,__STAT__) - deallocate(IntInOut,rgDim,ugDim,__STAT__) - deallocate(BufInp,BufInOut,BufOut,__STAT__) - call MAPL_TimerOn(MAPL,"--DESTROY") - if (LoadBalance) call MAPL_BalanceDestroy(Handle=SolarBalanceHandle, __RC__) - call MAPL_TimerOff(MAPL,"--DESTROY") +#ifdef SOLAR_RADVAL + ! REFRESH super-layer diagnostics (after delta-scaling TAUs). + ! ** Calculated from subcolumn ensemble, so stochastic ** + ! ------------------------------------------------------- + call compute_sprlyr_diags_postdelta( & + colS, ncols_block, LM, ngpt, nbnd, LCLDLM, LCLDMH, & + include_aerosols, & + CL, toa_flux, band_lims_gpt, & + forwliq, forwice, & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + CDSDTP, CDSDHP, CDSDMP, CDSDLP, & + CDSNTP, CDSNHP, CDSNMP, CDSNLP, & + CDSLDTP, CDSLDHP, CDSLDMP, CDSLDLP, & + CDSLNTP, CDSLNHP, CDSLNMP, CDSLNLP, & + CDSIDTP, CDSIDHP, CDSIDMP, CDSIDLP, & + CDSINTP, CDSINHP, CDSINMP, CDSINLP, & + SDSLDTP, SDSLDHP, SDSLDMP, SDSLDLP, & + SDSLNTP, SDSLNHP, SDSLNMP, SDSLNLP, & + SDSIDTP, SDSIDHP, SDSIDMP, SDSIDLP, & + SDSINTP, SDSINHP, SDSINMP, SDSINLP, & + ADSLDTP, ADSLDHP, ADSLDMP, ADSLDLP, & + ADSLNTP, ADSLNHP, ADSLNMP, ADSLNLP, & + ADSIDTP, ADSIDHP, ADSIDMP, ADSIDLP, & + ADSINTP, ADSINHP, ADSINMP, ADSINLP, & + FORLDTP, FORLDHP, FORLDMP, FORLDLP, & + FORLNTP, FORLNHP, FORLNMP, FORLNLP, & + FORIDTP, FORIDHP, FORIDMP, FORIDLP, & + FORINTP, FORINHP, FORINMP, FORINLP, & + MAPL, __RC__) +#endif - call MAPL_TimerOff(MAPL,"-BALANCE") + ! add in aerosol optical properties if requested and available + ! (done here rather than inside compute_rte_sw to avoid passing + ! an unallocated polymorphic aer_props when aerosols are inactive) + if (need_aer_optical_props) then + TEST_(aer_props%increment(optical_props)) + end if - RETURN_(ESMF_SUCCESS) - end subroutine SORADCORE + call compute_rte_sw( & + colS, colE, ngpt, & + tsi, toa_flux, optical_props, & + top_at_1, mu0, sfc_alb_dir, sfc_alb_dif, & + fluxes_clrsky, flux_up_clrsky(colS:colE,:), flux_net_clrsky(colS:colE,:), & + fluxes_allsky, flux_up_allsky(colS:colE,:), flux_net_allsky(colS:colE,:), & + bnd_flux_dn_allsky(colS:colE,:,:), bnd_flux_dir_allsky(colS:colE,:,:), & + bnd_flux_net_allsky(colS:colE,:,:), & + cloud_props_gpt_liq, cloud_props_gpt_ice, & + MAPL, __RC__) + + ! deallocate per-block arrays + deallocate(toa_flux, __STAT__) + deallocate(cld_mask, __STAT__) + deallocate(forwliq, __STAT__) + deallocate(forwice, __STAT__) + if (include_aerosols) & + deallocate(ClearCounts, __STAT__) + ! cloud_props_*, aer_props, optical_props are local allocatables; + ! they are automatically finalized/deallocated on return. + + end subroutine PROCESS_RRTMGP_BLOCK +#undef TEST_ subroutine SHRTWAVE(PLhPa,TA,WA,OA,CO2,COSZ , &