From 3fec480eb9e957fe6dbc1f122e40a20d5b3fb6d7 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 14 Dec 2023 22:44:10 -0500 Subject: [PATCH] WIP: compiles and runs tests; clear-sky pasess but all-sky not --- examples/all-sky/rrtmgp_allsky.F90 | 14 ++- examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 | 9 +- examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 | 11 +-- extensions/mo_compute_bc.F90 | 4 +- extensions/mo_rrtmgp_clr_all_sky.F90 | 32 ++----- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 31 ++++--- rte-frontend/mo_optical_props.F90 | 36 +++++++- rte-frontend/mo_rte_lw.F90 | 20 ++-- rte-frontend/mo_rte_sw.F90 | 29 +++--- tests/check_equivalence.F90 | 91 ++++++++----------- tests/check_variants.F90 | 77 ++++++++-------- ...test_zenith_angle_spherical_correction.F90 | 5 +- 12 files changed, 176 insertions(+), 183 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index e25ec0bf3..33432f097 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -89,7 +89,7 @@ program rte_rrtmgp_allsky ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: nbnd, ngpt integer :: icol, ilay, ibnd, iloop, igas @@ -191,7 +191,6 @@ program rte_rrtmgp_allsky ! nbnd = k_dist%get_nband() ngpt = k_dist%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! LW calculations neglect scattering; SW calculations use the 2-stream approximation @@ -251,7 +250,7 @@ program rte_rrtmgp_allsky ! Surface temperature !$acc kernels !$omp target - t_sfc = t_lev(1, merge(nlay+1, 1, top_at_1)) + t_sfc = t_lev(1, merge(nlay+1, 1, atmos%top_is_at_1())) emis_sfc = 0.98_wp !$acc end kernels !$omp end target @@ -322,9 +321,9 @@ program rte_rrtmgp_allsky tlev = t_lev)) if(do_clouds) call stop_on_err(clouds%increment(atmos)) if(do_aerosols) call stop_on_err(aerosols%increment(atmos)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - emis_sfc, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + emis_sfc, & fluxes)) !$acc end data !$omp end target data @@ -347,8 +346,7 @@ program rte_rrtmgp_allsky call stop_on_err(aerosols%delta_scale()) call stop_on_err(aerosols%increment(atmos)) end if - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) !$acc end data diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 index 607e372ef..07fac9c26 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 @@ -90,7 +90,6 @@ program rrtmgp_rfmip_lw character(len=132) :: kdist_file = 'coefficients_lw.nc' character(len=132) :: flxdn_file, flxup_file integer :: nargs, ncol, nlay, nbnd, nexp, nblocks, block_size, forcing_index, physics_index, n_quad_angles = 1 - logical :: top_at_1 integer :: b, icol, ibnd character(len=4) :: block_size_char, forcing_index_char = '1', physics_index_char = '1' @@ -168,10 +167,6 @@ program rrtmgp_rfmip_lw ! Allocation on assignment within reading routines ! call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev) - ! - ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? - ! - top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1) ! ! Read the gas concentrations and surface properties @@ -193,7 +188,8 @@ program rrtmgp_rfmip_lw ! is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep. ! This introduces an error but shows input sanitizing. ! - if(top_at_1) then + ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? + if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) else p_lev(:,nlay+1,:) & @@ -256,7 +252,6 @@ program rrtmgp_rfmip_lw ! via ty_fluxes_broadband ! call stop_on_err(rte_lw(optical_props, & - top_at_1, & source, & sfc_emis_spec, & fluxes, n_gauss_angles = n_quad_angles)) diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 index 547d9b4e3..e914f3f93 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 @@ -90,7 +90,6 @@ program rrtmgp_rfmip_sw character(len=132) :: kdist_file = 'coefficients_sw.nc' character(len=132) :: flxdn_file, flxup_file integer :: nargs, ncol, nlay, nbnd, ngpt, nexp, nblocks, block_size, forcing_index - logical :: top_at_1 integer :: b, icol, ibnd, igpt character(len=4) :: block_size_char, forcing_index_char = '1' @@ -166,10 +165,6 @@ program rrtmgp_rfmip_sw ! Allocation on assignment within reading routines ! call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev) - ! - ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? - ! - top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1) ! ! Read the gas concentrations and surface properties @@ -193,7 +188,10 @@ program rrtmgp_rfmip_sw ! is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep. ! This introduces an error but shows input sanitizing. ! - if(top_at_1) then + ! + ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? + ! + if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) else p_lev(:,nlay+1,:) & @@ -296,7 +294,6 @@ program rrtmgp_rfmip_sw ! via ty_fluxes_broadband ! call stop_on_err(rte_sw(optical_props, & - top_at_1, & mu0, & toa_flux, & sfc_alb_spec, & diff --git a/extensions/mo_compute_bc.F90 b/extensions/mo_compute_bc.F90 index ef2a50b29..fbec1d974 100644 --- a/extensions/mo_compute_bc.F90 +++ b/extensions/mo_compute_bc.F90 @@ -166,7 +166,6 @@ function compute_bc(k_dist, & ! Compute fluxes ! error_msg = rte_lw(optical_props_1lay, & - top_at_1, & lw_sources_1lay, & lower_bc, fluxes_1lev) else @@ -191,8 +190,7 @@ function compute_bc(k_dist, & optical_props_1lay, & solar_src) error_msg = rte_sw(optical_props_1lay, & - top_at_1, mu0, & - solar_src, & + mu0, solar_src, & lower_bc, lower_bc, fluxes_1lev) endif end function diff --git a/extensions/mo_rrtmgp_clr_all_sky.F90 b/extensions/mo_rrtmgp_clr_all_sky.F90 index 755681cbb..47ab1e0d2 100644 --- a/extensions/mo_rrtmgp_clr_all_sky.F90 +++ b/extensions/mo_rrtmgp_clr_all_sky.F90 @@ -85,11 +85,6 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ngpt = k_dist%get_ngpt() nband = k_dist%get_nband() - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) - !$acc end kernels - !$omp end target ! ------------------------------------------------------------------------------------ ! Error checking @@ -161,8 +156,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & if(present(aer_props)) error_msg = aer_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_lw(optical_props, top_at_1, sources, & - sfc_emis, clrsky_fluxes, & + error_msg = base_rte_lw(optical_props, sources, & + sfc_emis, clrsky_fluxes, & inc_flux, n_gauss_angles) if(error_msg /= '') return ! ------------------------------------------------------------------------------------ @@ -171,8 +166,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & error_msg = cloud_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_lw(optical_props, top_at_1, sources, & - sfc_emis, allsky_fluxes, & + error_msg = base_rte_lw(optical_props, sources, & + sfc_emis, allsky_fluxes, & inc_flux, n_gauss_angles) call sources%finalize() @@ -207,7 +202,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & class(ty_optical_props_arry), allocatable :: optical_props real(wp), dimension(:,:), allocatable :: toa_flux integer :: ncol, nlay, ngpt, nband, nstr - logical :: top_at_1 ! -------------------------------- ! Problem sizes ! @@ -218,12 +212,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ngpt = k_dist%get_ngpt() nband = k_dist%get_nband() - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) - !$acc end kernels - !$omp end target - ! ------------------------------------------------------------------------------------ ! Error checking ! @@ -276,7 +264,7 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ! Gas optical depth -- pressure need to be expressed as Pa ! error_msg = k_dist%gas_optics(p_lay, p_lev, t_lay, gas_concs, & - optical_props, toa_flux, & + optical_props, toa_flux, & col_dry) if (error_msg /= '') return ! @@ -289,9 +277,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & if(present(aer_props)) error_msg = aer_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_sw(optical_props, top_at_1, & - mu0, toa_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = base_rte_sw(optical_props, mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & clrsky_fluxes) if(error_msg /= '') return @@ -301,9 +288,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & error_msg = cloud_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_sw(optical_props, top_at_1, & - mu0, toa_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = base_rte_sw(optical_props, mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & allsky_fluxes) call optical_props%finalize() diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 7583fe579..e169f850d 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -252,12 +252,17 @@ function gas_optics_int(this, & integer, dimension(2, size(play,dim=1),size(play,dim=2), get_nflav(this)) :: jeta integer :: ncol, nlay, ngpt, nband + logical :: temp ! ---------------------------------------------------------- ncol = size(play,dim=1) nlay = size(play,dim=2) ngpt = this%get_ngpt() nband = this%get_nband() ! + ! Vertical orientation + ! + call optical_props%set_top_at_1(play(1,1) < play(1, nlay)) + ! ! Gas optics ! !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta) @@ -311,7 +316,7 @@ function gas_optics_int(this, & ! but isn't with PGI 19.10 ! error_msg = source(this, & - ncol, nlay, nband, ngpt, & + ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources, & @@ -320,7 +325,7 @@ function gas_optics_int(this, & !$omp target exit data map(release:tlev) else error_msg = source(this, & - ncol, nlay, nband, ngpt, & + ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources) @@ -329,6 +334,7 @@ function gas_optics_int(this, & !$omp target exit data map(release:tsfc) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta) + end function gas_optics_int !------------------------------------------------------------------------------------------ ! @@ -372,6 +378,10 @@ function gas_optics_ext(this, & ngas = this%get_ngas() nflav = get_nflav(this) ! + ! Vertical orientation + ! + call optical_props%set_top_at_1(play(1,1) < play(1, nlay)) + ! ! Gas optics ! !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta) @@ -384,6 +394,7 @@ function gas_optics_ext(this, & col_dry) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta) + if(error_msg /= '') return ! ---------------------------------------------------------- @@ -824,7 +835,7 @@ end function set_tsi ! Compute Planck source functions at layer centers and levels ! function source(this, & - ncol, nlay, nbnd, ngpt, & + ncol, nlay, nbnd, ngpt, top_at_1, & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources, & ! Planck sources @@ -833,6 +844,7 @@ function source(this, & ! inputs class(ty_gas_optics_rrtmgp), intent(in ) :: this integer, intent(in ) :: ncol, nlay, nbnd, ngpt + logical, intent(in ) :: top_at_1 real(wp), dimension(ncol,nlay), intent(in ) :: play ! layer pressures [Pa, mb] real(wp), dimension(ncol,nlay+1), intent(in ) :: plev ! level pressures [Pa, mb] real(wp), dimension(ncol,nlay), intent(in ) :: tlay ! layer temperatures [K] @@ -849,7 +861,6 @@ function source(this, & optional, target :: tlev ! level temperatures [K] character(len=128) :: error_msg ! ---------------------------------------------------------- - logical(wl) :: top_at_1 integer :: icol, ilay ! Variables for temperature at layer edges [K] (ncol, nlay+1) real(wp), dimension( ncol,nlay+1), target :: tlev_arr @@ -895,18 +906,12 @@ function source(this, & !$omp target data map(from:sources%lay_source, sources%lev_source) & !$omp map(from:sources%sfc_source, sources%sfc_source_Jac) - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = play(1,1) < play(1, nlay) - !$acc end kernels - !$omp end target - call compute_Planck_source(ncol, nlay, nbnd, ngpt, & get_nflav(this), this%get_neta(), this%get_npres(), this%get_ntemp(), this%get_nPlanckTemp(), & - tlay, tlev_wk, tsfc, merge(nlay, 1, top_at_1), & - fmajor, jeta, tropo, jtemp, jpress, & + tlay, tlev_wk, tsfc, merge(nlay, 1, logical(top_at_1, wl)), & + fmajor, jeta, tropo, jtemp, jpress, & this%get_gpoint_bands(), this%get_band_lims_gpoint(), this%planck_frac, this%temp_ref_min,& - this%totplnk_delta, this%totplnk, this%gpoint_flavor, & + this%totplnk_delta, this%totplnk, this%gpoint_flavor, & sources%sfc_source, sources%lay_source, sources%lev_source, & sources%sfc_source_Jac) !$acc end data diff --git a/rte-frontend/mo_optical_props.F90 b/rte-frontend/mo_optical_props.F90 index c49fd5dc8..83627bb23 100644 --- a/rte-frontend/mo_optical_props.F90 +++ b/rte-frontend/mo_optical_props.F90 @@ -33,6 +33,7 @@ !> !> These classes must be allocated before use. Initialization and allocation can be combined. !> The classes have a validate() function that checks all arrays for valid values (e.g. tau > 0.) +!> The vertical orientation can be specified via this%set_top_at_1() or obtained via this%top_at_1(). !> !> Optical properties can be delta-scaled (though this is currently implemented only for two-stream arrays) !> @@ -51,7 +52,7 @@ !@endnote !> ------------------------------------------------------------------------------------------------- module mo_optical_props - use mo_rte_kind, only: wp + use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are @@ -109,6 +110,7 @@ module mo_optical_props ! ------------------------------------------------------------------------------------------------- type, extends(ty_optical_props), abstract, public :: ty_optical_props_arry real(wp), dimension(:,:,:), allocatable :: tau !! optical depth (ncol, nlay, ngpt) + logical, private :: top_at_1 ! No default - maybe uninitialized values will get caught? contains procedure, public :: get_ncol procedure, public :: get_nlay @@ -116,6 +118,11 @@ module mo_optical_props !> Increment another set of values !> procedure, public :: increment + !> + !> + !> + procedure, public :: top_is_at_1 + procedure, public :: set_top_at_1 !> !> Deferred procedures -- each must be implemented in each child class with @@ -757,6 +764,7 @@ function subset_1scl_range(full, start, n, subset) result(err_message) subset%p(:,1:n,:,:) = 0._wp end select call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau) + call subset%set_top_at_1(full%top_is_at_1()) end function subset_1scl_range ! ------------------------------------------------------------------------------------------ @@ -810,6 +818,7 @@ function subset_2str_range(full, start, n, subset) result(err_message) subset%p(1,1:n,:,:) = full%g (start:start+n-1,:,:) subset%p(2:,:, :,:) = 0._wp end select + call subset%set_top_at_1(full%top_is_at_1()) end function subset_2str_range ! ------------------------------------------------------------------------------------------ @@ -858,6 +867,7 @@ function subset_nstr_range(full, start, n, subset) result(err_message) call extract_subset( ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa) call extract_subset(nmom, ncol, nlay, ngpt, full%p , start, start+n-1, subset%p ) end select + call subset%set_top_at_1(full%top_is_at_1()) end function subset_nstr_range !> ------------------------------------------------------------------------------------------ @@ -1057,6 +1067,30 @@ pure function get_nmom(this) get_nmom = 0 end if end function get_nmom + !> ----------------------------------------------------------------------------------------------- + !> + !> Routines for array classes: vertical orientation + !> + ! ------------------------------------------------------------------------------------------ + pure function top_is_at_1(this) + ! + ! Vertical orientation - .true. if array index 1 is top of atmosphere + ! + class(ty_optical_props_arry), intent(in) :: this + logical :: top_is_at_1 + + top_is_at_1 = this%top_at_1 + end function top_is_at_1 + ! ------------------------------------------------------------------------------------------ + subroutine set_top_at_1(this, top_at_1) + ! + !> Set vertical orientation of class - .true. if array index 1 is top of atmosphere + ! + class(ty_optical_props_arry), intent(inout) :: this + logical, intent(in ) :: top_at_1 + + this%top_at_1 = top_at_1 + end subroutine set_top_at_1 ! ----------------------------------------------------------------------------------------------- ! ! Routines for base class: spectral discretization diff --git a/rte-frontend/mo_rte_lw.F90 b/rte-frontend/mo_rte_lw.F90 index c36c55f8c..85b4a5b47 100644 --- a/rte-frontend/mo_rte_lw.F90 +++ b/rte-frontend/mo_rte_lw.F90 @@ -67,15 +67,13 @@ module mo_rte_lw ! Interface using only optical properties and source functions as inputs; fluxes as outputs. ! ! -------------------------------------------------- - function rte_lw(optical_props, top_at_1, & - sources, sfc_emis, & - fluxes, & + function rte_lw(optical_props, & + sources, sfc_emis, & + fluxes, & inc_flux, n_gauss_angles, use_2stream, & lw_Ds, flux_up_Jac) result(error_msg) class(ty_optical_props_arry), intent(in ) :: optical_props !! Set of optical properties as one or more arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) type(ty_source_func_lw), intent(in ) :: sources !! Derived type with Planck source functions real(wp), dimension(:,:), intent(in ) :: sfc_emis @@ -350,8 +348,8 @@ function rte_lw(optical_props, top_at_1, & end do end if call lw_solver_noscat(ncol, nlay, ngpt, & - logical(top_at_1, wl), n_quad_angs, & - secants, gauss_wts(1:n_quad_angs,n_quad_angs), & + logical(optical_props%top_is_at_1(), wl), & + n_quad_angs, secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & sources%lay_source, & sources%lev_source, & @@ -370,7 +368,8 @@ function rte_lw(optical_props, top_at_1, & ! ! two-stream calculation with scattering ! - call lw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & + call lw_solver_2stream(ncol, nlay, ngpt, & + logical(optical_props%top_is_at_1(), wl), & optical_props%tau, optical_props%ssa, optical_props%g, & sources%lay_source, sources%lev_source, & sfc_emis_gpt, sources%sfc_source, & @@ -393,7 +392,8 @@ function rte_lw(optical_props, top_at_1, & ! Re-scaled solution to account for scattering ! call lw_solver_noscat(ncol, nlay, ngpt, & - logical(top_at_1, wl), n_quad_angs, & + logical(optical_props%top_is_at_1(), wl), & + n_quad_angs, & secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & sources%lay_source, & @@ -435,7 +435,7 @@ function rte_lw(optical_props, top_at_1, & ! ! ...or reduce spectral fluxes to desired output quantities ! - error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, top_at_1) + error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, optical_props%top_is_at_1()) end select end if ! no error message from validation !$acc end data diff --git a/rte-frontend/mo_rte_sw.F90 b/rte-frontend/mo_rte_sw.F90 index 961209543..144da0514 100644 --- a/rte-frontend/mo_rte_sw.F90 +++ b/rte-frontend/mo_rte_sw.F90 @@ -53,14 +53,12 @@ module mo_rte_sw contains ! ------------------------------------------------------------------------------------------------- - function rte_sw_mu0_bycol(atmos, top_at_1, & + function rte_sw_mu0_bycol(atmos, & mu0, inc_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes, inc_flux_dif) result(error_msg) class(ty_optical_props_arry), intent(in ) :: atmos !! Optical properties provided as arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) real(wp), dimension(:), intent(in ) :: mu0 !! cosine of solar zenith angle (ncol) - will be assumed constant with height real(wp), dimension(:,:), intent(in ) :: inc_flux @@ -94,22 +92,19 @@ function rte_sw_mu0_bycol(atmos, top_at_1, & end do end do - error_msg = rte_sw_mu0_full(atmos, top_at_1, & - mu0_bylay, inc_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = rte_sw_mu0_full(atmos, & + mu0_bylay, inc_flux, & + sfc_alb_dir, sfc_alb_dif, & fluxes, inc_flux_dif) !$acc end data !$omp end target data end function rte_sw_mu0_bycol ! ------------------------------------------------------------------------------------------------- - function rte_sw_mu0_full(atmos, top_at_1, & - mu0, inc_flux, & - sfc_alb_dir, sfc_alb_dif, & - fluxes, inc_flux_dif) result(error_msg) + function rte_sw_mu0_full(atmos, mu0, inc_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes, inc_flux_dif) result(error_msg) class(ty_optical_props_arry), intent(in ) :: atmos !! Optical properties provided as arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) real(wp), dimension(:,:), intent(in ) :: mu0 !! cosine of solar zenith angle (ncol, nlay) real(wp), dimension(:,:), intent(in ) :: inc_flux @@ -293,8 +288,9 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! Direct beam only - for completeness, unlikely to be used in practice ! - call sw_solver_noscat(ncol, nlay, ngpt, logical(top_at_1, wl), & - atmos%tau, mu0, inc_flux, & + call sw_solver_noscat(ncol, nlay, ngpt, & + logical(atmos%top_is_at_1(), wl), & + atmos%tau, mu0, inc_flux, & gpt_flux_dir) call zero_array(ncol, nlay+1, ngpt, gpt_flux_up) ! @@ -308,7 +304,8 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! two-stream calculation with scattering ! - call sw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & + call sw_solver_2stream(ncol, nlay, ngpt, & + logical(atmos%top_is_at_1(), wl), & atmos%tau, atmos%ssa, atmos%g, mu0, & sfc_alb_dir_gpt, sfc_alb_dif_gpt, & inc_flux, & @@ -348,7 +345,7 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! ...or reduce spectral fluxes to desired output quantities ! - error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, atmos, top_at_1, gpt_flux_dir) + error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, atmos, atmos%top_is_at_1(), gpt_flux_dir) end select end if ! In case of an error we exit here diff --git a/tests/check_equivalence.F90 b/tests/check_equivalence.F90 index 8ac814965..9b057f7bc 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -100,7 +100,7 @@ program rte_check_equivalence ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: ncol, nlay, nbnd, ngpt, nexp integer :: icol, ilay, ibnd, iloop, igas @@ -169,7 +169,6 @@ program rte_check_equivalence ! nbnd = k_dist%get_nband() ngpt = k_dist%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! ! Boundary conditions @@ -237,9 +236,9 @@ program rte_check_equivalence atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) print *, " Default calculation" ! @@ -254,17 +253,17 @@ program rte_check_equivalence nullify(fluxes%flux_up) nullify(fluxes%flux_dn) allocate(fluxes%flux_net(ncol,nlay+1)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call stop_on_err("Net fluxes don't match when computed alone") fluxes%flux_up => tst_flux_up(:,:) fluxes%flux_dn => tst_flux_dn(:,:) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call report_err("Net fluxes don't match when computed in tandem") @@ -294,36 +293,36 @@ program rte_check_equivalence ! atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" halving/doubling fails") call increment_with_1scl - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 1scl fails") call increment_with_2str - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 2str fails") call increment_with_nstr - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & @@ -333,10 +332,10 @@ program rte_check_equivalence ! ! Computing Jacobian shouldn't change net fluxes ! - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & - fluxes, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & + fluxes, & flux_up_Jac = jFluxUp)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & @@ -350,9 +349,9 @@ program rte_check_equivalence atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) ! ! Comparision of fluxes with increased surface T aren't expected to match @@ -386,8 +385,7 @@ program rte_check_equivalence gas_concs, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) print *, " Default calculation" @@ -428,8 +426,7 @@ program rte_check_equivalence toa_flux)) atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & @@ -438,8 +435,7 @@ program rte_check_equivalence call report_err(" halving/doubling fails") call increment_with_1scl - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & @@ -464,8 +460,7 @@ program rte_check_equivalence atmos, & toa_flux)) call increment_with_nstr - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & @@ -496,7 +491,6 @@ subroutine lw_clear_sky_vr t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) t_lev (:,:) = t_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 ! ! No direct access to gas concentrations so use the classes ! This also tests otherwise uncovered routines for ty_gas_concs @@ -515,9 +509,9 @@ subroutine lw_clear_sky_vr atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) tst_flux_up(:,:) = tst_flux_up(:,(nlay+1):1:-1) tst_flux_dn(:,:) = tst_flux_dn(:,(nlay+1):1:-1) @@ -525,7 +519,6 @@ subroutine lw_clear_sky_vr t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) t_lev (:,:) = t_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 end subroutine lw_clear_sky_vr ! ---------------------------------------------------------------------------- ! @@ -555,8 +548,8 @@ subroutine lw_clear_sky_subset colE = i * ncol/2 call stop_on_err(atmos%get_subset (colS, ncol/2, atmos_subset)) call stop_on_err(lw_sources%get_subset(colS, ncol/2, sources_subset)) - call stop_on_err(rte_lw(atmos_subset, top_at_1, & - sources_subset, & + call stop_on_err(rte_lw(atmos_subset, & + sources_subset, & sfc_emis(:, colS:colE), & fluxes)) tst_flux_up(colS:colE,:) = up @@ -583,7 +576,6 @@ subroutine sw_clear_sky_vr p_lay (:,:) = p_lay (:, nlay :1:-1) t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 ! ! No direct access to gas concentrations so use the classes ! This also tests otherwise uncovered routines for ty_gas_concs @@ -601,8 +593,7 @@ subroutine sw_clear_sky_vr gas_concs_vr, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! @@ -615,7 +606,6 @@ subroutine sw_clear_sky_vr p_lay (:,:) = p_lay (:, nlay :1:-1) t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 end subroutine sw_clear_sky_vr ! ---------------------------------------------------------------------------- ! @@ -634,8 +624,7 @@ subroutine sw_clear_sky_tsi gas_concs, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) tst_flux_up (:,:) = tst_flux_up (:,:) / tsi_scale diff --git a/tests/check_variants.F90 b/tests/check_variants.F90 index 72564c2d1..12ed082ef 100644 --- a/tests/check_variants.F90 +++ b/tests/check_variants.F90 @@ -102,7 +102,7 @@ program rte_clear_sky_regression ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: ncol, nlay, nbnd, ngpt, nexp integer :: icol, ilay, ibnd, iloop, igas @@ -168,7 +168,6 @@ program rte_clear_sky_regression ! nbnd = k_dist%get_nband() ngpt = k_dist%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! ! Boundary conditions @@ -260,9 +259,9 @@ subroutine lw_clear_sky_default atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_up, "lw_flux_up", "LW flux up") call write_broadband_field(input_file, flux_dn, "lw_flux_dn", "LW flux dn") @@ -276,9 +275,9 @@ subroutine lw_clear_sky_default ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_net, "lw_flux_net_2", "LW flux net, direct") fluxes%flux_up => flux_up @@ -296,9 +295,9 @@ subroutine lw_clear_sky_notlev gas_concs, & atmos, & lw_sources)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_up, "lw_flux_up_notlev", "LW flux up, no level temperatures") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_notlev", "LW flux dn, no level temperatures") @@ -314,9 +313,9 @@ subroutine lw_clear_sky_3ang atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, n_gauss_angles=3)) call write_broadband_field(input_file, flux_up, "lw_flux_up_3ang", "LW flux up, three quadrature angles") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_3ang", "LW flux dn, three quadrature angles") @@ -334,9 +333,9 @@ subroutine lw_clear_sky_optangle lw_sources, & tlev = t_lev)) call stop_on_err(k_dist%compute_optimal_angles(atmos, lw_Ds)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, lw_Ds=lw_Ds)) call write_broadband_field(input_file, flux_up, "lw_flux_up_optang", "LW flux up, single optimal angles") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_optang", "LW flux dn, single optimal angles") @@ -354,10 +353,10 @@ subroutine lw_clear_sky_jaco atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & - fluxes, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & + fluxes, & flux_up_Jac = jFluxUp)) call write_broadband_field(input_file, flux_up, "lw_flux_up_jaco", "LW flux up, computing Jaobians") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_jaco", "LW flux dn, computing Jaobians") @@ -369,9 +368,9 @@ subroutine lw_clear_sky_jaco atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_up, "lw_flux_up_stp1", "LW flux up, surface T+1K") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_stp1", "LW flux dn, surface T+1K") @@ -390,16 +389,16 @@ subroutine lw_clear_sky_2str atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_up, "lw_flux_up_1rescl", "LW flux up, clear-sky _1rescl") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_1rescl", "LW flux dn, clear-sky _1rescl") - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, use_2stream=.true.)) call write_broadband_field(input_file, flux_up, "lw_flux_up_2str", "LW flux up, clear-sky _2str") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_2str", "LW flux dn, clear-sky _2str") @@ -418,9 +417,9 @@ subroutine lw_clear_sky_alt atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) call write_broadband_field(input_file, flux_up, "lw_flux_up_alt", "LW flux up, fewer g-points") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_alt", "LW flux dn, fewer g-points") @@ -431,9 +430,9 @@ subroutine lw_clear_sky_alt vert_dim_name = "layer") call stop_on_err(k_dist_2%compute_optimal_angles(atmos, lw_Ds)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, lw_Ds=lw_Ds)) call write_broadband_field(input_file, flux_up, "lw_flux_up_alt_oa", "LW flux up, fewer g-points, opt. angle") call write_broadband_field(input_file, flux_dn, "lw_flux_dn_alt_oa", "LW flux dn, fewer g-points, opt. angle") @@ -469,8 +468,7 @@ subroutine sw_clear_sky_default rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=ngpt) toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! @@ -500,8 +498,7 @@ subroutine sw_clear_sky_alt rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=k_dist_2%get_ngpt()) toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! diff --git a/tests/test_zenith_angle_spherical_correction.F90 b/tests/test_zenith_angle_spherical_correction.F90 index 851281521..aed77e1ab 100644 --- a/tests/test_zenith_angle_spherical_correction.F90 +++ b/tests/test_zenith_angle_spherical_correction.F90 @@ -40,7 +40,6 @@ program test_solar_zenith_angle p_lev(ncol, nlay+1) real(wp), dimension(ncol, nlay+1), target & :: flux_up, flux_dn, flux_dir - logical :: top_at_1 character(len=128) :: k_dist_file = "rrtmgp-gas-sw-g112.nc" real(wp), dimension(:,:), allocatable & :: toa_flux, sfc_alb_dir, sfc_alb_dif @@ -97,9 +96,7 @@ program test_solar_zenith_angle ! Set small solar zenith angles, varying by column; compute fluxes and heating rates ! call stop_on_err(zenith_angle_with_height(z_lay(:,1), 3 * [0.01_wp, 0.02_wp, 0.03_wp, 0.04_wp], z_lay, mu0)) - top_at_1 = p_lay(1, 1) < p_lay(1,nlay) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) call stop_on_err(compute_heating_rate(flux_up, flux_dn, flux_dir, p_lev, mu0, heating_rate))