diff --git a/.gitmodules b/.gitmodules index 8c794079dc..3492317cc6 100644 --- a/.gitmodules +++ b/.gitmodules @@ -35,8 +35,8 @@ [submodule "atmos_phys"] path = src/atmos_phys - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_11_000 + url = https://github.com/peverwhee/atmospheric_physics + fxtag = 49e6ec240f53dad382602d4b325d9198d8b399fc fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics diff --git a/bld/configure b/bld/configure index 20f5142d10..a7df51e824 100755 --- a/bld/configure +++ b/bld/configure @@ -2144,6 +2144,8 @@ sub write_filepath } print $fh "$camsrcdir/src/physics/rrtmgp/ext/rrtmgp-kernels\n"; print $fh "$camsrcdir/src/physics/rrtmgp/ext/rte-kernels\n"; + print $fh "$camsrcdir/src/atmos_phys/schemes/rrtmgp\n"; + print $fh "$camsrcdir/src/atmos_phys/schemes/rrtmgp/objects\n"; } if ($clubb_sgs) { @@ -2175,6 +2177,7 @@ sub write_filepath print $fh "$camsrcdir/src/atmos_phys/schemes/conservation_adjust/check_energy\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/conservation_adjust/dme_adjust\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/hack_shallow\n"; + print $fh "$camsrcdir/src/atmos_phys/schemes/rrtmgp/utils\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/utilities\n"; print $fh "$camsrcdir/src/atmos_phys/schemes/cloud_fraction\n"; diff --git a/src/physics/cam/cloud_rad_props.F90 b/src/physics/cam/cloud_rad_props.F90 index 257138e7b5..894f4a1356 100644 --- a/src/physics/cam/cloud_rad_props.F90 +++ b/src/physics/cam/cloud_rad_props.F90 @@ -37,6 +37,7 @@ module cloud_rad_props get_snow_optics_sw, & snow_cloud_get_rad_props_lw, & get_grau_optics_sw, & + get_mu_lambda_weights, & grau_cloud_get_rad_props_lw @@ -77,8 +78,9 @@ module cloud_rad_props contains !============================================================================== -subroutine cloud_rad_props_init() - +subroutine cloud_rad_props_init(nmu_out, nlambda_out, n_g_d_out, & + abs_lw_liq_out, abs_lw_ice_out, g_mu_out, g_lambda_out, & + g_d_eff_out, tiny_out) use netcdf use spmd_utils, only: masterproc use ioFileMod, only: getfil @@ -86,6 +88,15 @@ subroutine cloud_rad_props_init() #if ( defined SPMD ) use mpishorthand #endif + integer, intent(out) :: nmu_out + integer, intent(out) :: nlambda_out + integer, intent(out) :: n_g_d_out + real(r8), allocatable, intent(out) :: abs_lw_liq_out(:,:,:) + real(r8), allocatable, intent(out) :: abs_lw_ice_out(:,:) + real(r8), allocatable, intent(out) :: g_mu_out(:) + real(r8), allocatable, intent(out) :: g_lambda_out(:,:) + real(r8), allocatable, intent(out) :: g_d_eff_out(:) + real(r8), intent(out) :: tiny_out character(len=256) :: liquidfile character(len=256) :: icefile @@ -103,6 +114,7 @@ subroutine cloud_rad_props_init() integer :: err character(len=*), parameter :: sub = 'cloud_rad_props_init' + character(len=512) :: errmsg liquidfile = liqopticsfile icefile = iceopticsfile @@ -278,6 +290,36 @@ subroutine cloud_rad_props_init() call mpibcast(abs_lw_ice, n_g_d*nlwbands, mpir8, 0, mpicom, ierr) #endif + ! Set output variables + tiny_out = tiny + nmu_out = nmu + nlambda_out = nlambda + n_g_d_out = n_g_d + allocate(abs_lw_liq_out(nmu,nlambda,nlwbands), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + call endrun(sub//': Failed to allocate abs_lw_liq_out - message: '//errmsg) + end if + abs_lw_liq_out = abs_lw_liq + allocate(abs_lw_ice_out(n_g_d,nlwbands), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + call endrun(sub//': Failed to allocate abs_lw_ice_out - message: '//errmsg) + end if + abs_lw_ice_out = abs_lw_ice + allocate(g_mu_out(nmu), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + call endrun(sub//': Failed to allocate g_mu_out - message: '//errmsg) + end if + g_mu_out = g_mu + allocate(g_lambda_out(nmu,nlambda), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + call endrun(sub//': Failed to allocate g_lambda_out - message: '//errmsg) + end if + g_lambda_out = g_lambda + allocate(g_d_eff_out(n_g_d), stat=ierr, errmsg=errmsg) + if (ierr /= 0) then + call endrun(sub//': Failed to allocate g_d_eff_out - message: '//errmsg) + end if + g_d_eff_out = g_d_eff return end subroutine cloud_rad_props_init @@ -728,28 +770,21 @@ end subroutine gam_liquid_sw !============================================================================== subroutine get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts) + use radiation_utils, only: get_mu_lambda_weights_ccpp real(r8), intent(in) :: lamc ! prognosed value of lambda for cloud real(r8), intent(in) :: pgam ! prognosed value of mu for cloud ! Output interpolation weights. Caller is responsible for freeing these. type(interp_type), intent(out) :: mu_wgts type(interp_type), intent(out) :: lambda_wgts - integer :: ilambda - real(r8) :: g_lambda_interp(nlambda) - - ! Make interpolation weights for mu. - ! (Put pgam in a temporary array for this purpose.) - call lininterp_init(g_mu, nmu, [pgam], 1, extrap_method_bndry, mu_wgts) - - ! Use mu weights to interpolate to a row in the lambda table. - do ilambda = 1, nlambda - call lininterp(g_lambda(:,ilambda), nmu, & - g_lambda_interp(ilambda:ilambda), 1, mu_wgts) - end do + character(len=512) :: errmsg + integer :: errflg - ! Make interpolation weights for lambda. - call lininterp_init(g_lambda_interp, nlambda, [lamc], 1, & - extrap_method_bndry, lambda_wgts) + call get_mu_lambda_weights_ccpp(nmu, nlambda, g_mu, g_lambda, lamc, pgam, mu_wgts, & + lambda_wgts, errmsg, errflg) + if (errflg /= 0) then + call endrun('get_mu_lambda_weights: ERROR message: '//errmsg) + end if end subroutine get_mu_lambda_weights diff --git a/src/physics/cam/radheat.F90 b/src/physics/cam/radheat.F90 index 37f8127931..15ab38a843 100644 --- a/src/physics/cam/radheat.F90 +++ b/src/physics/cam/radheat.F90 @@ -82,6 +82,7 @@ subroutine radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & #if ( defined OFFLINE_DYN ) use metdata, only: met_rlx, met_srf_feedback #endif + use calculate_net_heating, only: calculate_net_heating_run !----------------------------------------------------------------------- ! Compute net radiative heating from qrs and qrl, and the associated net ! boundary flux. @@ -91,7 +92,7 @@ subroutine radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & type(physics_state), intent(in) :: state ! Physics state variables type(physics_buffer_desc), pointer :: pbuf(:) - type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencie + type(physics_ptend), intent(out) :: ptend ! individual parameterization tendencies real(r8), intent(in) :: qrl(pcols,pver) ! longwave heating real(r8), intent(in) :: qrs(pcols,pver) ! shortwave heating real(r8), intent(in) :: fsns(pcols) ! Surface solar absorbed flux @@ -105,12 +106,23 @@ subroutine radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & ! Local variables integer :: i, k integer :: ncol + character(len=512) :: errmsg + integer :: errflg !----------------------------------------------------------------------- + ! Set error variables + errmsg = '' + errflg = 0 + ncol = state%ncol call physics_ptend_init(ptend,state%psetcols, 'radheat', ls=.true.) + ! REMOVECAM no longer need once CAM is retired and pcols doesn't exist + ptend%s(:,:) = 0._r8 + net_flx(:) = 0._r8 + ! END_REMOVECAM + #if ( defined OFFLINE_DYN ) ptend%s(:ncol,:) = 0._r8 do k = 1,pver @@ -118,14 +130,13 @@ subroutine radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & ptend%s(:ncol,k) = (qrs(:ncol,k) + qrl(:ncol,k)) endif enddo + call calculate_net_heating_run(ncol, ptend%s(:ncol,:), qrl(:ncol,:), qrs(:ncol,:), fsns, fsnt, flns, & + flnt, .true., net_flx(:ncol), errmsg, errflg) #else - ptend%s(:ncol,:) = (qrs(:ncol,:) + qrl(:ncol,:)) + call calculate_net_heating_run(ncol, ptend%s(:ncol,:), qrl(:ncol,:), qrs(:ncol,:), fsns, fsnt, flns, & + flnt, .false., net_flx(:ncol), errmsg, errflg) #endif - do i = 1, ncol - net_flx(i) = fsnt(i) - fsns(i) - flnt(i) + flns(i) - end do - end subroutine radheat_tend !================================================================================================ diff --git a/src/physics/rrtmg/radiation.F90 b/src/physics/rrtmg/radiation.F90 index a4c0cae8f8..ff79c7dd71 100644 --- a/src/physics/rrtmg/radiation.F90 +++ b/src/physics/rrtmg/radiation.F90 @@ -382,6 +382,15 @@ subroutine radiation_init(pbuf2d) integer :: history_budget_histfile_num ! output history file number for budget fields integer :: err + ! Cloud optics variables + integer :: nmu, n_g_d, nlambda + real(kind=r8), allocatable :: abs_lw_ice(:,:) + real(kind=r8), allocatable :: abs_lw_liq(:,:,:) + real(kind=r8), allocatable :: g_lambda(:,:) + real(kind=r8), allocatable :: g_mu(:) + real(kind=r8), allocatable :: g_d_eff(:) + real(kind=r8) :: tiny + integer :: dtime !----------------------------------------------------------------------- @@ -390,7 +399,8 @@ subroutine radiation_init(pbuf2d) call rad_data_init(pbuf2d) ! initialize output fields for offline driver call radsw_init() call radlw_init() - call cloud_rad_props_init() + call cloud_rad_props_init(nmu, nlambda, n_g_d, abs_lw_liq, abs_lw_ice, & + g_mu, g_lambda, g_d_eff, tiny) cld_idx = pbuf_get_index('CLD') cldfsnow_idx = pbuf_get_index('CLDFSNOW',errcode=err) diff --git a/src/physics/rrtmgp/mcica_subcol_gen.F90 b/src/physics/rrtmgp/mcica_subcol_gen.F90 index 4c592646e2..8af483a917 100644 --- a/src/physics/rrtmgp/mcica_subcol_gen.F90 +++ b/src/physics/rrtmgp/mcica_subcol_gen.F90 @@ -36,126 +36,12 @@ module mcica_subcol_gen private save -public :: mcica_subcol_lw, mcica_subcol_sw +public :: mcica_subcol_sw !======================================================================================== contains !======================================================================================== -subroutine mcica_subcol_lw( & - kdist, nbnd, ngpt, ncol, nver, & - changeseed, pmid, cldfrac, tauc, taucmcl ) - - ! Arrays use CAM vertical index convention: index increases from top to bottom. - ! This index ordering is assumed in the maximum-random overlap algorithm which starts - ! at the top of a column and marches down, with each layer depending on the state - ! of the layer above it. - ! - ! For GCM mode, changeseed must be offset between LW and SW by at least the - ! number of subcolumns - - ! arguments - class(ty_gas_optics_rrtmgp), intent(in) :: kdist ! spectral information - integer, intent(in) :: nbnd ! number of spectral bands - integer, intent(in) :: ngpt ! number of subcolumns (g-point intervals) - integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: nver ! number of layers - integer, intent(in) :: changeseed ! if the subcolumn generator is called multiple times, - ! permute the seed between each call. - real(r8), intent(in) :: pmid(pcols,pver) ! layer pressures (Pa) - real(r8), intent(in) :: cldfrac(ncol,nver) ! layer cloud fraction - real(r8), intent(in) :: tauc(nbnd,ncol,nver) ! cloud optical depth - real(r8), intent(out) :: taucmcl(ngpt,ncol,nver) ! subcolumn cloud optical depth [mcica] - - ! Local variables - - integer :: i, isubcol, k, n - - real(r8), parameter :: cldmin = 1.0e-80_r8 ! min cloud fraction - real(r8) :: cldf(ncol,nver) ! cloud fraction clipped to cldmin - - type(ShrKissRandGen) :: kiss_gen ! KISS RNG object - integer :: kiss_seed(ncol,4) - real(r8) :: rand_num_1d(ncol,1) ! random number (kissvec) - real(r8) :: rand_num(ncol,nver) ! random number (kissvec) - - real(r8) :: cdf(ngpt,ncol,nver) ! random numbers - logical :: iscloudy(ngpt,ncol,nver) ! flag that says whether a gridbox is cloudy - !------------------------------------------------------------------------------------------ - - ! clip cloud fraction - cldf(:,:) = cldfrac(:ncol,:) - where (cldf(:,:) < cldmin) - cldf(:,:) = 0._r8 - end where - - ! Create a seed that depends on the state of the columns. - ! Use pmid from bottom four layers. - do i = 1, ncol - kiss_seed(i,1) = (pmid(i,pver) - int(pmid(i,pver))) * 1000000000 - kiss_seed(i,2) = (pmid(i,pver-1) - int(pmid(i,pver-1))) * 1000000000 - kiss_seed(i,3) = (pmid(i,pver-2) - int(pmid(i,pver-2))) * 1000000000 - kiss_seed(i,4) = (pmid(i,pver-3) - int(pmid(i,pver-3))) * 1000000000 - end do - - ! create the RNG object - kiss_gen = ShrKissRandGen(kiss_seed) - - ! Advance randum number generator by changeseed values - do i = 1, changeSeed - call kiss_gen%random(rand_num_1d) - end do - - ! Generate random numbers in each subcolumn at every level - do isubcol = 1,ngpt - call kiss_gen%random(rand_num) - cdf(isubcol,:,:) = rand_num(:,:) - enddo - - ! Maximum-Random overlap - ! i) pick a random number for top layer. - ! ii) walk down the column: - ! - if the layer above is cloudy, use the same random number as in the layer above - ! - if the layer above is clear, use a new random number - - do k = 2, nver - do i = 1, ncol - do isubcol = 1, ngpt - if (cdf(isubcol,i,k-1) > 1._r8 - cldf(i,k-1) ) then - cdf(isubcol,i,k) = cdf(isubcol,i,k-1) - else - cdf(isubcol,i,k) = cdf(isubcol,i,k) * (1._r8 - cldf(i,k-1)) - end if - end do - end do - end do - - do k = 1, nver - iscloudy(:,:,k) = (cdf(:,:,k) >= 1._r8 - spread(cldf(:,k), dim=1, nCopies=ngpt) ) - end do - - ! -- generate subcolumns for homogeneous clouds ----- - ! where there is a cloud, set the subcolumn cloud properties; - ! incoming tauc should be in-cloud quantites and not grid-averaged quantities - do k = 1,nver - do i = 1,ncol - do isubcol = 1,ngpt - if (iscloudy(isubcol,i,k) .and. (cldf(i,k) > 0._r8) ) then - n = kdist%convert_gpt2band(isubcol) - taucmcl(isubcol,i,k) = tauc(n,i,k) - else - taucmcl(isubcol,i,k) = 0._r8 - end if - end do - end do - end do - - call kiss_gen%finalize() - -end subroutine mcica_subcol_lw - -!======================================================================================== - subroutine mcica_subcol_sw( & kdist, nbnd, ngpt, ncol, nlay, nver, changeseed, & pmid, cldfrac, tauc, ssac, asmc, & diff --git a/src/physics/rrtmgp/rad_solar_var.F90 b/src/physics/rrtmgp/rad_solar_var.F90 index ab608db7f9..de09ad84a4 100644 --- a/src/physics/rrtmgp/rad_solar_var.F90 +++ b/src/physics/rrtmgp/rad_solar_var.F90 @@ -7,7 +7,7 @@ module rad_solar_var use shr_kind_mod , only : r8 => shr_kind_r8 - use radconstants, only : nswbands, get_sw_spectral_boundaries, band2gpt_sw + use radiation_utils, only : get_sw_spectral_boundaries_ccpp use solar_irrad_data, only : sol_irrad, we, nbins, has_spectrum, sol_tsi use solar_irrad_data, only : do_spctrl_scaling use cam_abortutils, only : endrun @@ -29,10 +29,12 @@ module rad_solar_var contains !------------------------------------------------------------------------------- - subroutine rad_solar_var_init( ) + subroutine rad_solar_var_init(nswbands) + integer, intent(in) :: nswbands - integer :: ierr + integer :: ierr, errflg integer :: radmax_loc + character(len=512) :: errmsg if ( do_spctrl_scaling ) then @@ -55,7 +57,7 @@ subroutine rad_solar_var_init( ) call endrun('rad_solar_var_init: Error allocating space for irrad') end if - call get_sw_spectral_boundaries(radbinmin, radbinmax, 'nm') + call get_sw_spectral_boundaries_ccpp(radbinmin, radbinmax, 'nm', errmsg, errflg) ! Make sure that the far-IR is included, even if radiation grid does not ! extend that far down. 10^5 nm corresponds to a wavenumber of @@ -70,11 +72,13 @@ end subroutine rad_solar_var_init !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- - subroutine get_variability(toa_flux, sfac) + subroutine get_variability(toa_flux, sfac, band2gpt_sw, nswbands) ! Arguments real(r8), intent(in) :: toa_flux(:,:) ! TOA flux to be scaled (columns,gpts) real(r8), intent(out) :: sfac(:,:) ! scaling factors (columns,gpts) + integer, intent(in) :: band2gpt_sw(:,:) + integer, intent(in) :: nswbands ! Local variables integer :: i, j, istat, gpt_start, gpt_end, ncols diff --git a/src/physics/rrtmgp/radconstants.F90 b/src/physics/rrtmgp/radconstants.F90 index 3d4b47d09e..0edf9772e2 100644 --- a/src/physics/rrtmgp/radconstants.F90 +++ b/src/physics/rrtmgp/radconstants.F90 @@ -4,8 +4,9 @@ module radconstants ! code used in the RRTMGP model. use shr_kind_mod, only: r8 => shr_kind_r8 -use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use cam_abortutils, only: endrun +use radiation_utils, only: get_sw_spectral_boundaries_ccpp +use radiation_utils, only: get_lw_spectral_boundaries_ccpp implicit none private @@ -24,21 +25,16 @@ module radconstants real(r8), target :: wavenumber_low_longwave(nlwbands) real(r8), target :: wavenumber_high_longwave(nlwbands) -logical :: wavenumber_boundaries_set = .false. +logical :: wavenumber_boundaries_set = .true. ! First and last g-point for each band. integer, public, protected :: band2gpt_sw(2,nswbands) -integer, public, protected :: nswgpts ! number of SW g-points -integer, public, protected :: nlwgpts ! number of LW g-points - ! These are indices to specific bands for diagnostic output and COSP input. integer, public, protected :: idx_sw_diag = -1 ! band contains 500-nm wave integer, public, protected :: idx_nir_diag = -1 ! band contains 1000-nm wave integer, public, protected :: idx_uv_diag = -1 ! band contains 400-nm wave integer, public, protected :: idx_lw_diag = -1 ! band contains 1000 cm-1 wave (H20 window) -integer, public, protected :: idx_sw_cloudsim = -1 ! band contains 670-nm wave (for COSP) -integer, public, protected :: idx_lw_cloudsim = -1 ! band contains 10.5 micron wave (for COSP) ! GASES TREATED BY RADIATION (line spectra) ! These names are recognized by RRTMGP. They are in the coefficients files as @@ -53,7 +49,7 @@ module radconstants real(r8), public, parameter :: minmmr(nradgas) = epsilon(1._r8) public :: & - set_wavenumber_bands, & + radconstants_init, & get_sw_spectral_boundaries, & get_lw_spectral_boundaries, & get_band_index_by_value, & @@ -62,116 +58,38 @@ module radconstants !========================================================================================= contains !========================================================================================= - -subroutine set_wavenumber_bands(kdist_sw, kdist_lw) - - ! Set the low and high limits of the wavenumber grid for sw and lw. - ! Values come from RRTMGP coefficients datasets, and are stored in the - ! kdist objects. - ! - ! Set band indices for bands containing specific wavelengths. - - ! Arguments - type(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw - type(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw - - ! Local variables - integer :: istat - real(r8), allocatable :: values(:,:) - - character(len=128) :: errmsg - character(len=*), parameter :: sub = 'set_wavenumber_bands' - !---------------------------------------------------------------------------- - - ! Check that number of sw/lw bands in gas optics files matches the parameters. - if (kdist_sw%get_nband() /= nswbands) then - write(errmsg,'(a,i4,a,i4)') 'number of sw bands in file, ', kdist_sw%get_nband(), & - ", doesn't match parameter nswbands= ", nswbands - call endrun(sub//': ERROR: '//trim(errmsg)) - end if - if (kdist_lw%get_nband() /= nlwbands) then - write(errmsg,'(a,i4,a,i4)') 'number of lw bands in file, ', kdist_lw%get_nband(), & - ", doesn't match parameter nlwbands= ", nlwbands - call endrun(sub//': ERROR: '//trim(errmsg)) - end if - - nswgpts = kdist_sw%get_ngpt() - nlwgpts = kdist_lw%get_ngpt() - - ! SW band bounds in cm^-1 - allocate( values(2,nswbands), stat=istat ) - if (istat/=0) then - call endrun(sub//': ERROR allocating array: values(2,nswbands)') - end if - values = kdist_sw%get_band_lims_wavenumber() - wavenumber_low_shortwave = values(1,:) - wavenumber_high_shortwave = values(2,:) - - ! First and last g-point for each SW band: - band2gpt_sw = kdist_sw%get_band_lims_gpoint() - - ! Indices into specific bands - idx_sw_diag = get_band_index_by_value('sw', 500.0_r8, 'nm') - idx_nir_diag = get_band_index_by_value('sw', 1000.0_r8, 'nm') - idx_uv_diag = get_band_index_by_value('sw', 400._r8, 'nm') - idx_sw_cloudsim = get_band_index_by_value('sw', 0.67_r8, 'micron') - - deallocate(values) - - ! LW band bounds in cm^-1 - allocate( values(2,nlwbands), stat=istat ) - if (istat/=0) then - call endrun(sub//': ERROR allocating array: values(2,nlwbands)') - end if - values = kdist_lw%get_band_lims_wavenumber() - wavenumber_low_longwave = values(1,:) - wavenumber_high_longwave = values(2,:) - - ! Indices into specific bands - idx_lw_diag = get_band_index_by_value('lw', 1000.0_r8, 'cm^-1') - idx_lw_cloudsim = get_band_index_by_value('lw', 10.5_r8, 'micron') - - wavenumber_boundaries_set = .true. - -end subroutine set_wavenumber_bands - +subroutine radconstants_init(idx_sw_diag_in, idx_nir_diag_in, idx_uv_diag_in, idx_lw_diag_in) + integer, intent(in) :: idx_sw_diag_in + integer, intent(in) :: idx_nir_diag_in + integer, intent(in) :: idx_uv_diag_in + integer, intent(in) :: idx_lw_diag_in + + idx_sw_diag = idx_sw_diag_in + idx_nir_diag = idx_nir_diag_in + idx_uv_diag = idx_uv_diag_in + idx_lw_diag = idx_lw_diag_in + +end subroutine radconstants_init !========================================================================================= -subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units) + subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units) ! provide spectral boundaries of each shortwave band - real(r8), intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands) + real(r8), dimension(:), intent(out) :: low_boundaries + real(r8), dimension(:), intent(out) :: high_boundaries character(*), intent(in) :: units ! requested units - character(len=*), parameter :: sub = 'get_sw_spectral_boundaries' + character(len=512) :: errmsg + integer :: errflg !---------------------------------------------------------------------------- - if (.not. wavenumber_boundaries_set) then - call endrun(sub//': ERROR, wavenumber boundaries not set. ') + call get_sw_spectral_boundaries_ccpp(low_boundaries, high_boundaries, units, errmsg, errflg) + if (errflg /= 0) then + call endrun(errmsg) end if - select case (units) - case ('inv_cm','cm^-1','cm-1') - low_boundaries = wavenumber_low_shortwave - high_boundaries = wavenumber_high_shortwave - case('m','meter','meters') - low_boundaries = 1.e-2_r8/wavenumber_high_shortwave - high_boundaries = 1.e-2_r8/wavenumber_low_shortwave - case('nm','nanometer','nanometers') - low_boundaries = 1.e7_r8/wavenumber_high_shortwave - high_boundaries = 1.e7_r8/wavenumber_low_shortwave - case('um','micrometer','micrometers','micron','microns') - low_boundaries = 1.e4_r8/wavenumber_high_shortwave - high_boundaries = 1.e4_r8/wavenumber_low_shortwave - case('cm','centimeter','centimeters') - low_boundaries = 1._r8/wavenumber_high_shortwave - high_boundaries = 1._r8/wavenumber_low_shortwave - case default - call endrun(sub//': ERROR, requested spectral units not recognized: '//units) - end select - -end subroutine get_sw_spectral_boundaries + end subroutine get_sw_spectral_boundaries !========================================================================================= @@ -182,35 +100,17 @@ subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units) real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands) character(*), intent(in) :: units ! requested units - character(len=*), parameter :: sub = 'get_lw_spectral_boundaries' + character(len=512) :: errmsg + integer :: errflg !---------------------------------------------------------------------------- - if (.not. wavenumber_boundaries_set) then - call endrun(sub//': ERROR, wavenumber boundaries not set. ') + call get_lw_spectral_boundaries_ccpp(low_boundaries, high_boundaries, units, errmsg, errflg) + if (errflg /= 0) then + call endrun(errmsg) end if - select case (units) - case ('inv_cm','cm^-1','cm-1') - low_boundaries = wavenumber_low_longwave - high_boundaries = wavenumber_high_longwave - case('m','meter','meters') - low_boundaries = 1.e-2_r8/wavenumber_high_longwave - high_boundaries = 1.e-2_r8/wavenumber_low_longwave - case('nm','nanometer','nanometers') - low_boundaries = 1.e7_r8/wavenumber_high_longwave - high_boundaries = 1.e7_r8/wavenumber_low_longwave - case('um','micrometer','micrometers','micron','microns') - low_boundaries = 1.e4_r8/wavenumber_high_longwave - high_boundaries = 1.e4_r8/wavenumber_low_longwave - case('cm','centimeter','centimeters') - low_boundaries = 1._r8/wavenumber_high_longwave - high_boundaries = 1._r8/wavenumber_low_longwave - case default - call endrun(sub//': ERROR, requested spectral units not recognized: '//units) - end select - end subroutine get_lw_spectral_boundaries - + !========================================================================================= integer function rad_gas_index(gasname) diff --git a/src/physics/rrtmgp/radiation.F90 b/src/physics/rrtmgp/radiation.F90 index 79238fdb69..5f840e9a60 100644 --- a/src/physics/rrtmgp/radiation.F90 +++ b/src/physics/rrtmgp/radiation.F90 @@ -15,17 +15,15 @@ module radiation use physics_buffer, only: physics_buffer_desc, pbuf_add_field, dtype_r8, pbuf_get_index, & pbuf_set_field, pbuf_get_field, pbuf_old_tim_idx use camsrfexch, only: cam_out_t, cam_in_t -use physconst, only: cappa, cpair, gravit +use physconst, only: cappa, cpair, gravit, stebol use time_manager, only: get_nstep, is_first_step, is_first_restart_step, & get_curr_calday, get_step_size -use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_gas, rad_cnst_out +use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_out -use rrtmgp_inputs, only: rrtmgp_inputs_init - -use radconstants, only: nradgas, gasnamelength, gaslist, nswbands, nlwbands, & - nswgpts, set_wavenumber_bands +use radconstants, only: nradgas, gasnamelength, nswbands, nlwbands, & + gaslist, radconstants_init use rad_solar_var, only: rad_solar_var_init, get_variability use cloud_rad_props, only: cloud_rad_props_init @@ -48,16 +46,17 @@ module radiation pio_def_var, pio_put_var, pio_get_var, & pio_put_att, PIO_NOWRITE, pio_closefile -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_1scl, ty_optical_props_2str -use mo_source_functions, only: ty_source_func_lw -use mo_fluxes, only: ty_fluxes_broadband -use mo_fluxes_byband, only: ty_fluxes_byband +use ccpp_gas_concentrations, only: ty_gas_concs_ccpp +use ccpp_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp_ccpp +use ccpp_optical_props, only: ty_optical_props_1scl_ccpp, ty_optical_props_2str_ccpp +use ccpp_source_functions, only: ty_source_func_lw_ccpp +use ccpp_fluxes, only: ty_fluxes_broadband_ccpp +use ccpp_fluxes_byband, only: ty_fluxes_byband_ccpp use string_utils, only: to_lower use cam_abortutils, only: endrun, handle_allocate_error use cam_logfile, only: iulog +use rrtmgp_pre, only: radiation_do_ccpp implicit none @@ -66,8 +65,8 @@ module radiation public :: & radiation_readnl, &! read namelist variables - radiation_register, &! registers radiation physics buffer fields radiation_do, &! query which radiation calcs are done this timestep + radiation_register, &! registers radiation physics buffer fields radiation_init, &! initialization radiation_define_restart, &! define variables for restart radiation_write_restart, &! write variables to restart @@ -151,6 +150,17 @@ module radiation logical :: graupel_in_rad = .false. ! graupel in radiation code logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation +! Gathered indices of day and night columns +! chunk_column_index = IdxDay(daylight_column_index) +integer :: nday ! Number of daylight columns +integer :: nnite ! Number of night columns +integer :: idxday(pcols) ! chunk indices of daylight columns +integer :: idxnite(pcols) ! chunk indices of night columns +real(r8) :: coszrs(pcols) ! Cosine solar zenith angle +real(r8) :: eccf ! Earth orbit eccentricity factor + +integer :: band2gpt_sw(2,nswbands) + ! active_calls is set by a rad_constituents method after parsing namelist input ! for the rad_climate and rad_diag_N entries. logical :: active_calls(0:N_DIAG) @@ -170,6 +180,15 @@ module radiation integer :: cld_idx = 0 integer :: cldfsnow_idx = 0 integer :: cldfgrau_idx = 0 +integer :: dei_idx +integer :: mu_idx +integer :: lambda_idx +integer :: iciwp_idx +integer :: iclwp_idx +integer :: des_idx +integer :: icswp_idx +integer :: icgrauwp_idx +integer :: degrau_idx character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ',& '_d6 ','_d7 ','_d8 ','_d9 ','_d10'/) @@ -180,6 +199,8 @@ module radiation ! Number of layers in radiation calculations. integer :: nlay +! Number of interfaces in radiation calculations. +integer :: nlayp ! Number of CAM layers in radiation calculations. Is either equal to nlay, or is ! 1 less than nlay if "extra layer" is used in the radiation calculations. @@ -200,9 +221,27 @@ module radiation ! Note: for CAM's top to bottom indexing, the index of a given layer ! (midpoint) and the upper interface of that layer, are the same. +integer :: nlwgpts +integer :: nswgpts + +! Band indices for bands containing specific wavelengths +integer :: idx_sw_diag +integer :: idx_nir_diag +integer :: idx_uv_diag +integer :: idx_sw_cloudsim +integer :: idx_lw_diag +integer :: idx_lw_cloudsim + +real(r8) :: sw_low_bounds(nswbands) +real(r8) :: sw_high_bounds(nswbands) + +! Flag to perform shortwave or longwave on current timestep +logical :: dosw +logical :: dolw + ! Gas optics objects contain the data read from the coefficients files. -type(ty_gas_optics_rrtmgp) :: kdist_sw -type(ty_gas_optics_rrtmgp) :: kdist_lw +type(ty_gas_optics_rrtmgp_ccpp) :: kdist_sw +type(ty_gas_optics_rrtmgp_ccpp) :: kdist_lw ! lower case version of gaslist for RRTMGP character(len=gasnamelength) :: gaslist_lc(nradgas) @@ -343,33 +382,26 @@ end subroutine radiation_register !================================================================================================ -function radiation_do(op, timestep) +function radiation_do(op) ! Return true if the specified operation is done this timestep. character(len=*), intent(in) :: op ! name of operation - integer, intent(in), optional:: timestep logical :: radiation_do ! return value ! Local variables integer :: nstep ! current timestep number + integer :: errcode + character(len=512) :: errmsg !----------------------------------------------------------------------- - if (present(timestep)) then - nstep = timestep - else - nstep = get_nstep() - end if + nstep = get_nstep() select case (op) case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + call radiation_do_ccpp(op, nstep, iradsw, irad_always, radiation_do, errmsg, errcode) case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always + call radiation_do_ccpp(op, nstep, iradlw, irad_always, radiation_do, errmsg, errcode) case default call endrun('radiation_do: unknown operation:'//op) end select @@ -378,48 +410,11 @@ end function radiation_do !================================================================================================ -real(r8) function radiation_nextsw_cday() - - ! If a SW radiation calculation will be done on the next time-step, then return - ! the calendar day of that time-step. Otherwise return -1.0 - - ! Local variables - integer :: nstep ! timestep counter - logical :: dosw ! true => do shosrtwave calc - integer :: offset ! offset for calendar day calculation - integer :: dtime ! integer timestep size - real(r8):: caldayp1 ! calendar day of next time-step - - !----------------------------------------------------------------------- - - radiation_nextsw_cday = -1._r8 - dosw = .false. - nstep = get_nstep() - dtime = get_step_size() - offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') - end if - - ! determine if next radiation time-step not equal to next time-step - if (get_nstep() >= 1) then - caldayp1 = get_curr_calday(offset=int(dtime)) - if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8 - end if - -end function radiation_nextsw_cday - -!================================================================================================ - subroutine radiation_init(pbuf2d) + use rrtmgp_pre, only: rrtmgp_pre_init + use rrtmgp_inputs, only: rrtmgp_inputs_init + use rrtmgp_inputs_cam, only: rrtmgp_inputs_cam_init + use rrtmgp_lw_cloud_optics, only: rrtmgp_lw_cloud_optics_init ! Initialize the radiation and cloud optics. ! Add fields to the history buffer. @@ -428,11 +423,13 @@ subroutine radiation_init(pbuf2d) type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! local variables - character(len=128) :: errmsg + character(len=512) :: errmsg ! names of gases that are available in the model ! -- needed for the kdist initialization routines - type(ty_gas_concs) :: available_gases + type(ty_gas_concs_ccpp) :: available_gases + + real(r8) :: qrl_unused(1,1) integer :: i, icall integer :: nstep ! current timestep number @@ -442,68 +439,64 @@ subroutine radiation_init(pbuf2d) ! temperature, water vapor, cloud ice and cloud ! liquid budgets. integer :: history_budget_histfile_num ! history file number for budget fields - integer :: ierr, istat + integer :: ierr, istat, errflg + + ! Cloud optics variables + integer :: nmu, n_g_d, nlambda + real(kind=r8), allocatable :: abs_lw_ice(:,:) + real(kind=r8), allocatable :: abs_lw_liq(:,:,:) + real(kind=r8), allocatable :: g_lambda(:,:) + real(kind=r8), allocatable :: g_mu(:) + real(kind=r8), allocatable :: g_d_eff(:) + real(kind=r8) :: tiny integer :: dtime character(len=*), parameter :: sub = 'radiation_init' !----------------------------------------------------------------------- - ! Number of layers in radiation calculation is capped by the number of - ! pressure interfaces below 1 Pa. When the entire model atmosphere is - ! below 1 Pa then an extra layer is added to the top of the model for - ! the purpose of the radiation calculation. - - nlay = count( pref_edge(:) > 1._r8 ) ! pascals (0.01 mbar) - - if (nlay == pverp) then - ! Top model interface is below 1 Pa. RRTMGP is active in all model layers plus - ! 1 extra layer between model top and 1 Pa. - ktopcam = 1 - ktoprad = 2 - nlaycam = pver - else if (nlay == (pverp-1)) then - ! Special case nlay == (pverp-1) -- topmost interface outside bounds (CAM MT config), treat as if it is ok. - ktopcam = 1 - ktoprad = 2 - nlaycam = pver - nlay = nlay+1 ! reassign the value so later code understands to treat this case like nlay==pverp - write(iulog,*) 'RADIATION_INIT: Special case of 1 model interface at p < 1Pa. Top layer will be INCLUDED in radiation calculation.' - write(iulog,*) 'RADIATION_INIT: nlay = ',nlay, ' same as pverp: ',nlay==pverp - else - ! nlay < pverp. nlay layers are used in radiation calcs, and they are - ! all CAM layers. - ktopcam = pver - nlay + 1 - ktoprad = 1 - nlaycam = nlay + ! Initialize available_gases object + call rrtmgp_pre_init(nradgas, gaslist, available_gases, gaslist_lc, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) end if - - ! Create lowercase version of the gaslist for RRTMGP. The ty_gas_concs objects - ! work with CAM's uppercase names, but other objects that get input from the gas - ! concs objects don't work. - do i = 1, nradgas - gaslist_lc(i) = to_lower(gaslist(i)) - end do - - errmsg = available_gases%init(gaslist_lc) - call stop_on_err(errmsg, sub, 'available_gases%init') ! Read RRTMGP coefficients files and initialize kdist objects. call coefs_init(coefs_sw_file, available_gases, kdist_sw) call coefs_init(coefs_lw_file, available_gases, kdist_lw) - ! Set the sw/lw band boundaries in radconstants. Also sets - ! indicies of specific bands for diagnostic output and COSP input. - call set_wavenumber_bands(kdist_sw, kdist_lw) - call rad_solar_var_init() + ! Set up inputs to RRTMGP + call rrtmgp_inputs_init(ktopcam, ktoprad, nlaycam, sw_low_bounds, sw_high_bounds, nswbands, & + pref_edge, nlay, pver, pverp, kdist_sw, kdist_lw, qrl_unused, is_first_step(), use_rad_dt_cosz, & + get_step_size(), get_nstep(), iradsw, dt_avg, irad_always, is_first_restart_step(), masterproc, & + nlwbands, nradgas, gasnamelength, iulog, idx_sw_diag, idx_nir_diag, idx_uv_diag, & + idx_sw_cloudsim, idx_lw_diag, idx_lw_cloudsim, gaslist, nswgpts, nlwgpts, nlayp, & + nextsw_cday, get_curr_calday(), band2gpt_sw, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if + + ! Set up CAM-side RRTMGP inputs - will go away once SW radiation is CCPPized + call rrtmgp_inputs_cam_init(ktopcam, ktoprad, idx_sw_diag, idx_nir_diag, idx_uv_diag, idx_sw_cloudsim, idx_lw_diag, & + idx_lw_cloudsim) - ! The spectral band boundaries need to be set before this init is called. - call rrtmgp_inputs_init(ktopcam, ktoprad) + ! Set radconstants module-level index variables that we're setting in CCPP-ized scheme now + call radconstants_init(idx_sw_diag, idx_nir_diag, idx_uv_diag, idx_lw_diag) + + call rad_solar_var_init(nswbands) ! initialize output fields for offline driver call rad_data_init(pbuf2d) - call cloud_rad_props_init() + call cloud_rad_props_init(nmu, nlambda, n_g_d, abs_lw_liq, abs_lw_ice, & + g_mu, g_lambda, g_d_eff, tiny) + + call rrtmgp_lw_cloud_optics_init(nmu, nlambda, n_g_d, & + abs_lw_liq, abs_lw_ice, nlwbands, g_mu, g_lambda, & + g_d_eff, tiny, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if cld_idx = pbuf_get_index('CLD') cldfsnow_idx = pbuf_get_index('CLDFSNOW', errcode=ierr) @@ -513,30 +506,11 @@ subroutine radiation_init(pbuf2d) call pbuf_set_field(pbuf2d, qrl_idx, 0._r8) end if - ! Set the radiation timestep for cosz calculations if requested using - ! the adjusted iradsw value from radiation - if (use_rad_dt_cosz) then - dtime = get_step_size() - dt_avg = iradsw*dtime - end if - - ! Surface components to get radiation computed today - if (.not. is_first_restart_step()) then - nextsw_cday = get_curr_calday() - end if - call phys_getopts(history_amwg_out = history_amwg, & history_vdiag_out = history_vdiag, & history_budget_out = history_budget, & history_budget_histfile_num_out = history_budget_histfile_num) - ! "irad_always" is number of time steps to execute radiation continuously from - ! start of initial OR restart run - nstep = get_nstep() - if (irad_always > 0) then - irad_always = irad_always + nstep - end if - if (docosp) call cospsimulator_intr_init() allocate(cosp_cnt(begchunk:endchunk), stat=istat) @@ -833,25 +807,36 @@ subroutine radiation_tend( & !----------------------------------------------------------------------- ! Location/Orbital Parameters for cosine zenith angle - use phys_grid, only: get_rlat_all_p, get_rlon_all_p - use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr - use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz - - use rrtmgp_inputs, only: rrtmgp_set_state, rrtmgp_set_gases_lw, rrtmgp_set_cloud_lw, & - rrtmgp_set_aer_lw, rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, & - rrtmgp_set_aer_sw + use phys_grid, only: get_rlat_all_p, get_rlon_all_p + use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr + use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz + + ! CCPPized schemes + use rrtmgp_inputs, only: rrtmgp_inputs_run + use rrtmgp_pre, only: rrtmgp_pre_run, rrtmgp_pre_timestep_init + use rrtmgp_lw_cloud_optics, only: rrtmgp_lw_cloud_optics_run + use rrtmgp_lw_mcica_subcol_gen, only: rrtmgp_lw_mcica_subcol_gen_run + use rrtmgp_lw_gas_optics_pre, only: rrtmgp_lw_gas_optics_pre_run + use rrtmgp_lw_gas_optics, only: rrtmgp_lw_gas_optics_run + use rrtmgp_lw_main, only: rrtmgp_lw_main_run + use rrtmgp_dry_static_energy_tendency, only: rrtmgp_dry_static_energy_tendency_run + use rrtmgp_post, only: rrtmgp_post_run + + use rrtmgp_inputs_cam, only: rrtmgp_get_gas_mmrs, rrtmgp_set_aer_lw, & + rrtmgp_set_gases_sw, rrtmgp_set_cloud_sw, & + rrtmgp_set_aer_sw ! RRTMGP drivers for flux calculations. - use mo_rte_lw, only: rte_lw - use mo_rte_sw, only: rte_sw + use mo_rte_lw, only: rte_lw + use mo_rte_sw, only: rte_sw - use radheat, only: radheat_tend + use radheat, only: radheat_tend - use radiation_data, only: rad_data_write + use radiation_data, only: rad_data_write - use interpolate_data, only: vertinterp - use tropopause, only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE - use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps + use interpolate_data, only: vertinterp + use tropopause, only: tropopause_find_cam, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE + use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps ! Arguments @@ -870,7 +855,7 @@ subroutine radiation_tend( & ! if the argument is not present logical :: write_output - integer :: i, k, istat + integer :: i, k, gas_idx, istat integer :: lchnk, ncol logical :: dosw, dolw integer :: icall ! loop index for climate/diagnostic radiation calls @@ -882,19 +867,20 @@ subroutine radiation_tend( & real(r8) :: clon(pcols) ! current longitudes(radians) real(r8) :: coszrs(pcols) ! Cosine solar zenith angle - ! Gathered indices of day and night columns - ! chunk_column_index = IdxDay(daylight_column_index) - integer :: Nday ! Number of daylight columns - integer :: Nnite ! Number of night columns - integer :: IdxDay(pcols) ! chunk indices of daylight columns - integer :: IdxNite(pcols) ! chunk indices of night columns - integer :: itim_old + integer :: nextsw_nstep + integer :: offset + real(r8) :: next_cday real(r8), pointer :: cld(:,:) ! cloud fraction real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds" real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds" + real(r8), pointer :: cldfsnow_in(:,:) ! Cloud fraction of just "snow clouds", subset + real(r8), pointer :: cldfgrau_in(:,:) ! Cloud fraction of just "graupel clouds", subset real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction + real(r8) :: cld_lw_abs(nlwbands,state%ncol,pver) ! Cloud absorption optics depth + real(r8) :: snow_lw_abs(nlwbands,state%ncol,pver) ! Snow absorption optics depth + real(r8) :: grau_lw_abs(nlwbands,state%ncol,pver) ! Graupel absorption optics depth real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate real(r8), pointer :: fsds(:) ! Surface solar down flux @@ -903,6 +889,16 @@ subroutine radiation_tend( & real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top + real(r8), pointer :: dei(:,:) + real(r8), pointer :: mu(:,:) + real(r8), pointer :: lambda(:,:) + real(r8), pointer :: iciwp(:,:) + real(r8), pointer :: iclwp(:,:) + real(r8), pointer :: des(:,:) + real(r8), pointer :: icswp(:,:) + real(r8), pointer :: icgrauwp(:,:) + real(r8), pointer :: degrau(:,:) + real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up @@ -924,6 +920,10 @@ subroutine radiation_tend( & real(r8), allocatable :: coszrs_day(:) real(r8), allocatable :: alb_dir(:,:) real(r8), allocatable :: alb_dif(:,:) + real(r8), allocatable :: tauc(:,:,:) + real(r8), allocatable :: cldf(:,:) + + real(r8), allocatable :: gas_mmrs(:,:,:) ! in-cloud optical depths for COSP real(r8) :: cld_tau_cloudsim(pcols,pver) ! liq + ice @@ -936,42 +936,44 @@ subroutine radiation_tend( & ! Set vertical indexing in RRTMGP to be the same as CAM (top to bottom). logical, parameter :: top_at_1 = .true. + logical :: do_graupel, do_snow + ! TOA solar flux on RRTMGP g-points real(r8), allocatable :: toa_flux(:,:) ! Scale factors based on spectral distribution from input irradiance dataset real(r8), allocatable :: sfac(:,:) ! Planck sources for LW. - type(ty_source_func_lw) :: sources_lw + type(ty_source_func_lw_ccpp) :: sources_lw ! Gas volume mixing ratios. Use separate objects for LW and SW because SW only does ! calculations for daylight columns. ! These objects have a final method which deallocates the internal memory when they ! go out of scope (i.e., when radiation_tend returns), so no need for explicit deallocation. - type(ty_gas_concs) :: gas_concs_lw - type(ty_gas_concs) :: gas_concs_sw + type(ty_gas_concs_ccpp) :: gas_concs_lw + type(ty_gas_concs_ccpp) :: gas_concs_sw ! Atmosphere optics. This object is initialized with gas optics, then is incremented ! by the aerosol optics for the clear-sky radiative flux calculations, and then ! incremented again by the cloud optics for the all-sky radiative flux calculations. - type(ty_optical_props_1scl) :: atm_optics_lw - type(ty_optical_props_2str) :: atm_optics_sw + type(ty_optical_props_1scl_ccpp) :: atm_optics_lw + type(ty_optical_props_2str_ccpp) :: atm_optics_sw ! Cloud optical properties objects (McICA sampling of cloud optical properties). - type(ty_optical_props_1scl) :: cloud_lw - type(ty_optical_props_2str) :: cloud_sw + type(ty_optical_props_1scl_ccpp) :: cloud_lw + type(ty_optical_props_2str_ccpp) :: cloud_sw ! Aerosol optical properties objects. - type(ty_optical_props_1scl) :: aer_lw - type(ty_optical_props_2str) :: aer_sw + type(ty_optical_props_1scl_ccpp) :: aer_lw + type(ty_optical_props_2str_ccpp) :: aer_sw ! Flux objects contain all fluxes computed by RRTMGP. ! SW allsky fluxes always include spectrally resolved fluxes needed for surface models. - type(ty_fluxes_byband) :: fsw + type(ty_fluxes_byband_ccpp) :: fsw ! LW allsky fluxes only need spectrally resolved fluxes when spectralflux=.true. - type(ty_fluxes_byband) :: flw + type(ty_fluxes_byband_ccpp) :: flw ! Only broadband fluxes needed for clear sky (diagnostics). - type(ty_fluxes_broadband) :: fswc, flwc + type(ty_fluxes_broadband_ccpp) :: fswc, flwc ! Arrays for output diagnostics on CAM grid. real(r8) :: fns(pcols,pverp) ! net shortwave flux @@ -979,6 +981,10 @@ subroutine radiation_tend( & real(r8) :: fnl(pcols,pverp) ! net longwave flux real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux + ! Unused variables for rte_lw + real(r8) :: fluxlwup_jac(1,1) + real(r8) :: lw_ds(1,1) + ! for COSP real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau @@ -987,6 +993,7 @@ subroutine radiation_tend( & real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables character(len=128) :: errmsg + integer :: errflg, err character(len=*), parameter :: sub = 'radiation_tend' !-------------------------------------------------------------------------------------- @@ -1002,9 +1009,6 @@ subroutine radiation_tend( & write_output = .true. end if - dosw = radiation_do('sw', get_nstep()) ! do shortwave radiation calc this timestep? - dolw = radiation_do('lw', get_nstep()) ! do longwave radiation calc this timestep? - ! Cosine solar zenith angle for current time step calday = get_curr_calday() call get_rlat_all_p(lchnk, ncol, clat) @@ -1025,18 +1029,24 @@ subroutine radiation_tend( & end do end if - ! Gather night/day column indices. - Nday = 0 - Nnite = 0 - do i = 1, ncol - if ( coszrs(i) > 0.0_r8 ) then - Nday = Nday + 1 - IdxDay(Nday) = i - else - Nnite = Nnite + 1 - IdxNite(Nnite) = i - end if - end do + ! Get next SW radiation timestep + call rrtmgp_pre_timestep_init(get_nstep(), get_step_size(), iradsw, irad_always, offset, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if + + ! Calculate next calendar day and next radiation calendar day + nextsw_cday = get_curr_calday(offset=offset) + next_cday = get_curr_calday(offset=int(get_step_size())) + + ! Determine if we're running radiation (sw and/or lw) this timestep, + ! find daylight and nighttime indices, and initialize fluxes + call rrtmgp_pre_run(coszrs(:ncol), get_nstep(), get_step_size(), iradsw, iradlw, irad_always, & + ncol, next_cday, idxday, nday, idxnite, nnite, dosw, dolw, nlay, nlwbands, & + nswbands, spectralflux, nextsw_cday, fsw, fswc, flw, flwc, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if ! Associate pointers to physics buffer fields itim_old = pbuf_old_tim_idx() @@ -1067,12 +1077,6 @@ subroutine radiation_tend( & call pbuf_get_field(pbuf, ld_idx, ld) end if - ! Allocate the flux arrays and init to zero. - call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fsw, do_direct=.true.) - call initialize_rrtmgp_fluxes(nday, nlay+1, nswbands, fswc, do_direct=.true.) - call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flw) - call initialize_rrtmgp_fluxes(ncol, nlay+1, nlwbands, flwc) - ! For CRM, make cloud equal to input observations: if (scm_crm_mode .and. have_cld) then do k = 1, pver @@ -1080,6 +1084,11 @@ subroutine radiation_tend( & end do end if + !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists + flns(:) = 0._r8 + flnt(:) = 0._r8 + !REMOVECAM_END + ! Find tropopause height if needed for diagnostic output if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists @@ -1090,10 +1099,6 @@ subroutine radiation_tend( & backup=TROP_ALG_CLIMATE) end if - ! Get time of next radiation calculation - albedos will need to be - ! calculated by each surface model at this time - nextsw_cday = radiation_nextsw_cday() - if (dosw .or. dolw) then allocate( & @@ -1102,23 +1107,45 @@ subroutine radiation_tend( & t_rad(ncol,nlay), pmid_rad(ncol,nlay), pint_rad(ncol,nlay+1), & t_day(nday,nlay), pmid_day(nday,nlay), pint_day(nday,nlay+1), & coszrs_day(nday), alb_dir(nswbands,nday), alb_dif(nswbands,nday), & - stat=istat) + cldf(ncol,nlaycam), tauc(nlwbands,ncol,nlaycam), stat=istat) call handle_allocate_error(istat, sub, 't_sfc,..,alb_dif') + allocate(gas_mmrs(ncol, pver, nradgas), stat=istat, errmsg=errmsg) + if (errflg /= 0) then + call handle_allocate_error(istat, sub, 'gas_mmrs, message: '//errmsg) + end if + + if (associated(cldfgrau)) then + cldfgrau_in => cldfgrau(:ncol,:) + else + cldfgrau_in => null() + end if - ! Prepares state variables, daylit columns, albedos for RRTMGP - call rrtmgp_set_state( & - state, cam_in, ncol, nlay, nday, & - idxday, coszrs, kdist_sw, t_sfc, emis_sfc, & - t_rad, pmid_rad, pint_rad, t_day, pmid_day, & - pint_day, coszrs_day, alb_dir, alb_dif) + if (associated(cldfsnow)) then + cldfsnow_in => cldfsnow(:ncol,:) + else + cldfsnow_in => null() + end if + ! Prepare state variables, daylit columns, albedos for RRTMGP + ! Also calculate modified cloud fraction + call rrtmgp_inputs_run(dosw, dolw, associated(cldfsnow), associated(cldfgrau), & + state%pmid(:ncol,:), state%pint(:ncol,:), state%t(:ncol,:), & + nday, idxday, cldfprime(:ncol,:), coszrs(:ncol), kdist_sw, t_sfc, & + emis_sfc, t_rad, pmid_rad, pint_rad, t_day, pmid_day, & + pint_day, coszrs_day, alb_dir, alb_dif, cam_in%lwup(:ncol), stebol, & + ncol, ktopcam, ktoprad, nswbands, cam_in%asdir(:ncol), cam_in%asdif(:ncol), & + sw_low_bounds, sw_high_bounds, cam_in%aldir(:ncol), cam_in%aldif(:ncol), nlay, & + pverp, pver, cld(:ncol,:), cldfsnow_in, cldfgrau_in, & + graupel_in_rad, gasnamelength, gaslist, gas_concs_lw, aer_lw, atm_optics_lw, & + kdist_lw, sources_lw, aer_sw, atm_optics_sw, gas_concs_sw, & + errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if ! Output the mass per layer, and total column burdens for gas and aerosol ! constituents in the climate list. call rad_cnst_out(0, state, pbuf) - ! Modified cloud fraction accounts for radiatively active snow and/or graupel - call modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime) - !========================! ! SHORTWAVE calculations ! !========================! @@ -1127,7 +1154,7 @@ subroutine radiation_tend( & ! Set cloud optical properties in cloud_sw object. call rrtmgp_set_cloud_sw( & - state, pbuf, nlay, nday, idxday, & + state, pbuf, nlay, nday, idxday, nswgpts, & nnite, idxnite, pmid_day, cld, cldfsnow, & cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, & rd%tot_cld_vistau, rd%tot_icld_vistau, rd%liq_icld_vistau, & @@ -1138,25 +1165,6 @@ subroutine radiation_tend( & call radiation_output_cld(lchnk, rd) end if - ! If no daylight columns, can't create empty RRTMGP objects - if (nday > 0) then - - ! Initialize object for gas concentrations. - errmsg = gas_concs_sw%init(gaslist_lc) - call stop_on_err(errmsg, sub, 'gas_concs_sw%init') - - ! Initialize object for combined gas + aerosol + cloud optics. - ! Allocates arrays for properties represented on g-points. - errmsg = atm_optics_sw%alloc_2str(nday, nlay, kdist_sw) - call stop_on_err(errmsg, sub, 'atm_optics_sw%alloc_2str') - - ! Initialize object for SW aerosol optics. Allocates arrays - ! for properties represented by band. - errmsg = aer_sw%alloc_2str(nday, nlay, kdist_sw%get_band_lims_wavenumber()) - call stop_on_err(errmsg, sub, 'aer_sw%alloc_2str') - - end if - ! The climate (icall==0) calculation must occur last. do icall = N_DIAG, 0, -1 if (active_calls(icall)) then @@ -1170,17 +1178,17 @@ subroutine radiation_tend( & ! Compute the gas optics (stored in atm_optics_sw). ! toa_flux is the reference solar source from RRTMGP data. - !$acc data copyin(kdist_sw,pmid_day,pint_day,t_day,gas_concs_sw) & - !$acc copy(atm_optics_sw) & + !$acc data copyin(kdist_sw%gas_props,pmid_day,pint_day,t_day,gas_concs_sw%gas_concs) & + !$acc copy(atm_optics_sw%optical_props) & !$acc copyout(toa_flux) - errmsg = kdist_sw%gas_optics( & - pmid_day, pint_day, t_day, gas_concs_sw, atm_optics_sw, & + errmsg = kdist_sw%gas_props%gas_optics( & + pmid_day, pint_day, t_day, gas_concs_sw%gas_concs, atm_optics_sw%optical_props, & toa_flux) + call stop_on_err(errmsg, sub, 'kdist_sw%gas_props%gas_optics') !$acc end data - call stop_on_err(errmsg, sub, 'kdist_sw%gas_optics') ! Scale the solar source - call get_variability(toa_flux, sfac) + call get_variability(toa_flux, sfac, band2gpt_sw, nswbands) toa_flux = toa_flux * sfac * eccf end if @@ -1195,31 +1203,31 @@ subroutine radiation_tend( & ! Increment the gas optics (in atm_optics_sw) by the aerosol optics in aer_sw. !$acc data copyin(coszrs_day, toa_flux, alb_dir, alb_dif, & - !$acc atm_optics_sw, atm_optics_sw%tau, & - !$acc atm_optics_sw%ssa, atm_optics_sw%g, & - !$acc aer_sw, aer_sw%tau, & - !$acc aer_sw%ssa, aer_sw%g, & - !$acc cloud_sw, cloud_sw%tau, & - !$acc cloud_sw%ssa, cloud_sw%g) & - !$acc copy(fswc, fswc%flux_net,fswc%flux_up,fswc%flux_dn, & - !$acc fsw, fsw%flux_net, fsw%flux_up, fsw%flux_dn) - errmsg = aer_sw%increment(atm_optics_sw) - call stop_on_err(errmsg, sub, 'aer_sw%increment') + !$acc atm_optics_sw%optical_props, atm_optics_sw%optical_props%tau, & + !$acc atm_optics_sw%optical_props%ssa, atm_optics_sw%optical_props%g, & + !$acc aer_sw%optical_props, aer_sw%optical_props%tau, & + !$acc aer_sw%optical_props%ssa, aer_sw%optical_props%g, & + !$acc cloud_sw%optical_props, cloud_sw%optical_props%tau, & + !$acc cloud_sw%optical_props%ssa, cloud_sw%optical_props%g) & + !$acc copy(fswc%fluxes, fswc%fluxes%flux_net,fswc%fluxes%flux_up,fswc%fluxes%flux_dn, & + !$acc fsw%fluxes, fsw%fluxes%flux_net, fsw%fluxes%flux_up, fsw%fluxes%flux_dn) + errmsg = aer_sw%optical_props%increment(atm_optics_sw%optical_props) + call stop_on_err(errmsg, sub, 'aer_sw%optical_props%increment') ! Compute clear-sky fluxes. errmsg = rte_sw(& - atm_optics_sw, top_at_1, coszrs_day, toa_flux, & - alb_dir, alb_dif, fswc) + atm_optics_sw%optical_props, top_at_1, coszrs_day, toa_flux, & + alb_dir, alb_dif, fswc%fluxes) call stop_on_err(errmsg, sub, 'clear-sky rte_sw') ! Increment the aerosol+gas optics (in atm_optics_sw) by the cloud optics in cloud_sw. - errmsg = cloud_sw%increment(atm_optics_sw) - call stop_on_err(errmsg, sub, 'cloud_sw%increment') + errmsg = cloud_sw%optical_props%increment(atm_optics_sw%optical_props) + call stop_on_err(errmsg, sub, 'cloud_sw%optical_props%increment') ! Compute all-sky fluxes. errmsg = rte_sw(& - atm_optics_sw, top_at_1, coszrs_day, toa_flux, & - alb_dir, alb_dif, fsw) + atm_optics_sw%optical_props, top_at_1, coszrs_day, toa_flux, & + alb_dir, alb_dif, fsw%fluxes) call stop_on_err(errmsg, sub, 'all-sky rte_sw') !$acc end data end if @@ -1233,11 +1241,6 @@ subroutine radiation_tend( & end if ! (active_calls(icall)) end do ! loop over diagnostic calcs (icall) - - else - ! SW calc not done. pbuf carries Q*dp across timesteps. - ! Convert to Q before calling radheat_tend. - qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:) end if ! if (dosw) !=======================! @@ -1246,76 +1249,107 @@ subroutine radiation_tend( & if (dolw) then - ! Initialize object for Planck sources. - errmsg = sources_lw%alloc(ncol, nlay, kdist_lw) - call stop_on_err(errmsg, sub, 'sources_lw%alloc') - - ! Set cloud optical properties in cloud_lw object. - call rrtmgp_set_cloud_lw( & - state, pbuf, ncol, nlay, nlaycam, & - cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, & - kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim) + ! Grab additional pbuf fields for LW cloud optics + dei_idx = pbuf_get_index('DEI',errcode=err) + mu_idx = pbuf_get_index('MU',errcode=err) + lambda_idx = pbuf_get_index('LAMBDAC',errcode=err) + iciwp_idx = pbuf_get_index('ICIWP',errcode=err) + iclwp_idx = pbuf_get_index('ICLWP',errcode=err) + des_idx = pbuf_get_index('DES',errcode=err) + icswp_idx = pbuf_get_index('ICSWP',errcode=err) + icgrauwp_idx = pbuf_get_index('ICGRAUWP',errcode=err) ! Available when using MG3 + degrau_idx = pbuf_get_index('DEGRAU',errcode=err) ! Available when using MG3 + call pbuf_get_field(pbuf, lambda_idx, lambda) + call pbuf_get_field(pbuf, mu_idx, mu) + call pbuf_get_field(pbuf, iclwp_idx, iclwp) + call pbuf_get_field(pbuf, iciwp_idx, iciwp) + call pbuf_get_field(pbuf, dei_idx, dei) + call pbuf_get_field(pbuf, icswp_idx, icswp) + call pbuf_get_field(pbuf, des_idx, des) + if (icgrauwp_idx > 0) then + call pbuf_get_field(pbuf, icgrauwp_idx, icgrauwp) + end if + if (degrau_idx > 0) then + call pbuf_get_field(pbuf, degrau_idx, degrau) + end if - ! Initialize object for gas concentrations - errmsg = gas_concs_lw%init(gaslist_lc) - call stop_on_err(errmsg, sub, 'gas_concs_lw%init') + do_graupel = ((icgrauwp_idx > 0) .and. (degrau_idx > 0) .and. associated(cldfgrau)) .and. graupel_in_rad + do_snow = associated(cldfsnow) - ! Initialize object for combined gas + aerosol + cloud optics. - errmsg = atm_optics_lw%alloc_1scl(ncol, nlay, kdist_lw) - call stop_on_err(errmsg, sub, 'atm_optics_lw%alloc_1scl') + ! Set cloud optical properties in cloud_lw object. + call rrtmgp_lw_cloud_optics_run(dolw, ncol, nlay, nlaycam, cld(:ncol,:), cldfsnow_in, & + cldfgrau_in, cldfprime(:ncol,:), kdist_lw, cloud_lw, lambda(:ncol,:), & + mu(:ncol,:), iclwp(:ncol,:), iciwp(:ncol,:), dei(:ncol,:), icswp(:ncol,:), des(:ncol,:), & + icgrauwp(:ncol,:), degrau(:ncol,:), nlwbands, do_snow, do_graupel, pver, & + ktopcam, tauc, cldf, cld_lw_abs, snow_lw_abs, grau_lw_abs, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if - ! Initialize object for LW aerosol optics. - errmsg = aer_lw%alloc_1scl(ncol, nlay, kdist_lw%get_band_lims_wavenumber()) - call stop_on_err(errmsg, sub, 'aer_lw%alloc_1scl') + ! Cloud optics for COSP + cld_lw_abs_cloudsim(:ncol,:) = cld_lw_abs(idx_lw_cloudsim,:,:) + snow_lw_abs_cloudsim(:ncol,:) = snow_lw_abs(idx_lw_cloudsim,:,:) + grau_lw_abs_cloudsim(:ncol,:) = grau_lw_abs(idx_lw_cloudsim,:,:) + + ! Create McICA stochastic arrays for lw cloud optical properties + call rrtmgp_lw_mcica_subcol_gen_run(dolw, ktoprad, & + kdist_lw, nlwbands, nlwgpts, ncol, pver, nlaycam, nlwgpts, & + state%pmid(:ncol,:), cldf, tauc, cloud_lw, errmsg, errflg ) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if ! The climate (icall==0) calculation must occur last. do icall = N_DIAG, 0, -1 if (active_calls(icall)) then - ! Set gas volume mixing ratios for this call in gas_concs_lw. - call rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs_lw) + ! Grab the gas mass mixing ratios from rad_constituents + gas_mmrs = 0._r8 + call rrtmgp_get_gas_mmrs(icall, state, pbuf, nlay, gas_mmrs) + + ! Set gas volume mixing ratios for this call in gas_concs_lw + call rrtmgp_lw_gas_optics_pre_run(gas_mmrs, state%pmid(:ncol,:), state%pint(:ncol,:), nlay, ncol, gaslist, & + idxday, pverp, ktoprad, ktopcam, dolw, nradgas, gas_concs_lw, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if ! Compute the gas optics and Planck sources. - !$acc data copyin(kdist_lw, pmid_rad, pint_rad, & - !$acc t_rad, t_sfc, gas_concs_lw) & - !$acc copy(atm_optics_lw, atm_optics_lw%tau, & - !$acc sources_lw, sources_lw%lay_source, & - !$acc sources_lw%sfc_source, sources_lw%lev_source_inc, & - !$acc sources_lw%lev_source_dec, sources_lw%sfc_source_jac) - errmsg = kdist_lw%gas_optics( & - pmid_rad, pint_rad, t_rad, t_sfc, gas_concs_lw, & - atm_optics_lw, sources_lw) + !$acc data copyin(kdist_lw%gas_props, pmid_rad, pint_rad, & + !$acc t_rad, t_sfc, gas_concs_lw%gas_concs) & + !$acc copy(atm_optics_lw%optical_props, atm_optics_lw%optical_props%tau, & + !$acc sources_lw%sources, sources_lw%sources%lay_source, & + !$acc sources_lw%sources%sfc_source, sources_lw%sources%lev_source_inc, & + !$acc sources_lw%sources%lev_source_dec, sources_lw%sources%sfc_source_jac) + call rrtmgp_lw_gas_optics_run(dolw, 1, ncol, ncol, pmid_rad, pint_rad, t_rad, & + t_sfc, gas_concs_lw, atm_optics_lw, sources_lw, t_rad, .false., kdist_lw, errmsg, & + errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if !$acc end data - call stop_on_err(errmsg, sub, 'kdist_lw%gas_optics') ! Set LW aerosol optical properties in the aer_lw object. call rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw) - - ! Increment the gas optics by the aerosol optics. - !$acc data copyin(atm_optics_lw, atm_optics_lw%tau, & - !$acc aer_lw, aer_lw%tau, & - !$acc cloud_lw, cloud_lw%tau, & - !$acc sources_lw, sources_lw%lay_source, & - !$acc sources_lw%sfc_source, sources_lw%lev_source_inc, & - !$acc sources_lw%lev_source_dec, sources_lw%sfc_source_Jac, & + + ! Call the main rrtmgp_lw driver + !$acc data copyin(atm_optics_lw%optical_props, atm_optics_lw%optical_props%tau, & + !$acc aer_lw%optical_props, aer_lw%optical_props%tau, & + !$acc cloud_lw%optical_props, cloud_lw%optical_props%tau, & + !$acc sources_lw%sources, sources_lw%sources%lay_source, & + !$acc sources_lw%sources%sfc_source, sources_lw%sources%lev_source_inc, & + !$acc sources_lw%sources%lev_source_dec, sources_lw%sources%sfc_source_Jac, & !$acc emis_sfc) & - !$acc copy(flwc, flwc%flux_net, flwc%flux_up, flwc%flux_dn, & - !$acc flw, flw%flux_net, flw%flux_up, flw%flux_dn) - errmsg = aer_lw%increment(atm_optics_lw) - call stop_on_err(errmsg, sub, 'aer_lw%increment') - - ! Compute clear-sky LW fluxes - errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flwc) - call stop_on_err(errmsg, sub, 'clear-sky rte_lw') - - ! Increment the gas+aerosol optics by the cloud optics. - errmsg = cloud_lw%increment(atm_optics_lw) - call stop_on_err(errmsg, sub, 'cloud_lw%increment') - - ! Compute all-sky LW fluxes - errmsg = rte_lw(atm_optics_lw, top_at_1, sources_lw, emis_sfc, flw) - call stop_on_err(errmsg, sub, 'all-sky rte_lw') + !$acc copy(flwc%fluxes, flwc%fluxes%flux_net, flwc%fluxes%flux_up, flwc%fluxes%flux_dn, & + !$acc flw, flw%fluxes%flux_net, flw%fluxes%flux_up, flw%fluxes%flux_dn) + call rrtmgp_lw_main_run(dolw, dolw, .false., .false., .false., & + 0, ncol, 1, ncol, atm_optics_lw, & + cloud_lw, top_at_1, sources_lw, emis_sfc, kdist_lw, & + aer_lw, fluxlwup_jac, lw_ds, flwc, flw, errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if !$acc end data ! Transform RRTMGP outputs to CAM outputs and compute heating rates. @@ -1327,11 +1361,6 @@ subroutine radiation_tend( & end if ! (active_calls(icall)) end do ! loop over diagnostic calcs (icall) - - else - ! LW calc not done. pbuf carries Q*dp across timesteps. - ! Convert to Q before calling radheat_tend. - qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:) end if ! if (dolw) deallocate( & @@ -1386,16 +1415,15 @@ subroutine radiation_tend( & cosp_cnt(lchnk) = 0 end if end if ! docosp - - else - ! Radiative flux calculations not done. The quantity Q*dp is carried by the - ! physics buffer across timesteps. It must be converted to Q (dry static energy - ! tendency) before being passed to radheat_tend. - qrs(:ncol,:) = qrs(:ncol,:) / state%pdel(:ncol,:) - qrl(:ncol,:) = qrl(:ncol,:) / state%pdel(:ncol,:) - end if ! if (dosw .or. dolw) then + ! Calculate dry static energy if LW calc or SW calc wasn't done; needed before calling radheat_run + call rrtmgp_dry_static_energy_tendency_run(state%pdel(:ncol,:), (.not. dosw), (.not. dolw), & + qrs(:ncol,:), qrl(:ncol,:), errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if + ! Output for PORT: Parallel Offline Radiative Transport call rad_data_write(pbuf, state, cam_in, coszrs) @@ -1415,27 +1443,18 @@ subroutine radiation_tend( & call outfld('HR', ftem, pcols, lchnk) end if - ! The radiative heating rates are carried in the physics buffer across timesteps - ! as Q*dp (for energy conservation). - qrs(:ncol,:) = qrs(:ncol,:) * state%pdel(:ncol,:) - qrl(:ncol,:) = qrl(:ncol,:) * state%pdel(:ncol,:) - - cam_out%netsw(:ncol) = fsns(:ncol) - if (.not. present(rd_out)) then deallocate(rd) end if - call free_optics_sw(atm_optics_sw) - call free_optics_sw(cloud_sw) - call free_optics_sw(aer_sw) - call free_fluxes(fsw) - call free_fluxes(fswc) - - call sources_lw%finalize() - call free_optics_lw(cloud_lw) - call free_optics_lw(aer_lw) - call free_fluxes(flw) - call free_fluxes(flwc) + + cam_out%netsw(:) = 0._r8 + + ! Calculate radiative heating (Q*dp), set netsw flux, and do object cleanup + call rrtmgp_post_run(qrs(:ncol,:), qrl(:ncol,:), fsns(:ncol), state%pdel(:ncol,:), atm_optics_sw, cloud_sw, aer_sw, & + fsw, fswc, sources_lw, cloud_lw, aer_lw, flw, flwc, cam_out%netsw(:ncol), errmsg, errflg) + if (errflg /= 0) then + call endrun(sub//': '//errmsg) + end if !------------------------------------------------------------------------------- contains @@ -1448,9 +1467,9 @@ subroutine set_sw_diags() ! full chunks for output to CAM history. integer :: i - real(r8), dimension(size(fsw%bnd_flux_dn,1), & - size(fsw%bnd_flux_dn,2), & - size(fsw%bnd_flux_dn,3)) :: flux_dn_diffuse + real(r8), dimension(size(fsw%fluxes%bnd_flux_dn,1), & + size(fsw%fluxes%bnd_flux_dn,2), & + size(fsw%fluxes%bnd_flux_dn,3)) :: flux_dn_diffuse !------------------------------------------------------------------------- ! Initialize to provide 0.0 values for night columns. @@ -1475,18 +1494,18 @@ subroutine set_sw_diags() rd%fsntc = 0._r8 do i = 1, nday - fns(idxday(i),ktopcam:) = fsw%flux_net(i, ktoprad:) - fcns(idxday(i),ktopcam:) = fswc%flux_net(i,ktoprad:) - fsds(idxday(i)) = fsw%flux_dn(i, nlay+1) - rd%fsdsc(idxday(i)) = fswc%flux_dn(i, nlay+1) - rd%fsutoa(idxday(i)) = fsw%flux_up(i, 1) - rd%fsntoa(idxday(i)) = fsw%flux_net(i, 1) - rd%fsntoac(idxday(i)) = fswc%flux_net(i, 1) - rd%solin(idxday(i)) = fswc%flux_dn(i, 1) - rd%flux_sw_up(idxday(i),ktopcam:) = fsw%flux_up(i,ktoprad:) - rd%flux_sw_dn(idxday(i),ktopcam:) = fsw%flux_dn(i,ktoprad:) - rd%flux_sw_clr_up(idxday(i),ktopcam:) = fswc%flux_up(i,ktoprad:) - rd%flux_sw_clr_dn(idxday(i),ktopcam:) = fswc%flux_dn(i,ktoprad:) + fns(idxday(i),ktopcam:) = fsw%fluxes%flux_net(i, ktoprad:) + fcns(idxday(i),ktopcam:) = fswc%fluxes%flux_net(i,ktoprad:) + fsds(idxday(i)) = fsw%fluxes%flux_dn(i, nlay+1) + rd%fsdsc(idxday(i)) = fswc%fluxes%flux_dn(i, nlay+1) + rd%fsutoa(idxday(i)) = fsw%fluxes%flux_up(i, 1) + rd%fsntoa(idxday(i)) = fsw%fluxes%flux_net(i, 1) + rd%fsntoac(idxday(i)) = fswc%fluxes%flux_net(i, 1) + rd%solin(idxday(i)) = fswc%fluxes%flux_dn(i, 1) + rd%flux_sw_up(idxday(i),ktopcam:) = fsw%fluxes%flux_up(i,ktoprad:) + rd%flux_sw_dn(idxday(i),ktopcam:) = fsw%fluxes%flux_dn(i,ktoprad:) + rd%flux_sw_clr_up(idxday(i),ktopcam:) = fswc%fluxes%flux_up(i,ktoprad:) + rd%flux_sw_clr_dn(idxday(i),ktopcam:) = fswc%fluxes%flux_dn(i,ktoprad:) end do ! Compute heating rate as a dry static energy tendency. @@ -1511,8 +1530,8 @@ subroutine set_sw_diags() su = 0._r8 sd = 0._r8 do i = 1, nday - su(idxday(i),ktopcam:,:) = fsw%bnd_flux_up(i,ktoprad:,:) - sd(idxday(i),ktopcam:,:) = fsw%bnd_flux_dn(i,ktoprad:,:) + su(idxday(i),ktopcam:,:) = fsw%fluxes%bnd_flux_up(i,ktoprad:,:) + sd(idxday(i),ktopcam:,:) = fsw%fluxes%bnd_flux_dn(i,ktoprad:,:) end do end if @@ -1531,14 +1550,14 @@ subroutine set_sw_diags() cam_out%solld = 0.0_r8 ! Calculate diffuse flux from total and direct - flux_dn_diffuse = fsw%bnd_flux_dn - fsw%bnd_flux_dn_dir + flux_dn_diffuse = fsw%fluxes%bnd_flux_dn - fsw%fluxes%bnd_flux_dn_dir do i = 1, nday - cam_out%soll(idxday(i)) = sum(fsw%bnd_flux_dn_dir(i,nlay+1,1:9)) & - + 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10) + cam_out%soll(idxday(i)) = sum(fsw%fluxes%bnd_flux_dn_dir(i,nlay+1,1:9)) & + + 0.5_r8 * fsw%fluxes%bnd_flux_dn_dir(i,nlay+1,10) - cam_out%sols(idxday(i)) = 0.5_r8 * fsw%bnd_flux_dn_dir(i,nlay+1,10) & - + sum(fsw%bnd_flux_dn_dir(i,nlay+1,11:14)) + cam_out%sols(idxday(i)) = 0.5_r8 * fsw%fluxes%bnd_flux_dn_dir(i,nlay+1,10) & + + sum(fsw%fluxes%bnd_flux_dn_dir(i,nlay+1,11:14)) cam_out%solld(idxday(i)) = sum(flux_dn_diffuse(i,nlay+1,1:9)) & + 0.5_r8 * flux_dn_diffuse(i,nlay+1,10) @@ -1560,13 +1579,13 @@ subroutine set_lw_diags() fcnl = 0._r8 ! RTE-RRTMGP convention for net is (down - up) **CAM assumes (up - down) !! - fnl(:ncol,ktopcam:) = -1._r8 * flw%flux_net( :, ktoprad:) - fcnl(:ncol,ktopcam:) = -1._r8 * flwc%flux_net( :, ktoprad:) + fnl(:ncol,ktopcam:) = -1._r8 * flw%fluxes%flux_net( :, ktoprad:) + fcnl(:ncol,ktopcam:) = -1._r8 * flwc%fluxes%flux_net( :, ktoprad:) - rd%flux_lw_up(:ncol,ktopcam:) = flw%flux_up( :, ktoprad:) - rd%flux_lw_clr_up(:ncol,ktopcam:) = flwc%flux_up(:, ktoprad:) - rd%flux_lw_dn(:ncol,ktopcam:) = flw%flux_dn( :, ktoprad:) - rd%flux_lw_clr_dn(:ncol,ktopcam:) = flwc%flux_dn(:, ktoprad:) + rd%flux_lw_up(:ncol,ktopcam:) = flw%fluxes%flux_up( :, ktoprad:) + rd%flux_lw_clr_up(:ncol,ktopcam:) = flwc%fluxes%flux_up(:, ktoprad:) + rd%flux_lw_dn(:ncol,ktopcam:) = flw%fluxes%flux_dn( :, ktoprad:) + rd%flux_lw_clr_dn(:ncol,ktopcam:) = flwc%fluxes%flux_dn(:, ktoprad:) call heating_rate('LW', ncol, fnl, qrl) call heating_rate('LW', ncol, fcnl, rd%qrlc) @@ -1577,11 +1596,11 @@ subroutine set_lw_diags() rd%flnsc(:ncol) = fcnl(:ncol, pverp) rd%flntc(:ncol) = fcnl(:ncol, ktopcam) ! net lw flux at top-of-model - cam_out%flwds(:ncol) = flw%flux_dn(:, nlay+1) - rd%fldsc(:ncol) = flwc%flux_dn(:, nlay+1) + cam_out%flwds(:ncol) = flw%fluxes%flux_dn(:, nlay+1) + rd%fldsc(:ncol) = flwc%fluxes%flux_dn(:, nlay+1) - rd%flut(:ncol) = flw%flux_up(:, ktoprad) - rd%flutc(:ncol) = flwc%flux_up(:, ktoprad) + rd%flut(:ncol) = flw%fluxes%flux_up(:, ktoprad) + rd%flutc(:ncol) = flwc%fluxes%flux_up(:, ktoprad) ! Output fluxes at 200 mb call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200) @@ -1595,8 +1614,8 @@ subroutine set_lw_diags() if (spectralflux) then lu = 0._r8 ld = 0._r8 - lu(:ncol, ktopcam:, :) = flw%bnd_flux_up(:, ktoprad:, :) - ld(:ncol, ktopcam:, :) = flw%bnd_flux_dn(:, ktoprad:, :) + lu(:ncol, ktopcam:, :) = flw%fluxes%bnd_flux_up(:, ktoprad:, :) + ld(:ncol, ktopcam:, :) = flw%fluxes%bnd_flux_dn(:, ktoprad:, :) end if end subroutine set_lw_diags @@ -1796,14 +1815,15 @@ end subroutine radiation_output_lw !=============================================================================== subroutine coefs_init(coefs_file, available_gases, kdist) + use rrtmgp_lw_gas_optics_data, only: rrtmgp_lw_gas_optics_data_init ! Read data from coefficients file. Initialize the kdist object. ! available_gases object provides the gas names that CAM provides. ! arguments character(len=*), intent(in) :: coefs_file - class(ty_gas_concs), intent(in) :: available_gases - class(ty_gas_optics_rrtmgp), intent(out) :: kdist + class(ty_gas_concs_ccpp), intent(in) :: available_gases + class(ty_gas_optics_rrtmgp_ccpp), intent(inout) :: kdist ! local variables type(file_desc_t) :: fh ! pio file handle @@ -1868,6 +1888,7 @@ subroutine coefs_init(coefs_file, available_gases, kdist) fit_coeffs character(len=128) :: error_msg + character(len=512) :: errmsg character(len=*), parameter :: sub = 'coefs_init' !---------------------------------------------------------------------------- @@ -2275,25 +2296,22 @@ subroutine coefs_init(coefs_file, available_gases, kdist) ! on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave) if (allocated(totplnk) .and. allocated(planck_frac)) then - error_msg = kdist%load( & - available_gases, gas_names, key_species, & - band2gpt, band_lims_wavenum, & - press_ref, press_ref_trop, temp_ref, & - temp_ref_p, temp_ref_t, vmr_ref, & - kmajor, kminor_lower, kminor_upper, & - gas_minor, identifier_minor, & - minor_gases_lower, minor_gases_upper, & - minor_limits_gpt_lower, minor_limits_gpt_upper, & - minor_scales_with_density_lower, & - minor_scales_with_density_upper, & - scaling_gas_lower, scaling_gas_upper, & - scale_by_complement_lower, scale_by_complement_upper, & - kminor_start_lower, kminor_start_upper, & - totplnk, planck_frac, rayl_lower, rayl_upper, & - optimal_angle_fit) + call rrtmgp_lw_gas_optics_data_init(kdist, available_gases, gas_names, & + key_species, band2gpt, band_lims_wavenum, press_ref, press_ref_trop, & + temp_ref, temp_ref_p, temp_ref_t, vmr_ref, kmajor, kminor_lower, & + kminor_upper, gas_minor, identifier_minor, minor_gases_lower, & + minor_gases_upper, minor_limits_gpt_lower, minor_limits_gpt_upper, & + minor_scales_with_density_lower, minor_scales_with_density_upper, & + scaling_gas_lower, scaling_gas_upper, scale_by_complement_lower, & + scale_by_complement_upper, kminor_start_lower, kminor_start_upper, & + totplnk, planck_frac, rayl_lower, rayl_upper, optimal_angle_fit, & + errmsg, ierr) + if (ierr /= 0) then + call endrun(sub//': '//errmsg) + end if else if (allocated(solar_src_quiet)) then - error_msg = kdist%load( & - available_gases, gas_names, key_species, & + error_msg = kdist%gas_props%load( & + available_gases%gas_concs, gas_names, key_species, & band2gpt, band_lims_wavenum, & press_ref, press_ref_trop, temp_ref, & temp_ref_p, temp_ref_t, vmr_ref, & @@ -2315,7 +2333,7 @@ subroutine coefs_init(coefs_file, available_gases, kdist) error_msg = 'must supply either totplnk and planck_frac, or solar_src_[*]' end if - call stop_on_err(error_msg, sub, 'kdist%load') + call stop_on_err(error_msg, sub, 'kdist%gas_props%load') deallocate( & gas_names, key_species, & @@ -2346,176 +2364,6 @@ end subroutine coefs_init !========================================================================================= -subroutine initialize_rrtmgp_fluxes(ncol, nlevels, nbands, fluxes, do_direct) - - ! Allocate flux arrays and set values to zero. - - ! Arguments - integer, intent(in) :: ncol, nlevels, nbands - class(ty_fluxes_broadband), intent(inout) :: fluxes - logical, optional, intent(in) :: do_direct - - ! Local variables - logical :: do_direct_local - integer :: istat - character(len=*), parameter :: sub = 'initialize_rrtmgp_fluxes' - !---------------------------------------------------------------------------- - - if (present(do_direct)) then - do_direct_local = .true. - else - do_direct_local = .false. - end if - - ! Broadband fluxes - allocate(fluxes%flux_up(ncol, nlevels), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%flux_up') - allocate(fluxes%flux_dn(ncol, nlevels), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%flux_dn') - allocate(fluxes%flux_net(ncol, nlevels), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%flux_net') - if (do_direct_local) then - allocate(fluxes%flux_dn_dir(ncol, nlevels), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%flux_dn_dir') - end if - - select type (fluxes) - type is (ty_fluxes_byband) - ! Fluxes by band always needed for SW. Only allocate for LW - ! when spectralflux is true. - if (nbands == nswbands .or. spectralflux) then - allocate(fluxes%bnd_flux_up(ncol, nlevels, nbands), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_up') - allocate(fluxes%bnd_flux_dn(ncol, nlevels, nbands), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn') - allocate(fluxes%bnd_flux_net(ncol, nlevels, nbands), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_net') - if (do_direct_local) then - allocate(fluxes%bnd_flux_dn_dir(ncol, nlevels, nbands), stat=istat) - call handle_allocate_error(istat, sub, 'fluxes%bnd_flux_dn_dir') - end if - end if - end select - - ! Initialize - call reset_fluxes(fluxes) - -end subroutine initialize_rrtmgp_fluxes - -!========================================================================================= - -subroutine reset_fluxes(fluxes) - - ! Reset flux arrays to zero. - - class(ty_fluxes_broadband), intent(inout) :: fluxes - !---------------------------------------------------------------------------- - - ! Reset broadband fluxes - fluxes%flux_up(:,:) = 0._r8 - fluxes%flux_dn(:,:) = 0._r8 - fluxes%flux_net(:,:) = 0._r8 - if (associated(fluxes%flux_dn_dir)) fluxes%flux_dn_dir(:,:) = 0._r8 - - select type (fluxes) - type is (ty_fluxes_byband) - ! Reset band-by-band fluxes - if (associated(fluxes%bnd_flux_up)) fluxes%bnd_flux_up(:,:,:) = 0._r8 - if (associated(fluxes%bnd_flux_dn)) fluxes%bnd_flux_dn(:,:,:) = 0._r8 - if (associated(fluxes%bnd_flux_net)) fluxes%bnd_flux_net(:,:,:) = 0._r8 - if (associated(fluxes%bnd_flux_dn_dir)) fluxes%bnd_flux_dn_dir(:,:,:) = 0._r8 - end select - -end subroutine reset_fluxes - -!========================================================================================= - -subroutine free_optics_sw(optics) - - type(ty_optical_props_2str), intent(inout) :: optics - - if (allocated(optics%tau)) deallocate(optics%tau) - if (allocated(optics%ssa)) deallocate(optics%ssa) - if (allocated(optics%g)) deallocate(optics%g) - call optics%finalize() - -end subroutine free_optics_sw - -!========================================================================================= - -subroutine free_optics_lw(optics) - - type(ty_optical_props_1scl), intent(inout) :: optics - - if (allocated(optics%tau)) deallocate(optics%tau) - call optics%finalize() - -end subroutine free_optics_lw - -!========================================================================================= - -subroutine free_fluxes(fluxes) - - class(ty_fluxes_broadband), intent(inout) :: fluxes - - if (associated(fluxes%flux_up)) deallocate(fluxes%flux_up) - if (associated(fluxes%flux_dn)) deallocate(fluxes%flux_dn) - if (associated(fluxes%flux_net)) deallocate(fluxes%flux_net) - if (associated(fluxes%flux_dn_dir)) deallocate(fluxes%flux_dn_dir) - - select type (fluxes) - type is (ty_fluxes_byband) - if (associated(fluxes%bnd_flux_up)) deallocate(fluxes%bnd_flux_up) - if (associated(fluxes%bnd_flux_dn)) deallocate(fluxes%bnd_flux_dn) - if (associated(fluxes%bnd_flux_net)) deallocate(fluxes%bnd_flux_net) - if (associated(fluxes%bnd_flux_dn_dir)) deallocate(fluxes%bnd_flux_dn_dir) - end select - -end subroutine free_fluxes - -!========================================================================================= - -subroutine modified_cloud_fraction(ncol, cld, cldfsnow, cldfgrau, cldfprime) - - ! Compute modified cloud fraction, cldfprime. - ! 1. initialize as cld - ! 2. modify for snow if cldfsnow is available. use max(cld, cldfsnow) - ! 3. modify for graupel if cldfgrau is available and graupel_in_rad is true. - ! use max(cldfprime, cldfgrau) - - ! Arguments - integer, intent(in) :: ncol - real(r8), pointer :: cld(:,:) ! cloud fraction - real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds" - real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds" - real(r8), intent(out) :: cldfprime(:,:) ! modified cloud fraction - - ! Local variables - integer :: i, k - !---------------------------------------------------------------------------- - - if (associated(cldfsnow)) then - do k = 1, pver - do i = 1, ncol - cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k)) - end do - end do - else - cldfprime(:ncol,:) = cld(:ncol,:) - end if - - if (associated(cldfgrau) .and. graupel_in_rad) then - do k = 1, pver - do i = 1, ncol - cldfprime(i,k) = max(cldfprime(i,k), cldfgrau(i,k)) - end do - end do - end if - -end subroutine modified_cloud_fraction - -!========================================================================================= - subroutine stop_on_err(errmsg, sub, info) ! call endrun if RRTMGP function returns non-empty error message. diff --git a/src/physics/rrtmgp/rrtmgp_inputs.F90 b/src/physics/rrtmgp/rrtmgp_inputs_cam.F90 similarity index 61% rename from src/physics/rrtmgp/rrtmgp_inputs.F90 rename to src/physics/rrtmgp/rrtmgp_inputs_cam.F90 index 4f73ae9029..4c65ffbb69 100644 --- a/src/physics/rrtmgp/rrtmgp_inputs.F90 +++ b/src/physics/rrtmgp/rrtmgp_inputs_cam.F90 @@ -1,4 +1,4 @@ -module rrtmgp_inputs +module rrtmgp_inputs_cam !-------------------------------------------------------------------------------- ! Transform data for inputs from CAM's data structures to those used by @@ -17,9 +17,7 @@ module rrtmgp_inputs use physics_buffer, only: physics_buffer_desc use camsrfexch, only: cam_in_t -use radconstants, only: nradgas, gaslist, nswbands, nlwbands, nswgpts, nlwgpts, & - get_sw_spectral_boundaries, idx_sw_diag, idx_sw_cloudsim, & - idx_lw_cloudsim +use radconstants, only: nradgas, gaslist, nswbands, nlwbands use rad_constituents, only: rad_cnst_get_gas @@ -28,29 +26,27 @@ module rrtmgp_inputs get_snow_optics_sw, snow_cloud_get_rad_props_lw, & get_grau_optics_sw, grau_cloud_get_rad_props_lw -use mcica_subcol_gen, only: mcica_subcol_sw, mcica_subcol_lw +use mcica_subcol_gen, only: mcica_subcol_sw use aer_rad_props, only: aer_rad_props_sw, aer_rad_props_lw -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_2str, ty_optical_props_1scl +use ccpp_gas_concentrations, only: ty_gas_concs_ccpp +use ccpp_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp_ccpp +use ccpp_optical_props, only: ty_optical_props_2str_ccpp, ty_optical_props_1scl_ccpp use cam_history_support, only: fillvalue -use cam_logfile, only: iulog use cam_abortutils, only: endrun use error_messages, only: alloc_err +use radiation_utils, only: get_sw_spectral_boundaries_ccpp implicit none private save public :: & - rrtmgp_inputs_init, & - rrtmgp_set_state, & - rrtmgp_set_gases_lw, & + rrtmgp_inputs_cam_init, & + rrtmgp_get_gas_mmrs, & rrtmgp_set_gases_sw, & - rrtmgp_set_cloud_lw, & rrtmgp_set_cloud_sw, & rrtmgp_set_aer_lw, & rrtmgp_set_aer_sw @@ -59,15 +55,16 @@ module rrtmgp_inputs ! This value is to match the arbitrary small value used in RRTMG to decide ! when a quantity is effectively zero. real(r8), parameter :: tiny = 1.0e-80_r8 - -! Indices for copying data between cam and rrtmgp arrays -integer :: ktopcam ! Index in CAM arrays of top level (layer or interface) at which - ! RRTMGP is active. -integer :: ktoprad ! Index in RRTMGP arrays of the layer or interface corresponding - ! to CAM's top layer or interface - -! wavenumber (cm^-1) boundaries of shortwave bands -real(r8) :: sw_low_bounds(nswbands), sw_high_bounds(nswbands) +real(r8) :: sw_low_bounds(nswbands) +real(r8) :: sw_high_bounds(nswbands) +integer :: ktopcam +integer :: ktoprad +integer :: idx_sw_diag +integer :: idx_nir_diag +integer :: idx_uv_diag +integer :: idx_sw_cloudsim +integer :: idx_lw_diag +integer :: idx_lw_cloudsim ! Mapping from RRTMG shortwave bands to RRTMGP. Currently needed to continue using ! the SW optics datasets from RRTMG (even thought there is a slight mismatch in the @@ -79,186 +76,37 @@ module rrtmgp_inputs contains !================================================================================================== -subroutine rrtmgp_inputs_init(ktcam, ktrad) +!================================================================================================== +subroutine rrtmgp_inputs_cam_init(ktcam, ktrad, idx_sw_diag_in, idx_nir_diag_in, idx_uv_diag_in, & + idx_sw_cloudsim_in, idx_lw_diag_in, idx_lw_cloudsim_in) ! Note that this routine must be called after the calls to set_wavenumber_bands which set ! the sw/lw band boundaries in the radconstants module. integer, intent(in) :: ktcam integer, intent(in) :: ktrad + integer, intent(in) :: idx_sw_diag_in + integer, intent(in) :: idx_nir_diag_in + integer, intent(in) :: idx_uv_diag_in + integer, intent(in) :: idx_sw_cloudsim_in + integer, intent(in) :: idx_lw_diag_in + integer, intent(in) :: idx_lw_cloudsim_in + character(len=512) :: errmsg + integer :: errflg ktopcam = ktcam ktoprad = ktrad + idx_sw_diag = idx_sw_diag_in + idx_nir_diag = idx_nir_diag_in + idx_uv_diag = idx_uv_diag_in + idx_sw_cloudsim = idx_sw_cloudsim_in + idx_lw_diag = idx_lw_diag_in + idx_lw_cloudsim = idx_lw_cloudsim_in ! Initialize the module data containing the SW band boundaries. - call get_sw_spectral_boundaries(sw_low_bounds, sw_high_bounds, 'cm^-1') + call get_sw_spectral_boundaries_ccpp(sw_low_bounds, sw_high_bounds, 'cm^-1', errmsg, errflg) -end subroutine rrtmgp_inputs_init - -!========================================================================================= - -subroutine rrtmgp_set_state( & - state, cam_in, ncol, nlay, nday, & - idxday, coszrs, kdist_sw, t_sfc, emis_sfc, & - t_rad, pmid_rad, pint_rad, t_day, pmid_day, & - pint_day, coszrs_day, alb_dir, alb_dif) - - ! arguments - type(physics_state), intent(in) :: state ! CAM physics state - type(cam_in_t), intent(in) :: cam_in ! CAM import state - integer, intent(in) :: ncol ! # cols in CAM chunk - integer, intent(in) :: nlay ! # layers in rrtmgp grid - integer, intent(in) :: nday ! # daylight columns - integer, intent(in) :: idxday(:) ! chunk indicies of daylight columns - real(r8), intent(in) :: coszrs(:) ! cosine of solar zenith angle - class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw ! spectral information - - real(r8), intent(out) :: t_sfc(ncol) ! surface temperature [K] - real(r8), intent(out) :: emis_sfc(nlwbands,ncol) ! emissivity at surface [] - real(r8), intent(out) :: t_rad(ncol,nlay) ! layer midpoint temperatures [K] - real(r8), intent(out) :: pmid_rad(ncol,nlay) ! layer midpoint pressures [Pa] - real(r8), intent(out) :: pint_rad(ncol,nlay+1) ! layer interface pressures [Pa] - real(r8), intent(out) :: t_day(nday,nlay) ! layer midpoint temperatures [K] - real(r8), intent(out) :: pmid_day(nday,nlay) ! layer midpoint pressure [Pa] - real(r8), intent(out) :: pint_day(nday,nlay+1) ! layer interface pressures [Pa] - real(r8), intent(out) :: coszrs_day(nday) ! cosine of solar zenith angle - real(r8), intent(out) :: alb_dir(nswbands,nday) ! surface albedo, direct radiation - real(r8), intent(out) :: alb_dif(nswbands,nday) ! surface albedo, diffuse radiation - - ! local variables - integer :: i, k, iband - - real(r8) :: tref_min, tref_max - - character(len=*), parameter :: sub='rrtmgp_set_state' - character(len=512) :: errmsg - !-------------------------------------------------------------------------------- - - t_sfc = sqrt(sqrt(cam_in%lwup(:ncol)/stebol)) ! Surface temp set based on longwave up flux. - - ! Set surface emissivity to 1.0. - ! The land model *does* have its own surface emissivity, but is not spectrally resolved. - ! The LW upward flux is calculated with that land emissivity, and the "radiative temperature" - ! t_sfc is derived from that flux. We assume, therefore, that the emissivity is unity - ! to be consistent with t_sfc. - emis_sfc(:,:) = 1._r8 - - ! Level ordering is the same for both CAM and RRTMGP (top to bottom) - t_rad(:,ktoprad:) = state%t(:ncol,ktopcam:) - pmid_rad(:,ktoprad:) = state%pmid(:ncol,ktopcam:) - pint_rad(:,ktoprad:) = state%pint(:ncol,ktopcam:) - - ! Add extra layer values if needed. - if (nlay == pverp) then - t_rad(:,1) = state%t(:ncol,1) - ! The top reference pressure from the RRTMGP coefficients datasets is 1.005183574463 Pa - ! Set the top of the extra layer just below that. - pint_rad(:,1) = 1.01_r8 - - ! next interface down in LT will always be > 1Pa - ! but in MT we apply adjustment to have it be 1.02 Pa if it was too high - where (pint_rad(:,2) <= pint_rad(:,1)) pint_rad(:,2) = pint_rad(:,1)+0.01_r8 - - ! set the highest pmid (in the "extra layer") to the midpoint (guarantees > 1Pa) - pmid_rad(:,1) = pint_rad(:,1) + 0.5_r8 * (pint_rad(:,2) - pint_rad(:,1)) - - ! For case of CAM MT, also ensure pint_rad(:,2) > pint_rad(:,1) & pmid_rad(:,2) > max(pmid_rad(:,1), min_pressure) - where (pmid_rad(:,2) <= kdist_sw%get_press_min()) pmid_rad(:,2) = pint_rad(:,2) + 0.01_r8 - else - ! nlay < pverp, thus the 1 Pa level is within a CAM layer. Assuming the top interface of - ! this layer is at a pressure < 1 Pa, we need to adjust the top of this layer so that it - ! is within the valid pressure range of RRTMGP (otherwise RRTMGP issues an error). Then - ! set the midpoint pressure halfway between the interfaces. - pint_rad(:,1) = 1.01_r8 - pmid_rad(:,1) = 0.5_r8 * (pint_rad(:,1) + pint_rad(:,2)) - end if - - ! Limit temperatures to be within the limits of RRTMGP validity. - tref_min = kdist_sw%get_temp_min() - tref_max = kdist_sw%get_temp_max() - t_rad = merge(t_rad, tref_min, t_rad > tref_min) - t_rad = merge(t_rad, tref_max, t_rad < tref_max) - - ! Construct arrays containing only daylight columns - do i = 1, nday - t_day(i,:) = t_rad(idxday(i),:) - pmid_day(i,:) = pmid_rad(idxday(i),:) - pint_day(i,:) = pint_rad(idxday(i),:) - coszrs_day(i) = coszrs(idxday(i)) - end do - - ! Assign albedos to the daylight columns (from E3SM implementation) - ! Albedos are imported from the surface models as broadband (visible, and near-IR), - ! and we need to map these to appropriate narrower bands used in RRTMGP. Bands - ! are categorized broadly as "visible/UV" or "infrared" based on wavenumber. - ! Loop over bands, and determine for each band whether it is broadly in the - ! visible or infrared part of the spectrum based on a dividing line of - ! 0.7 micron, or 14286 cm^-1 - do iband = 1,nswbands - if (is_visible(sw_low_bounds(iband)) .and. & - is_visible(sw_high_bounds(iband))) then - - ! Entire band is in the visible - do i = 1, nday - alb_dir(iband,i) = cam_in%asdir(idxday(i)) - alb_dif(iband,i) = cam_in%asdif(idxday(i)) - end do - - else if (.not.is_visible(sw_low_bounds(iband)) .and. & - .not.is_visible(sw_high_bounds(iband))) then - ! Entire band is in the longwave (near-infrared) - do i = 1, nday - alb_dir(iband,i) = cam_in%aldir(idxday(i)) - alb_dif(iband,i) = cam_in%aldif(idxday(i)) - end do - else - ! Band straddles the visible to near-infrared transition, so we take - ! the albedo to be the average of the visible and near-infrared - ! broadband albedos - do i = 1, nday - alb_dir(iband,i) = 0.5_r8 * (cam_in%aldir(idxday(i)) + cam_in%asdir(idxday(i))) - alb_dif(iband,i) = 0.5_r8 * (cam_in%aldif(idxday(i)) + cam_in%asdif(idxday(i))) - end do - end if - end do - - ! Strictly enforce albedo bounds - where (alb_dir < 0) - alb_dir = 0.0_r8 - end where - where (alb_dir > 1) - alb_dir = 1.0_r8 - end where - where (alb_dif < 0) - alb_dif = 0.0_r8 - end where - where (alb_dif > 1) - alb_dif = 1.0_r8 - end where - -end subroutine rrtmgp_set_state - -!========================================================================================= - -pure logical function is_visible(wavenumber) - - ! Wavenumber is in the visible if it is above the visible threshold - ! wavenumber, and in the infrared if it is below the threshold - ! This function doesn't distinquish between visible and UV. - - ! wavenumber in inverse cm (cm^-1) - real(r8), intent(in) :: wavenumber - - ! Set threshold between visible and infrared to 0.7 micron, or 14286 cm^-1 - real(r8), parameter :: visible_wavenumber_threshold = 14286._r8 ! cm^-1 - - if (wavenumber > visible_wavenumber_threshold) then - is_visible = .true. - else - is_visible = .false. - end if - -end function is_visible +end subroutine rrtmgp_inputs_cam_init !========================================================================================= @@ -319,7 +167,7 @@ subroutine rad_gas_get_vmr(icall, gas_name, state, pbuf, nlay, numactivecols, ga integer, intent(in) :: nlay ! number of layers in radiation calculation integer, intent(in) :: numactivecols ! number of columns, ncol for LW, nday for SW - type(ty_gas_concs), intent(inout) :: gas_concs ! the result is VRM inside gas_concs + type(ty_gas_concs_ccpp), intent(inout) :: gas_concs ! the result is VRM inside gas_concs integer, optional, intent(in) :: idxday(:) ! indices of daylight columns in a chunk @@ -402,7 +250,7 @@ subroutine rad_gas_get_vmr(icall, gas_name, state, pbuf, nlay, numactivecols, ga end do end if - errmsg = gas_concs%set_vmr(gas_name, gas_vmr) + errmsg = gas_concs%gas_concs%set_vmr(gas_name, gas_vmr) if (len_trim(errmsg) > 0) then call endrun(sub//': ERROR, gas_concs%set_vmr: '//trim(errmsg)) end if @@ -414,32 +262,29 @@ end subroutine rad_gas_get_vmr !================================================================================================== -subroutine rrtmgp_set_gases_lw(icall, state, pbuf, nlay, gas_concs) - - ! Set gas vmr for the gases in the radconstants module's gaslist. +subroutine rrtmgp_get_gas_mmrs(icall, state, pbuf, nlay, gas_mmrs) - ! The memory management for the gas_concs object is internal. The arrays passed to it - ! are copied to the internally allocated memory. Each call to the set_vmr method checks - ! whether the gas already has memory allocated, and if it does that memory is deallocated - ! and new memory is allocated. + ! Retrieve mass mixing ratios for radiatively active gases from rad_constituents ! arguments integer, intent(in) :: icall ! index of climate/diagnostic radiation call type(physics_state), target, intent(in) :: state type(physics_buffer_desc), pointer :: pbuf(:) integer, intent(in) :: nlay - type(ty_gas_concs), intent(inout) :: gas_concs + real(r8), intent(out) :: gas_mmrs(:,:,:) ! local variables integer :: i, ncol - character(len=*), parameter :: sub = 'rrtmgp_set_gases_lw' + real(r8), pointer :: gas_mmr(:,:) + character(len=*), parameter :: sub = 'rrtmgp_get_gas_mmrs' !-------------------------------------------------------------------------------- ncol = state%ncol do i = 1, nradgas - call rad_gas_get_vmr(icall, gaslist(i), state, pbuf, nlay, ncol, gas_concs) + call rad_cnst_get_gas(icall, gaslist(i), state, pbuf, gas_mmr) + gas_mmrs(:,:,i) = gas_mmr end do -end subroutine rrtmgp_set_gases_lw +end subroutine rrtmgp_get_gas_mmrs !================================================================================================== @@ -457,7 +302,7 @@ subroutine rrtmgp_set_gases_sw( & integer, intent(in) :: nlay integer, intent(in) :: nday integer, intent(in) :: idxday(:) - type(ty_gas_concs), intent(inout) :: gas_concs + type(ty_gas_concs_ccpp), intent(inout) :: gas_concs ! local variables integer :: i @@ -473,143 +318,8 @@ end subroutine rrtmgp_set_gases_sw !================================================================================================== -subroutine rrtmgp_set_cloud_lw( & - state, pbuf, ncol, nlay, nlaycam, & - cld, cldfsnow, cldfgrau, cldfprime, graupel_in_rad, & - kdist_lw, cloud_lw, cld_lw_abs_cloudsim, snow_lw_abs_cloudsim, grau_lw_abs_cloudsim) - - ! Compute combined cloud optical properties. - ! Create MCICA stochastic arrays for cloud LW optical properties. - ! Initialize optical properties object (cloud_lw) and load with MCICA columns. - - ! arguments - type(physics_state), intent(in) :: state - type(physics_buffer_desc), pointer :: pbuf(:) - integer, intent(in) :: ncol ! number of columns in CAM chunk - integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer") - integer, intent(in) :: nlaycam ! number of CAM layers in radiation calculation - real(r8), pointer :: cld(:,:) ! cloud fraction (liq+ice) - real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds" - real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds" - real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction - - logical, intent(in) :: graupel_in_rad ! use graupel in radiation code - class(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw - type(ty_optical_props_1scl), intent(out) :: cloud_lw - - ! Diagnostic outputs - real(r8), intent(out) :: cld_lw_abs_cloudsim(pcols,pver) ! in-cloud liq+ice optical depth (for COSP) - real(r8), intent(out) :: snow_lw_abs_cloudsim(pcols,pver)! in-cloud snow optical depth (for COSP) - real(r8), intent(out) :: grau_lw_abs_cloudsim(pcols,pver)! in-cloud Graupel optical depth (for COSP) - - ! Local variables - - integer :: i, k - - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: liq_lw_abs(nlwbands,pcols,pver) ! liquid absorption optics depth (LW) - real(r8) :: ice_lw_abs(nlwbands,pcols,pver) ! ice absorption optics depth (LW) - real(r8) :: cld_lw_abs(nlwbands,pcols,pver) ! cloud absorption optics depth (LW) - real(r8) :: snow_lw_abs(nlwbands,pcols,pver) ! snow absorption optics depth (LW) - real(r8) :: grau_lw_abs(nlwbands,pcols,pver) ! graupel absorption optics depth (LW) - real(r8) :: c_cld_lw_abs(nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW) - - ! Arrays for converting from CAM chunks to RRTMGP inputs. - real(r8) :: cldf(ncol,nlaycam) - real(r8) :: tauc(nlwbands,ncol,nlaycam) - real(r8) :: taucmcl(nlwgpts,ncol,nlaycam) - - character(len=128) :: errmsg - character(len=*), parameter :: sub = 'rrtmgp_set_cloud_lw' - !-------------------------------------------------------------------------------- - - ! Combine the cloud optical properties. These calculations are done on CAM "chunks". - - ! gammadist liquid optics - call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs) - ! Mitchell ice optics - call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs) - - cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:) - - if (associated(cldfsnow)) then - ! add in snow - call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs) - do i = 1, ncol - do k = 1, pver - if (cldfprime(i,k) > 0._r8) then - c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) & - + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k) - else - c_cld_lw_abs(:,i,k) = 0._r8 - end if - end do - end do - else - c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:) - end if - - ! add in graupel - if (associated(cldfgrau) .and. graupel_in_rad) then - call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs) - do i = 1, ncol - do k = 1, pver - if (cldfprime(i,k) > 0._r8) then - c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) & - + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k) - else - c_cld_lw_abs(:,i,k) = 0._r8 - end if - end do - end do - end if - - ! Cloud optics for COSP - cld_lw_abs_cloudsim = cld_lw_abs(idx_lw_cloudsim,:,:) - snow_lw_abs_cloudsim = snow_lw_abs(idx_lw_cloudsim,:,:) - grau_lw_abs_cloudsim = grau_lw_abs(idx_lw_cloudsim,:,:) - - ! Extract just the layers of CAM where RRTMGP does calculations. - - ! Subset "chunk" data so just the number of CAM layers in the - ! radiation calculation are used by MCICA to produce subcolumns. - cldf = cldfprime(:ncol, ktopcam:) - tauc = c_cld_lw_abs(:, :ncol, ktopcam:) - - ! Enforce tauc >= 0. - tauc = merge(tauc, 0.0_r8, tauc > 0.0_r8) - - ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point) - call mcica_subcol_lw( & - kdist_lw, nlwbands, nlwgpts, ncol, nlaycam, & - nlwgpts, state%pmid, cldf, tauc, taucmcl ) - - errmsg =cloud_lw%alloc_1scl(ncol, nlay, kdist_lw) - if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: cloud_lw%alloc_1scalar: '//trim(errmsg)) - end if - - ! If there is an extra layer in the radiation then this initialization - ! will provide zero optical depths there. - cloud_lw%tau = 0.0_r8 - - ! Set the properties on g-points. - do i = 1, nlwgpts - cloud_lw%tau(:,ktoprad:,i) = taucmcl(i,:,:) - end do - - ! validate checks that: tau > 0 - errmsg = cloud_lw%validate() - if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: cloud_lw%validate: '//trim(errmsg)) - end if - -end subroutine rrtmgp_set_cloud_lw - -!================================================================================================== - subroutine rrtmgp_set_cloud_sw( & - state, pbuf, nlay, nday, idxday, & + state, pbuf, nlay, nday, idxday, nswgpts, & nnite, idxnite, pmid, cld, cldfsnow, & cldfgrau, cldfprime, graupel_in_rad, kdist_sw, cloud_sw, & tot_cld_vistau, tot_icld_vistau, liq_icld_vistau, ice_icld_vistau, snow_icld_vistau, & @@ -625,6 +335,7 @@ subroutine rrtmgp_set_cloud_sw( & integer, intent(in) :: nlay ! number of layers in radiation calculation (may include "extra layer") integer, intent(in) :: nday ! number of daylight columns integer, intent(in) :: idxday(pcols) ! indices of daylight columns in the chunk + integer, intent(in) :: nswgpts integer, intent(in) :: nnite ! number of night columns integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk @@ -635,9 +346,9 @@ subroutine rrtmgp_set_cloud_sw( & real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "graupel clouds" real(r8), intent(in) :: cldfprime(pcols,pver) ! combined cloud fraction - logical, intent(in) :: graupel_in_rad ! graupel in radiation code - class(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw ! shortwave gas optics object - type(ty_optical_props_2str), intent(out) :: cloud_sw ! SW cloud optical properties object + logical, intent(in) :: graupel_in_rad ! graupel in radiation code + class(ty_gas_optics_rrtmgp_ccpp), intent(in) :: kdist_sw ! shortwave gas optics object + type(ty_optical_props_2str_ccpp), intent(out) :: cloud_sw ! SW cloud optical properties object ! Diagnostic outputs real(r8), intent(out) :: tot_cld_vistau(pcols,pver) ! gbx total cloud optical depth @@ -840,39 +551,39 @@ subroutine rrtmgp_set_cloud_sw( & ! MCICA uses spectral data (on bands) to construct subcolumns (one per g-point) call mcica_subcol_sw( & - kdist_sw, nswbands, nswgpts, nday, nlay, & + kdist_sw%gas_props, nswbands, nswgpts, nday, nlay, & nver, changeseed, pmid, cldf, tauc, & ssac, asmc, taucmcl, ssacmcl, asmcmcl) ! Initialize object for SW cloud optical properties. - errmsg = cloud_sw%alloc_2str(nday, nlay, kdist_sw) + errmsg = cloud_sw%optical_props%alloc_2str(nday, nlay, kdist_sw%gas_props) if (len_trim(errmsg) > 0) then - call endrun(trim(sub)//': ERROR: cloud_sw%alloc_2str: '//trim(errmsg)) + call endrun(trim(sub)//': ERROR: cloud_sw%optical_props%alloc_2str: '//trim(errmsg)) end if ! If there is an extra layer in the radiation then this initialization ! will provide the optical properties there. - cloud_sw%tau = 0.0_r8 - cloud_sw%ssa = 1.0_r8 - cloud_sw%g = 0.0_r8 + cloud_sw%optical_props%tau = 0.0_r8 + cloud_sw%optical_props%ssa = 1.0_r8 + cloud_sw%optical_props%g = 0.0_r8 ! Set the properties on g-points. do igpt = 1,nswgpts - cloud_sw%g (:,ktoprad:,igpt) = asmcmcl(igpt,:,:) - cloud_sw%ssa(:,ktoprad:,igpt) = ssacmcl(igpt,:,:) - cloud_sw%tau(:,ktoprad:,igpt) = taucmcl(igpt,:,:) + cloud_sw%optical_props%g (:,ktoprad:,igpt) = asmcmcl(igpt,:,:) + cloud_sw%optical_props%ssa(:,ktoprad:,igpt) = ssacmcl(igpt,:,:) + cloud_sw%optical_props%tau(:,ktoprad:,igpt) = taucmcl(igpt,:,:) end do ! validate checks that: tau > 0, ssa is in range [0,1], and g is in range [-1,1]. - errmsg = cloud_sw%validate() + errmsg = cloud_sw%optical_props%validate() if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: cloud_sw%validate: '//trim(errmsg)) + call endrun(sub//': ERROR: cloud_sw%optical_props%validate: '//trim(errmsg)) end if ! delta scaling adjusts for forward scattering - errmsg = cloud_sw%delta_scale() + errmsg = cloud_sw%optical_props%delta_scale() if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: cloud_sw%delta_scale: '//trim(errmsg)) + call endrun(sub//': ERROR: cloud_sw%optical_props%delta_scale: '//trim(errmsg)) end if ! All information is in cloud_sw, now deallocate local vars. @@ -896,7 +607,7 @@ subroutine rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw) type(physics_state), target, intent(in) :: state type(physics_buffer_desc), pointer :: pbuf(:) - type(ty_optical_props_1scl), intent(inout) :: aer_lw + type(ty_optical_props_1scl_ccpp), intent(inout) :: aer_lw ! Local variables integer :: ncol @@ -915,13 +626,13 @@ subroutine rrtmgp_set_aer_lw(icall, state, pbuf, aer_lw) ! If there is an extra layer in the radiation then this initialization ! will provide zero optical depths there. - aer_lw%tau = 0.0_r8 + aer_lw%optical_props%tau = 0.0_r8 - aer_lw%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :) + aer_lw%optical_props%tau(:ncol, ktoprad:, :) = aer_lw_abs(:ncol, ktopcam:, :) - errmsg = aer_lw%validate() + errmsg = aer_lw%optical_props%validate() if (len_trim(errmsg) > 0) then - call endrun(sub//': ERROR: aer_lw%validate: '//trim(errmsg)) + call endrun(sub//': ERROR: aer_lw%optical_props%validate: '//trim(errmsg)) end if end subroutine rrtmgp_set_aer_lw @@ -941,7 +652,7 @@ subroutine rrtmgp_set_aer_sw( & integer, intent(in) :: nnite ! number of night columns integer, intent(in) :: idxnite(pcols) ! indices of night columns in the chunk - type(ty_optical_props_2str), intent(inout) :: aer_sw + type(ty_optical_props_2str_ccpp), intent(inout) :: aer_sw ! local variables integer :: i @@ -977,31 +688,31 @@ subroutine rrtmgp_set_aer_sw( & ! If there is an extra layer in the radiation then this initialization ! will provide default values. - aer_sw%tau = 0.0_r8 - aer_sw%ssa = 1.0_r8 - aer_sw%g = 0.0_r8 + aer_sw%optical_props%tau = 0.0_r8 + aer_sw%optical_props%ssa = 1.0_r8 + aer_sw%optical_props%g = 0.0_r8 ! CAM fields are products tau, tau*ssa, tau*ssa*asy ! Fields expected by RRTMGP are computed by - ! aer_sw%tau = aer_tau - ! aer_sw%ssa = aer_tau_w / aer_tau - ! aer_sw%g = aer_tau_w_g / aer_taw_w + ! aer_sw%optical_props%tau = aer_tau + ! aer_sw%optical_props%ssa = aer_tau_w / aer_tau + ! aer_sw%optical_props%g = aer_tau_w_g / aer_taw_w ! aer_sw arrays have dimensions of (nday,nlay,nswbands) do i = 1, nday ! set aerosol optical depth, clip to zero - aer_sw%tau(i,ktoprad:,:) = max(aer_tau(idxday(i),ktopcam:,:), 0._r8) + aer_sw%optical_props%tau(i,ktoprad:,:) = max(aer_tau(idxday(i),ktopcam:,:), 0._r8) ! set value of single scattering albedo - aer_sw%ssa(i,ktoprad:,:) = merge(aer_tau_w(idxday(i),ktopcam:,:)/aer_tau(idxday(i),ktopcam:,:), & + aer_sw%optical_props%ssa(i,ktoprad:,:) = merge(aer_tau_w(idxday(i),ktopcam:,:)/aer_tau(idxday(i),ktopcam:,:), & 1._r8, aer_tau(idxday(i),ktopcam:,:) > 0._r8) ! set value of asymmetry - aer_sw%g(i,ktoprad:,:) = merge(aer_tau_w_g(idxday(i),ktopcam:,:)/aer_tau_w(idxday(i),ktopcam:,:), & + aer_sw%optical_props%g(i,ktoprad:,:) = merge(aer_tau_w_g(idxday(i),ktopcam:,:)/aer_tau_w(idxday(i),ktopcam:,:), & 0._r8, aer_tau_w(idxday(i),ktopcam:,:) > tiny) end do ! impose limits on the components - aer_sw%ssa = min(max(aer_sw%ssa, 0._r8), 1._r8) - aer_sw%g = min(max(aer_sw%g, -1._r8), 1._r8) + aer_sw%optical_props%ssa = min(max(aer_sw%optical_props%ssa, 0._r8), 1._r8) + aer_sw%optical_props%g = min(max(aer_sw%optical_props%g, -1._r8), 1._r8) end if @@ -1009,4 +720,4 @@ end subroutine rrtmgp_set_aer_sw !================================================================================================== -end module rrtmgp_inputs +end module rrtmgp_inputs_cam diff --git a/src/utils/cam_ccpp/machine.F90 b/src/utils/cam_ccpp/machine.F90 new file mode 100644 index 0000000000..4d1a37e4ad --- /dev/null +++ b/src/utils/cam_ccpp/machine.F90 @@ -0,0 +1,12 @@ +! This module is the CAM version of the CCPP generated module of the same name +module machine + + use ccpp_kinds, only: kind_phys => kind_phys + + + implicit none + private + + public kind_phys + +end module machine