From c45b8dc826af94942efee82805a1a2b0c1161607 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 16:47:01 -0400 Subject: [PATCH 01/38] Remove GPTL, reorganize --- examples/all-sky/rrtmgp_allsky.F90 | 92 ++++++++++---------- examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 | 41 --------- examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 | 38 -------- 3 files changed, 46 insertions(+), 125 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 65b983003..423661271 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -33,50 +33,6 @@ subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) !$acc end data !$omp end target data end subroutine vmr_2d_to_1d -! -------------------------------------------------------------------------------------- -! -! Calculate layer relative humidity for aerosol optics calculations -! -subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) - use mo_rte_kind, only: wp - use mo_gas_optics_constants, only: m_h2o, m_dry - - integer, intent(in) :: ncol, nlay - real(wp), intent(in) :: p_lay(ncol,nlay) ! layer pressure (Pa) - real(wp), intent(in) :: t_lay(ncol,nlay) ! layer temperature (K) - real(wp), intent(in) :: vmr_h2o(ncol,nlay) ! water volume mixing ratio - - real(wp), intent(inout) :: relhum(ncol,nlay) ! relative humidity (fraction, 0-1) - - ! Local variables - integer :: i, k - - real(wp) :: mmr_h2o ! water mass mixing ratio - real(wp) :: q_lay ! water specific humidity - real(wp) :: q_lay_min, q_tmp, es_tmp - real(wp) :: mwd, t_ref, rh - - ! Set constants - mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights - t_ref = 273.16_wp ! reference temperature (K) - q_lay_min = 1.e-7_wp ! minimum water mass mixing ratio - ! ------------------- - - ! Derive layer virtual temperature - do i = 1, ncol - do k = 1, nlay - ! Convert h2o vmr to mmr - mmr_h2o = vmr_h2o(i,k) * mwd - q_lay = mmr_h2o / (1 + mmr_h2o) - q_tmp = max(q_lay_min, q_lay) - es_tmp = exp( (17.67_wp * (t_lay(i,k)-t_ref)) / (t_lay(i,k)-29.65_wp) ) - rh = (0.263_wp * p_lay(i,k) * q_tmp) / es_tmp - ! Convert rh from percent to fraction - relhum(i,k) = 0.01_wp * rh - enddo - enddo - -end subroutine get_relhum ! ---------------------------------------------------------------------------------- program rte_rrtmgp_clouds_aerosols use mo_rte_kind, only: wp, i8, wl @@ -182,7 +138,7 @@ program rte_rrtmgp_clouds_aerosols real(wp) :: rel_val, rei_val character(len=8) :: char_input - integer :: nUserArgs=0, nloops + integer :: nUserArgs=0, nloops=1, ncol = 1, nlay=1 logical :: use_luts = .true., write_fluxes = .true. integer, parameter :: ngas = 8 character(len=3), dimension(ngas) & @@ -307,7 +263,7 @@ program rte_rrtmgp_clouds_aerosols end if ! Clouds optical props are defined by band call stop_on_err(clouds%init(k_dist%get_band_lims_wavenumber())) - ! Clouds optical props are defined by band + ! Aerosol optical props are defined by band call stop_on_err(aerosols%init(k_dist%get_band_lims_wavenumber())) ! ! Allocate arrays for the optical properties themselves. @@ -590,4 +546,48 @@ program rte_rrtmgp_clouds_aerosols end if !$acc enter data create(lwp, iwp, rel, rei) !$omp target enter data map(alloc:lwp, iwp, rel, rei) +contains +! -------------------------------------------------------------------------------------- +! +! Calculate layer relative humidity for aerosol optics calculations +! +subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) + use mo_rte_kind, only: wp + use mo_gas_optics_constants, only: m_h2o, m_dry + + integer, intent(in) :: ncol, nlay + real(wp), intent(in) :: p_lay(ncol,nlay) ! layer pressure (Pa) + real(wp), intent(in) :: t_lay(ncol,nlay) ! layer temperature (K) + real(wp), intent(in) :: vmr_h2o(ncol,nlay) ! water volume mixing ratio + + real(wp), intent(inout) :: relhum(ncol,nlay) ! relative humidity (fraction, 0-1) + + ! Local variables + integer :: i, k + + real(wp) :: mmr_h2o ! water mass mixing ratio + real(wp) :: q_lay ! water specific humidity + real(wp) :: q_lay_min, q_tmp, es_tmp + real(wp) :: mwd, t_ref, rh + + ! Set constants + mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights + t_ref = 273.16_wp ! reference temperature (K) + q_lay_min = 1.e-7_wp ! minimum water mass mixing ratio + ! ------------------- + + ! Derive layer virtual temperature + do i = 1, ncol + do k = 1, nlay + ! Convert h2o vmr to mmr + mmr_h2o = vmr_h2o(i,k) * mwd + q_lay = mmr_h2o / (1 + mmr_h2o) + q_tmp = max(q_lay_min, q_lay) + es_tmp = exp( (17.67_wp * (t_lay(i,k)-t_ref)) / (t_lay(i,k)-29.65_wp) ) + rh = (0.263_wp * p_lay(i,k) * q_tmp) / es_tmp + ! Convert rh from percent to fraction + relhum(i,k) = 0.01_wp * rh + enddo + enddo + end subroutine get_relhum end program rte_rrtmgp_clouds_aerosols diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 index b41fd82fc..352c9d728 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 @@ -81,13 +81,6 @@ program rrtmgp_rfmip_lw use mo_load_coefficients, only: load_and_init use mo_rfmip_io, only: read_size, read_and_block_pt, read_and_block_gases_ty, unblock_and_write, & read_and_block_lw_bc, determine_gas_names -#ifdef USE_TIMING - ! - ! Timing library - ! - use gptl, only: gptlstart, gptlstop, gptlinitialize, gptlpr, gptlfinalize, gptlsetoption, & - gptlpercent, gptloverhead -#endif implicit none ! -------------------------------------------------- ! @@ -121,9 +114,6 @@ program rrtmgp_rfmip_lw ! type(ty_gas_concs), dimension(:), allocatable :: gas_conc_array -#ifdef USE_TIMING - integer :: ret, i -#endif ! ------------------------------------------------------------------------------------------------- ! ! Code starts @@ -232,20 +222,9 @@ program rrtmgp_rfmip_lw !$acc enter data create(source, source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source) !$omp target enter data map(alloc:source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source) ! -------------------------------------------------- -#ifdef USE_TIMING - ! - ! Initialize timers - ! - ret = gptlsetoption (gptlpercent, 1) ! Turn on "% of" print - ret = gptlsetoption (gptloverhead, 0) ! Turn off overhead estimate - ret = gptlinitialize() -#endif ! ! Loop over blocks ! -#ifdef USE_TIMING - do i = 1, 4 -#endif do b = 1, nblocks fluxes%flux_up => flux_up(:,:,b) fluxes%flux_dn => flux_dn(:,:,b) @@ -264,9 +243,6 @@ program rrtmgp_rfmip_lw ! Compute the optical properties of the atmosphere and the Planck source functions ! from pressures, temperatures, and gas concentrations... ! -#ifdef USE_TIMING - ret = gptlstart('gas_optics (LW)') -#endif call stop_on_err(k_dist%gas_optics(p_lay(:,:,b), & p_lev(:,:,b), & t_lay(:,:,b), & @@ -275,33 +251,16 @@ program rrtmgp_rfmip_lw optical_props, & source, & tlev = t_lev(:,:,b))) -#ifdef USE_TIMING - ret = gptlstop('gas_optics (LW)') -#endif ! ! ... and compute the spectrally-resolved fluxes, providing reduced values ! via ty_fluxes_broadband ! -#ifdef USE_TIMING - ret = gptlstart('rte_lw') -#endif call stop_on_err(rte_lw(optical_props, & top_at_1, & source, & sfc_emis_spec, & fluxes, n_gauss_angles = n_quad_angles)) -#ifdef USE_TIMING - ret = gptlstop('rte_lw') -#endif end do -#ifdef USE_TIMING - end do - ! - ! End timers - ! - ret = gptlpr(block_size) - ret = gptlfinalize() -#endif !$acc exit data delete(sfc_emis_spec) !$omp target exit data map(release:sfc_emis_spec) !$acc exit data delete(optical_props%tau, optical_props) diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 index 52215ade3..e56298962 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 @@ -81,13 +81,6 @@ program rrtmgp_rfmip_sw use mo_load_coefficients, only: load_and_init use mo_rfmip_io, only: read_size, read_and_block_pt, read_and_block_gases_ty, unblock_and_write, & read_and_block_sw_bc, determine_gas_names -#ifdef USE_TIMING - ! - ! Timing library - ! - use gptl, only: gptlstart, gptlstop, gptlinitialize, gptlpr, gptlfinalize, gptlsetoption, & - gptlpercent, gptloverhead -#endif implicit none ! -------------------------------------------------- ! @@ -124,9 +117,6 @@ program rrtmgp_rfmip_sw ! type(ty_gas_concs), dimension(:), allocatable :: gas_conc_array real(wp), parameter :: deg_to_rad = acos(-1._wp)/180._wp -#ifdef USE_TIMING - integer :: ret, i -#endif ! ------------------------------------------------------------------------------------------------- ! ! Code starts @@ -235,20 +225,9 @@ program rrtmgp_rfmip_sw !$acc enter data create (sfc_alb_spec, mu0) !$omp target enter data map(alloc:sfc_alb_spec, mu0) ! -------------------------------------------------- -#ifdef USE_TIMING - ! - ! Initialize timers - ! - ret = gptlsetoption (gptlpercent, 1) ! Turn on "% of" print - ret = gptlsetoption (gptloverhead, 0) ! Turn off overhead estimate - ret = gptlinitialize() -#endif ! ! Loop over blocks ! -#ifdef USE_TIMING - do i = 1, 4 -#endif do b = 1, nblocks fluxes%flux_up => flux_up(:,:,b) fluxes%flux_dn => flux_dn(:,:,b) @@ -256,18 +235,12 @@ program rrtmgp_rfmip_sw ! Compute the optical properties of the atmosphere and the Planck source functions ! from pressures, temperatures, and gas concentrations... ! -#ifdef USE_TIMING - ret = gptlstart('gas_optics (SW)') -#endif call stop_on_err(k_dist%gas_optics(p_lay(:,:,b), & p_lev(:,:,b), & t_lay(:,:,b), & gas_conc_array(b), & optical_props, & toa_flux)) -#ifdef USE_TIMING - ret = gptlstop('gas_optics (SW)') -#endif ! Boundary conditions ! (This is partly to show how to keep work on GPUs using OpenACC in a host application) ! What's the total solar irradiance assumed by RRTMGP? @@ -322,9 +295,6 @@ program rrtmgp_rfmip_sw ! ... and compute the spectrally-resolved fluxes, providing reduced values ! via ty_fluxes_broadband ! -#ifdef USE_TIMING - ret = gptlstart('rte_sw') -#endif call stop_on_err(rte_sw(optical_props, & top_at_1, & mu0, & @@ -332,9 +302,6 @@ program rrtmgp_rfmip_sw sfc_alb_spec, & sfc_alb_spec, & fluxes)) -#ifdef USE_TIMING - ret = gptlstop('rte_sw') -#endif ! ! Zero out fluxes for which the original solar zenith angle is > 90 degrees. ! @@ -348,11 +315,6 @@ program rrtmgp_rfmip_sw ! ! End timers ! -#ifdef USE_TIMING - end do - ret = gptlpr(block_size) - ret = gptlfinalize() -#endif !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g, optical_props) !$omp target exit data map(release:optical_props%tau, optical_props%ssa, optical_props%g) !$acc exit data delete(sfc_alb_spec, mu0) From d2b404f19101dab70ba81101a44c9e4f907c9b26 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 19:10:48 -0400 Subject: [PATCH 02/38] Refactored all-sky example, WIP --- examples/all-sky/Makefile | 4 +- examples/all-sky/rrtmgp_allsky.F90 | 478 +++++++++++++++-------------- 2 files changed, 258 insertions(+), 224 deletions(-) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index 21e129a0d..119ed2b44 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -50,8 +50,8 @@ mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merr tests: cp ${RRTMGP_DATA}/examples/all-sky/inputs/garand-atmos-1.nc rrtmgp-allsky.nc - $(RUN_CMD) ./rrtmgp_allsky rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc 128 - $(RUN_CMD) ./rrtmgp_allsky rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc 128 + $(RUN_CMD) ./rrtmgp_allsky 128 54 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc + $(RUN_CMD) ./rrtmgp_allsky 128 54 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc check: python ./compare-to-reference.py diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 423661271..124be9230 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -89,19 +89,10 @@ program rte_rrtmgp_clouds_aerosols ! Aerosols ! logical :: cell_has_aerosols - integer, dimension(:,:), allocatable :: aero_type - ! MERRA2/GOCART aerosol type - ! 0: no aerosol - ! 1: dust - ! 2: sea salt - ! 3: sulfate - ! 4: black carbon, hydrophobic - ! 5: black carbon, hydrophilic - ! 6: organic carbon, hydrophobic - ! 7: organic carbon, hydrophilic + integer, dimension(:,:), allocatable :: aero_type + ! MERRA2/GOCART aerosol type real(wp), dimension(:,:), allocatable :: aero_size ! Aerosol size for dust and sea salt - ! Allowable range: 0 - 10 microns real(wp), dimension(:,:), allocatable :: aero_mass ! Aerosol mass column (kg/m2) real(wp), dimension(:,:), allocatable :: relhum @@ -133,13 +124,13 @@ program rte_rrtmgp_clouds_aerosols ! logical :: top_at_1, is_sw, is_lw - integer :: ncol, nlay, nbnd, ngpt + integer :: nbnd, ngpt integer :: icol, ilay, ibnd, iloop, igas - real(wp) :: rel_val, rei_val character(len=8) :: char_input - integer :: nUserArgs=0, nloops=1, ncol = 1, nlay=1 - logical :: use_luts = .true., write_fluxes = .true. + integer :: nUserArgs, nloops, ncol, nlay + logical :: do_clouds = .false., use_luts = .true., write_fluxes = .true. + logical :: do_aerosols = .false. integer, parameter :: ngas = 8 character(len=3), dimension(ngas) & :: gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] @@ -157,29 +148,32 @@ program rte_rrtmgp_clouds_aerosols ! Code ! ---------------------------------------------------------------------------------- ! - ! Parse command line for any file names, block size + ! Parse command line: rrtmgp_allsky ncol nlay nreps kdist [clouds [aerosols]] + ! nUserArgs = command_argument_count() - nloops = 1 - if (nUserArgs < 4) call stop_on_err("Need to supply input_file k_distribution_file ncol.") - if (nUserArgs >= 1) call get_command_argument(1,input_file) - if (nUserArgs >= 2) call get_command_argument(2,k_dist_file) - if (nUserArgs >= 3) call get_command_argument(3,cloud_optics_file) - if (nUserArgs >= 4) call get_command_argument(4,aerosol_optics_file) - if (nUserArgs >= 5) then - call get_command_argument(5, char_input) - read(char_input, '(i8)') ncol - if(ncol <= 0) call stop_on_err("Specify positive ncol.") - end if - if (nUserArgs >= 6) then - call get_command_argument(6, char_input) - read(char_input, '(i8)') nloops - if(nloops <= 0) call stop_on_err("Specify positive nloops.") - end if - if (nUserArgs > 7) print *, "Ignoring command line arguments beyond the first six..." - if(trim(input_file) == '-h' .or. trim(input_file) == "--help") then - call stop_on_err("rte_rrtmgp_clouds_aerosols input_file absorption_coefficients_file cloud_optics_file aerosol_optics_file ncol") - end if + print *, "Got so many arugments", nUserArgs + if (nUserArgs < 5) call stop_on_err("Need to supply ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") + call get_command_argument(1, char_input) + read(char_input, '(i8)') ncol + if(ncol <= 0) call stop_on_err("Specify positive ncol.") + call get_command_argument(2, char_input) + read(char_input, '(i8)') nlay + if(nlay <= 0) call stop_on_err("Specify positive nlay.") + call get_command_argument(3, char_input) + read(char_input, '(i8)') nloops + if(nloops <= 0) call stop_on_err("Specify positive nreps (number of times to do ncol examples.") + call get_command_argument(4,input_file) + call get_command_argument(5,k_dist_file) + if (nUserArgs >= 6) then + call get_command_argument(6,cloud_optics_file) + do_clouds = .true. + end if + if (nUserArgs >= 7) then + call get_command_argument(7,aerosol_optics_file) + do_aerosols = .true. + end if + if (nUserArgs > 7) print *, "Ignoring command line arguments beyond the first seven..." ! ! Read temperature, pressure, gas concentrations. @@ -218,27 +212,32 @@ program rte_rrtmgp_clouds_aerosols call load_and_init(k_dist, k_dist_file, gas_concs) is_sw = k_dist%source_is_external() is_lw = .not. is_sw - ! - ! Should also try with Pade calculations - ! call load_cld_padecoeff(cloud_optics, cloud_optics_file) - ! - if(use_luts) then - call load_cld_lutcoeff (cloud_optics, cloud_optics_file) - else - call load_cld_padecoeff(cloud_optics, cloud_optics_file) + if (do_clouds) then + ! + ! Should also try with Pade calculations + ! call load_cld_padecoeff(cloud_optics, cloud_optics_file) + ! + if(use_luts) then + call load_cld_lutcoeff (cloud_optics, cloud_optics_file) + else + call load_cld_padecoeff(cloud_optics, cloud_optics_file) + end if + call stop_on_err(cloud_optics%set_ice_roughness(2)) end if - call stop_on_err(cloud_optics%set_ice_roughness(2)) - ! - ! Load aerosol optics coefficients from lookup tables - ! - call load_aero_lutcoeff (aerosol_optics, aerosol_optics_file) - ! - ! Derive relative humidity from profile - ! - allocate(vmr_h2o(ncol, nlay)) - call stop_on_err(gas_concs%get_vmr(gas_names(1),vmr_h2o)) - allocate(relhum(ncol, nlay)) - call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) + + if (do_aerosols) then + ! + ! Load aerosol optics coefficients from lookup tables + ! + call load_aero_lutcoeff (aerosol_optics, aerosol_optics_file) + ! + ! Derive relative humidity from profile + ! + allocate(vmr_h2o(ncol, nlay)) + call stop_on_err(gas_concs%get_vmr(gas_names(1),vmr_h2o)) + allocate(relhum(ncol, nlay)) + call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) + end if ! ---------------------------------------------------------------------------- ! @@ -254,17 +253,9 @@ program rte_rrtmgp_clouds_aerosols ! if(is_sw) then allocate(ty_optical_props_2str::atmos) - allocate(ty_optical_props_2str::clouds) - allocate(ty_optical_props_2str::aerosols) else allocate(ty_optical_props_1scl::atmos) - allocate(ty_optical_props_1scl::clouds) - allocate(ty_optical_props_1scl::aerosols) end if - ! Clouds optical props are defined by band - call stop_on_err(clouds%init(k_dist%get_band_lims_wavenumber())) - ! Aerosol optical props are defined by band - call stop_on_err(aerosols%init(k_dist%get_band_lims_wavenumber())) ! ! Allocate arrays for the optical properties themselves. ! @@ -281,30 +272,7 @@ program rte_rrtmgp_clouds_aerosols class default call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") end select - select type(clouds) - class is (ty_optical_props_1scl) - call stop_on_err(clouds%alloc_1scl(ncol, nlay)) - !$acc enter data copyin(clouds) create(clouds%tau) - !$omp target enter data map(alloc:clouds%tau) - class is (ty_optical_props_2str) - call stop_on_err(clouds%alloc_2str(ncol, nlay)) - !$acc enter data copyin(clouds) create(clouds%tau, clouds%ssa, clouds%g) - !$omp target enter data map(alloc:clouds%tau, clouds%ssa, clouds%g) - class default - call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") - end select - select type(aerosols) - class is (ty_optical_props_1scl) - call stop_on_err(aerosols%alloc_1scl(ncol, nlay)) - !$acc enter data copyin(aerosols) create(aerosols%tau) - !$omp target enter data map(alloc:aerosols%tau) - class is (ty_optical_props_2str) - call stop_on_err(aerosols%alloc_2str(ncol, nlay)) - !$acc enter data copyin(aerosols) create(aerosols%tau, aerosols%ssa, aerosols%g) - !$omp target enter data map(alloc:aerosols%tau, aerosols%ssa, aerosols%g) - class default - call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") - end select + ! ---------------------------------------------------------------------------- ! Boundary conditions depending on whether the k-distribution being supplied ! is LW or SW @@ -357,84 +325,9 @@ program rte_rrtmgp_clouds_aerosols !$acc enter data create(flux_dir) !$omp target enter data map(alloc:flux_dir) end if - ! - ! Clouds - ! - allocate(lwp(ncol,nlay), iwp(ncol,nlay), & - rel(ncol,nlay), rei(ncol,nlay), cloud_mask(ncol,nlay)) - !$acc enter data create(cloud_mask, lwp, iwp, rel, rei) - !$omp target enter data map(alloc:cloud_mask, lwp, iwp, rel, rei) - - ! Restrict clouds to troposphere (> 100 hPa = 100*100 Pa) - ! and not very close to the ground (< 900 hPa), and - ! put them in 2/3 of the columns since that's roughly the - ! total cloudiness of earth - rel_val = 0.5 * (cloud_optics%get_min_radius_liq() + cloud_optics%get_max_radius_liq()) - rei_val = 0.5 * (cloud_optics%get_min_radius_ice() + cloud_optics%get_max_radius_ice()) - !$acc parallel loop collapse(2) copyin(t_lay) copyout(lwp, iwp, rel, rei) - !$omp target teams distribute parallel do simd collapse(2) map(to:t_lay) map(from:lwp, iwp, rel, rei) - do ilay=1,nlay - do icol=1,ncol - cloud_mask(icol,ilay) = p_lay(icol,ilay) > 100._wp * 100._wp .and. & - p_lay(icol,ilay) < 900._wp * 100._wp .and. & - mod(icol, 3) /= 0 - ! - ! Ice and liquid will overlap in a few layers - ! - lwp(icol,ilay) = merge(10._wp, 0._wp, cloud_mask(icol,ilay) .and. t_lay(icol,ilay) > 263._wp) - iwp(icol,ilay) = merge(10._wp, 0._wp, cloud_mask(icol,ilay) .and. t_lay(icol,ilay) < 273._wp) - rel(icol,ilay) = merge(rel_val, 0._wp, lwp(icol,ilay) > 0._wp) - rei(icol,ilay) = merge(rei_val, 0._wp, iwp(icol,ilay) > 0._wp) - end do - end do - !$acc exit data delete(cloud_mask) - !$omp target exit data map(release:cloud_mask) - ! - ! Aerosols - ! - allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & - aero_mass(ncol,nlay), aero_mask(ncol,nlay)) - !$acc enter data create( aero_mask, aero_type, aero_size, aero_mass) - !$omp target enter data map(alloc:aero_mask, aero_type, aero_size, aero_mass) - - ! Restrict sulfate aerosols to lower stratosphere (> 50 hPa = 50*100 Pa; < 100 hPa = 100*100 Pa) - ! and dust aerosols to the lower troposphere (> 700 hPa; < 900 hPa), and - ! put them in 1/2 of the columns - ! - ! - !$acc parallel loop collapse(2) copyin(p_lay) - !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay) - do ilay=1,nlay - do icol=1,ncol - cell_has_aerosols = ((p_lay(icol,ilay) > 50._wp * 100._wp .and. & - p_lay(icol,ilay) < 100._wp * 100._wp) .or. & - (p_lay(icol,ilay) > 700._wp * 100._wp .and. & - p_lay(icol,ilay) < 900._wp * 100._wp)) .and. & - mod(icol, 2) /= 0 - if (cell_has_aerosols) then - ! Sulfate aerosol - if (p_lay(icol,ilay) > 50._wp * 100._wp .and. & - p_lay(icol,ilay) < 100._wp * 100._wp) then - aero_type(icol,ilay) = merra_aero_sulf - aero_size(icol,ilay) = 0.2_wp - aero_mass(icol,ilay) = 1.e-6_wp - endif - ! Dust aerosol - if (p_lay(icol,ilay) > 700._wp * 100._wp .and. & - p_lay(icol,ilay) < 900._wp * 100._wp) then - aero_type(icol,ilay) = merra_aero_dust - aero_size(icol,ilay) = 0.5_wp - aero_mass(icol,ilay) = 3.e-5_wp - endif - else - aero_type(icol,ilay) = 0 - aero_size(icol,ilay) = 0._wp - aero_mass(icol,ilay) = 0._wp - end if - aero_mask(icol, ilay) = cell_has_aerosols - end do - end do + if (do_clouds) call compute_clouds + if (do_aerosols) call compute_aerosols ! ---------------------------------------------------------------------------- ! @@ -454,14 +347,14 @@ program rte_rrtmgp_clouds_aerosols ! ! Cloud optics ! - call stop_on_err( & - cloud_optics%cloud_optics(lwp, iwp, rel, rei, clouds)) + if(do_clouds) & + call stop_on_err(cloud_optics%cloud_optics(lwp, iwp, rel, rei, clouds)) ! ! Aerosol optics ! - call stop_on_err( & - aerosol_optics%aerosol_optics(aero_type, aero_size, & - aero_mass, relhum, aerosols)) + if(do_aerosols) & + call stop_on_err(aerosol_optics%aerosol_optics(aero_type, aero_size, & + aero_mass, relhum, aerosols)) ! ! Solvers ! @@ -478,8 +371,8 @@ program rte_rrtmgp_clouds_aerosols atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(clouds%increment(atmos)) - call stop_on_err(aerosols%increment(atmos)) + 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, & @@ -497,10 +390,14 @@ program rte_rrtmgp_clouds_aerosols gas_concs, & atmos, & toa_flux)) - call stop_on_err(clouds%delta_scale()) - call stop_on_err(clouds%increment(atmos)) - call stop_on_err(aerosols%delta_scale()) - call stop_on_err(aerosols%increment(atmos)) + if(do_clouds) then + call stop_on_err(clouds%delta_scale()) + call stop_on_err(clouds%increment(atmos)) + end if + if(do_aerosols) then + 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, & sfc_alb_dir, sfc_alb_dif, & @@ -515,8 +412,13 @@ program rte_rrtmgp_clouds_aerosols ! call system_clock(finish_all, clock_rate) ! - !$acc exit data delete(lwp, iwp, rel, rei) - !$omp target exit data map(release:lwp, iwp, rel, rei) + if(do_clouds) then + !$acc exit data delete( lwp, iwp, rel, rei) + !$omp target exit data map(release:lwp, iwp, rel, rei) + end if + ! + ! TK - release aerosol memory + ! !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) @@ -544,50 +446,182 @@ program rte_rrtmgp_clouds_aerosols !$acc exit data delete(sfc_alb_dir, sfc_alb_dif, mu0) !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if - !$acc enter data create(lwp, iwp, rel, rei) - !$omp target enter data map(alloc:lwp, iwp, rel, rei) contains -! -------------------------------------------------------------------------------------- -! -! Calculate layer relative humidity for aerosol optics calculations -! -subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) - use mo_rte_kind, only: wp - use mo_gas_optics_constants, only: m_h2o, m_dry - - integer, intent(in) :: ncol, nlay - real(wp), intent(in) :: p_lay(ncol,nlay) ! layer pressure (Pa) - real(wp), intent(in) :: t_lay(ncol,nlay) ! layer temperature (K) - real(wp), intent(in) :: vmr_h2o(ncol,nlay) ! water volume mixing ratio - - real(wp), intent(inout) :: relhum(ncol,nlay) ! relative humidity (fraction, 0-1) - - ! Local variables - integer :: i, k - - real(wp) :: mmr_h2o ! water mass mixing ratio - real(wp) :: q_lay ! water specific humidity - real(wp) :: q_lay_min, q_tmp, es_tmp - real(wp) :: mwd, t_ref, rh - - ! Set constants - mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights - t_ref = 273.16_wp ! reference temperature (K) - q_lay_min = 1.e-7_wp ! minimum water mass mixing ratio - ! ------------------- - - ! Derive layer virtual temperature - do i = 1, ncol - do k = 1, nlay - ! Convert h2o vmr to mmr - mmr_h2o = vmr_h2o(i,k) * mwd - q_lay = mmr_h2o / (1 + mmr_h2o) - q_tmp = max(q_lay_min, q_lay) - es_tmp = exp( (17.67_wp * (t_lay(i,k)-t_ref)) / (t_lay(i,k)-29.65_wp) ) - rh = (0.263_wp * p_lay(i,k) * q_tmp) / es_tmp - ! Convert rh from percent to fraction - relhum(i,k) = 0.01_wp * rh - enddo - enddo + ! -------------------------------------------------------------------------------------- + ! + subroutine compute_clouds + real(wp) :: rel_val, rei_val + ! + ! Variable and memory allocation + ! + if(is_sw) then + allocate(ty_optical_props_2str::clouds) + else + allocate(ty_optical_props_1scl::clouds) + end if + ! Clouds optical props are defined by band + call stop_on_err(clouds%init(k_dist%get_band_lims_wavenumber())) + + select type(clouds) + class is (ty_optical_props_1scl) + call stop_on_err(clouds%alloc_1scl(ncol, nlay)) + !$acc enter data copyin(clouds) create(clouds%tau) + !$omp target enter data map(alloc:clouds%tau) + class is (ty_optical_props_2str) + call stop_on_err(clouds%alloc_2str(ncol, nlay)) + !$acc enter data copyin(clouds) create(clouds%tau, clouds%ssa, clouds%g) + !$omp target enter data map(alloc:clouds%tau, clouds%ssa, clouds%g) + class default + call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") + end select + ! + ! Cloud physical properties + ! + allocate(lwp(ncol,nlay), iwp(ncol,nlay), & + rel(ncol,nlay), rei(ncol,nlay), cloud_mask(ncol,nlay)) + !$acc enter data create( cloud_mask, lwp, iwp, rel, rei) + !$omp target enter data map(alloc:cloud_mask, lwp, iwp, rel, rei) + + ! Restrict clouds to troposphere (> 100 hPa = 100*100 Pa) + ! and not very close to the ground (< 900 hPa), and + ! put them in 2/3 of the columns since that's roughly the + ! total cloudiness of earth + rel_val = 0.5 * (cloud_optics%get_min_radius_liq() + cloud_optics%get_max_radius_liq()) + rei_val = 0.5 * (cloud_optics%get_min_radius_ice() + cloud_optics%get_max_radius_ice()) + !$acc parallel loop collapse(2) copyin(t_lay) copyout( lwp, iwp, rel, rei) + !$omp target teams distribute parallel do simd collapse(2) map(to:t_lay) map(from:lwp, iwp, rel, rei) + do ilay=1,nlay + do icol=1,ncol + cloud_mask(icol,ilay) = p_lay(icol,ilay) > 100._wp * 100._wp .and. & + p_lay(icol,ilay) < 900._wp * 100._wp .and. & + mod(icol, 3) /= 0 + ! + ! Ice and liquid will overlap in a few layers + ! + lwp(icol,ilay) = merge(10._wp, 0._wp, cloud_mask(icol,ilay) .and. t_lay(icol,ilay) > 263._wp) + iwp(icol,ilay) = merge(10._wp, 0._wp, cloud_mask(icol,ilay) .and. t_lay(icol,ilay) < 273._wp) + rel(icol,ilay) = merge(rel_val, 0._wp, lwp(icol,ilay) > 0._wp) + rei(icol,ilay) = merge(rei_val, 0._wp, iwp(icol,ilay) > 0._wp) + end do + end do + !$acc exit data delete(cloud_mask) + !$omp target exit data map(release:cloud_mask) + + end subroutine compute_clouds + ! + ! -------------------------------------------------------------------------------------- + ! + subroutine compute_aerosols + ! + ! Variable and memory allocation + ! + if(is_sw) then + allocate(ty_optical_props_2str::aerosols) + else + allocate(ty_optical_props_1scl::aerosols) + end if + call stop_on_err(aerosols%init(k_dist%get_band_lims_wavenumber())) + select type(aerosols) + class is (ty_optical_props_1scl) + call stop_on_err(aerosols%alloc_1scl(ncol, nlay)) + !$acc enter data copyin(aerosols) create(aerosols%tau) + !$omp target enter data map(alloc:aerosols%tau) + class is (ty_optical_props_2str) + call stop_on_err(aerosols%alloc_2str(ncol, nlay)) + !$acc enter data copyin(aerosols) create(aerosols%tau, aerosols%ssa, aerosols%g) + !$omp target enter data map(alloc:aerosols%tau, aerosols%ssa, aerosols%g) + class default + call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") + end select + + ! + ! Aerosols + ! + allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & + aero_mass(ncol,nlay), aero_mask(ncol,nlay)) + !$acc enter data create( aero_mask, aero_type, aero_size, aero_mass) + !$omp target enter data map(alloc:aero_mask, aero_type, aero_size, aero_mass) + + ! Restrict sulfate aerosols to lower stratosphere (> 50 hPa = 50*100 Pa; < 100 hPa = 100*100 Pa) + ! and dust aerosols to the lower troposphere (> 700 hPa; < 900 hPa), and + ! put them in 1/2 of the columns + ! + ! + !$acc parallel loop collapse(2) copyin(p_lay) + !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay) + do ilay=1,nlay + do icol=1,ncol + cell_has_aerosols = ((p_lay(icol,ilay) > 50._wp * 100._wp .and. & + p_lay(icol,ilay) < 100._wp * 100._wp) .or. & + (p_lay(icol,ilay) > 700._wp * 100._wp .and. & + p_lay(icol,ilay) < 900._wp * 100._wp)) .and. & + mod(icol, 2) /= 0 + if (cell_has_aerosols) then + ! Sulfate aerosol + if (p_lay(icol,ilay) > 50._wp * 100._wp .and. & + p_lay(icol,ilay) < 100._wp * 100._wp) then + aero_type(icol,ilay) = merra_aero_sulf + aero_size(icol,ilay) = 0.2_wp + aero_mass(icol,ilay) = 1.e-6_wp + endif + ! Dust aerosol + if (p_lay(icol,ilay) > 700._wp * 100._wp .and. & + p_lay(icol,ilay) < 900._wp * 100._wp) then + aero_type(icol,ilay) = merra_aero_dust + aero_size(icol,ilay) = 0.5_wp + aero_mass(icol,ilay) = 3.e-5_wp + endif + else + aero_type(icol,ilay) = 0 + aero_size(icol,ilay) = 0._wp + aero_mass(icol,ilay) = 0._wp + end if + aero_mask(icol, ilay) = cell_has_aerosols + end do + end do + + end subroutine compute_aerosols + ! -------------------------------------------------------------------------------------- + ! + ! Calculate layer relative humidity for aerosol optics calculations + ! + subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) + use mo_rte_kind, only: wp + use mo_gas_optics_constants, only: m_h2o, m_dry + + integer, intent(in) :: ncol, nlay + real(wp), intent(in) :: p_lay(ncol,nlay) ! layer pressure (Pa) + real(wp), intent(in) :: t_lay(ncol,nlay) ! layer temperature (K) + real(wp), intent(in) :: vmr_h2o(ncol,nlay) ! water volume mixing ratio + + real(wp), intent(inout) :: relhum(ncol,nlay) ! relative humidity (fraction, 0-1) + + ! Local variables + integer :: i, k + + real(wp) :: mmr_h2o ! water mass mixing ratio + real(wp) :: q_lay ! water specific humidity + real(wp) :: q_lay_min, q_tmp, es_tmp + real(wp) :: mwd, t_ref, rh + + ! Set constants + mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights + t_ref = 273.16_wp ! reference temperature (K) + q_lay_min = 1.e-7_wp ! minimum water mass mixing ratio + ! ------------------- + + ! Derive layer virtual temperature + do i = 1, ncol + do k = 1, nlay + ! Convert h2o vmr to mmr + mmr_h2o = vmr_h2o(i,k) * mwd + q_lay = mmr_h2o / (1 + mmr_h2o) + q_tmp = max(q_lay_min, q_lay) + es_tmp = exp( (17.67_wp * (t_lay(i,k)-t_ref)) / (t_lay(i,k)-29.65_wp) ) + rh = (0.263_wp * p_lay(i,k) * q_tmp) / es_tmp + ! Convert rh from percent to fraction + relhum(i,k) = 0.01_wp * rh + enddo + enddo end subroutine get_relhum end program rte_rrtmgp_clouds_aerosols From 1f6b934525adaa15e7d78f6880f2e52067fb87d8 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 19:12:53 -0400 Subject: [PATCH 03/38] No GPTL in Makefile either --- examples/rfmip-clear-sky/Makefile | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index ee6fe7de7..a0c05c8dd 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -17,22 +17,6 @@ FCINCLUDE += -I$(NFHOME)/include LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib LIBS += -lnetcdff -lnetcdf -# -# General Purpose Timing Library https://jmrosinski.github.io/GPTL/ -# Set environment variable GPTL_DIR to the root of a GPTL installation to build -# the RFMIP example with timers -# -ifneq ($(origin GPTL_DIR),undefined) - # - # Timing library - # - FCINCLUDE += -I$(GPTL_DIR)/include - # Compiler specific - FCFLAGS += -DUSE_TIMING - LDFLAGS += -L$(GPTL_DIR)/lib - LIBS += -lgptl -endif - VPATH = ../ # Compilation rules From d159a391dc4843acc93855a4acce11d0ff6123ca Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 21:01:39 -0400 Subject: [PATCH 04/38] Refactor RH - why's nvfortran failing? --- examples/all-sky/rrtmgp_allsky.F90 | 58 +++++++++++++----------------- 1 file changed, 24 insertions(+), 34 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 124be9230..663decdce 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -99,8 +99,6 @@ program rte_rrtmgp_clouds_aerosols ! Relative humidity (fraction) logical, dimension(:,:), allocatable :: aero_mask ! Aerosol mask - real(wp), allocatable, dimension(:,:) :: vmr_h2o - ! h2o vmr ! ! Output variables @@ -152,7 +150,6 @@ program rte_rrtmgp_clouds_aerosols ! nUserArgs = command_argument_count() - print *, "Got so many arugments", nUserArgs if (nUserArgs < 5) call stop_on_err("Need to supply ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") call get_command_argument(1, char_input) read(char_input, '(i8)') ncol @@ -230,13 +227,6 @@ program rte_rrtmgp_clouds_aerosols ! Load aerosol optics coefficients from lookup tables ! call load_aero_lutcoeff (aerosol_optics, aerosol_optics_file) - ! - ! Derive relative humidity from profile - ! - allocate(vmr_h2o(ncol, nlay)) - call stop_on_err(gas_concs%get_vmr(gas_names(1),vmr_h2o)) - allocate(relhum(ncol, nlay)) - call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) end if ! ---------------------------------------------------------------------------- @@ -512,7 +502,9 @@ end subroutine compute_clouds ! -------------------------------------------------------------------------------------- ! subroutine compute_aerosols - ! + real(wp), dimension(ncol,nlay) :: vmr_h2o ! h2o vmr + logical :: is_sulfate, is_dust, is_even_column + ! ! Variable and memory allocation ! if(is_sw) then @@ -533,10 +525,15 @@ subroutine compute_aerosols class default call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") end select - ! - ! Aerosols + ! Derive relative humidity from profile + ! + call stop_on_err(gas_concs%get_vmr("h2o",vmr_h2o)) + allocate(relhum(ncol, nlay)) + call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! + ! Aerosol properties + ! allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & aero_mass(ncol,nlay), aero_mask(ncol,nlay)) !$acc enter data create( aero_mask, aero_type, aero_size, aero_mass) @@ -551,32 +548,25 @@ subroutine compute_aerosols !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay) do ilay=1,nlay do icol=1,ncol - cell_has_aerosols = ((p_lay(icol,ilay) > 50._wp * 100._wp .and. & - p_lay(icol,ilay) < 100._wp * 100._wp) .or. & - (p_lay(icol,ilay) > 700._wp * 100._wp .and. & - p_lay(icol,ilay) < 900._wp * 100._wp)) .and. & - mod(icol, 2) /= 0 - if (cell_has_aerosols) then - ! Sulfate aerosol - if (p_lay(icol,ilay) > 50._wp * 100._wp .and. & - p_lay(icol,ilay) < 100._wp * 100._wp) then - aero_type(icol,ilay) = merra_aero_sulf - aero_size(icol,ilay) = 0.2_wp - aero_mass(icol,ilay) = 1.e-6_wp - endif + is_sulfate = (p_lay(icol,ilay) > 50._wp * 100._wp .and. & + p_lay(icol,ilay) < 100._wp * 100._wp) + is_dust = (p_lay(icol,ilay) > 700._wp * 100._wp .and. & + p_lay(icol,ilay) < 900._wp * 100._wp) + is_even_column = mod(icol, 2) /= 0 + if (is_even_column .and. is_sulfate) then + aero_type(icol,ilay) = merra_aero_sulf + aero_size(icol,ilay) = 0.2_wp + aero_mass(icol,ilay) = 1.e-6_wp + else if(is_even_column .and. is_dust) then ! Dust aerosol - if (p_lay(icol,ilay) > 700._wp * 100._wp .and. & - p_lay(icol,ilay) < 900._wp * 100._wp) then - aero_type(icol,ilay) = merra_aero_dust - aero_size(icol,ilay) = 0.5_wp - aero_mass(icol,ilay) = 3.e-5_wp - endif - else + aero_type(icol,ilay) = merra_aero_dust + aero_size(icol,ilay) = 0.5_wp + aero_mass(icol,ilay) = 3.e-5_wp + else aero_type(icol,ilay) = 0 aero_size(icol,ilay) = 0._wp aero_mass(icol,ilay) = 0._wp end if - aero_mask(icol, ilay) = cell_has_aerosols end do end do From 1825269aec90d85ea8d21c72b35a6f1d91515157 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 22:34:33 -0400 Subject: [PATCH 05/38] Debugging prints --- examples/all-sky/rrtmgp_allsky.F90 | 3 ++- gas-optics/mo_gas_concentrations.F90 | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 663decdce..67684035c 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -22,6 +22,7 @@ subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) !$acc data create( tmp, tmp_col) !$omp target data map(alloc:tmp, tmp_col) + print *, "vmr_2d_to_1d " // name call stop_on_err(gas_concs_garand%get_vmr(name, tmp)) !$acc kernels !$omp target @@ -562,7 +563,7 @@ subroutine compute_aerosols aero_type(icol,ilay) = merra_aero_dust aero_size(icol,ilay) = 0.5_wp aero_mass(icol,ilay) = 3.e-5_wp - else + else aero_type(icol,ilay) = 0 aero_size(icol,ilay) = 0._wp aero_mass(icol,ilay) = 0._wp diff --git a/gas-optics/mo_gas_concentrations.F90 b/gas-optics/mo_gas_concentrations.F90 index 26de09998..faa7c5533 100644 --- a/gas-optics/mo_gas_concentrations.F90 +++ b/gas-optics/mo_gas_concentrations.F90 @@ -336,6 +336,7 @@ function get_vmr_1d(this, gas, array) result(error_msg) if(this%nlay > 0 .and. this%nlay /= size(array)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' + print *, this%nlay, size(array) end if if(error_msg /= "") return @@ -396,6 +397,7 @@ function get_vmr_2d(this, gas, array) result(error_msg) end if if(this%nlay > 0 .and. this%nlay /= size(array,2)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' + print *, this%nlay, size(array, 2) end if if(error_msg /= "") return From feeff71011afd0689af7f8f692a9bd9b7b47b4f2 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 22:53:15 -0400 Subject: [PATCH 06/38] More debuggin --- examples/all-sky/rrtmgp_allsky.F90 | 1 + gas-optics/mo_gas_concentrations.F90 | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 67684035c..3f357f1e3 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -23,6 +23,7 @@ subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) !$acc data create( tmp, tmp_col) !$omp target data map(alloc:tmp, tmp_col) print *, "vmr_2d_to_1d " // name + print *, " sizes", sz1, sz2 call stop_on_err(gas_concs_garand%get_vmr(name, tmp)) !$acc kernels !$omp target diff --git a/gas-optics/mo_gas_concentrations.F90 b/gas-optics/mo_gas_concentrations.F90 index faa7c5533..0ae9965b2 100644 --- a/gas-optics/mo_gas_concentrations.F90 +++ b/gas-optics/mo_gas_concentrations.F90 @@ -336,7 +336,7 @@ function get_vmr_1d(this, gas, array) result(error_msg) if(this%nlay > 0 .and. this%nlay /= size(array)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' - print *, this%nlay, size(array) + print *, "get_vmr_1d", this%nlay, size(array) end if if(error_msg /= "") return @@ -397,7 +397,7 @@ function get_vmr_2d(this, gas, array) result(error_msg) end if if(this%nlay > 0 .and. this%nlay /= size(array,2)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' - print *, this%nlay, size(array, 2) + print *, "get_vmr_2d", this%nlay, size(array, 2) end if if(error_msg /= "") return From f4fe1711d90f692403e4278ec4ce4d7a9fd08c56 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 23:04:12 -0400 Subject: [PATCH 07/38] Even more debuggin' --- examples/all-sky/rrtmgp_allsky.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 3f357f1e3..361c00eb6 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -183,6 +183,7 @@ program rte_rrtmgp_clouds_aerosols gas_concs_garand, col_dry) deallocate(col_dry) nlay = size(p_lay, 2) + print *, "In main, Garand atmos sizes : ", size(p_lay, 1), size(p_lay, 2) ! For clouds we'll use the first column, repeated over and over call stop_on_err(gas_concs%init(gas_names)) From 3e3685bd65ecff1d2da609acd95c07230786bcb4 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 23:23:33 -0400 Subject: [PATCH 08/38] Explicit allocation? --- examples/all-sky/mo_garand_atmos_io.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/all-sky/mo_garand_atmos_io.F90 b/examples/all-sky/mo_garand_atmos_io.F90 index 5b730c9ab..3ddcf7ba0 100644 --- a/examples/all-sky/mo_garand_atmos_io.F90 +++ b/examples/all-sky/mo_garand_atmos_io.F90 @@ -70,6 +70,7 @@ subroutine read_atmos(fileName, & ! allocating on assignment. This may require explicit compiler support ! e.g. -assume realloc_lhs flag for Intel ! + allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlev), t_lev(ncol, nlev)) p_lay = read_field(ncid, 'p_lay', ncol, nlay) t_lay = read_field(ncid, 't_lay', ncol, nlay) p_lev = read_field(ncid, 'p_lev', ncol, nlev) From 35d310b18c709e11cdb35f01e40e918feca846db Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 13 Jun 2023 23:37:28 -0400 Subject: [PATCH 09/38] EMD --- examples/all-sky/mo_garand_atmos_io.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/all-sky/mo_garand_atmos_io.F90 b/examples/all-sky/mo_garand_atmos_io.F90 index 3ddcf7ba0..dde3bfbf0 100644 --- a/examples/all-sky/mo_garand_atmos_io.F90 +++ b/examples/all-sky/mo_garand_atmos_io.F90 @@ -76,6 +76,9 @@ subroutine read_atmos(fileName, & p_lev = read_field(ncid, 'p_lev', ncol, nlev) t_lev = read_field(ncid, 't_lev', ncol, nlev) + print *, "in read_atmos: ncol, nlay are ", ncol, nlay + print *, "in read_atmos: p_lay is ", size(p_lay, 1), size(p_lay, 2) + call stop_on_err(gas_concs%init(gas_names)) do igas = 1, ngas vmr_name = 'vmr_' // trim(gas_names(igas)) From 89b6245acdff1407569f24e3b52b939ba9b7e079 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 14 Jun 2023 13:47:07 -0400 Subject: [PATCH 10/38] Put all subroutines with `contains`? --- examples/all-sky/rrtmgp_allsky.F90 | 81 ++++++++++++++++-------------- 1 file changed, 43 insertions(+), 38 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 361c00eb6..99a699992 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -1,41 +1,3 @@ -subroutine stop_on_err(error_msg) - use iso_fortran_env, only : error_unit - character(len=*), intent(in) :: error_msg - - if(error_msg /= "") then - write (error_unit,*) trim(error_msg) - write (error_unit,*) "rte_rrtmgp_clouds_aerosols stopping" - error stop 1 - end if -end subroutine stop_on_err - -subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) - use mo_gas_concentrations, only: ty_gas_concs - use mo_rte_kind, only: wp - - type(ty_gas_concs), intent(in) :: gas_concs_garand - type(ty_gas_concs), intent(inout) :: gas_concs - character(len=*), intent(in) :: name - integer, intent(in) :: sz1, sz2 - - real(wp) :: tmp(sz1, sz2), tmp_col(sz2) - - !$acc data create( tmp, tmp_col) - !$omp target data map(alloc:tmp, tmp_col) - print *, "vmr_2d_to_1d " // name - print *, " sizes", sz1, sz2 - call stop_on_err(gas_concs_garand%get_vmr(name, tmp)) - !$acc kernels - !$omp target - tmp_col(:) = tmp(1, :) - !$acc end kernels - !$omp end target - - call stop_on_err(gas_concs%set_vmr (name, tmp_col)) - !$acc end data - !$omp end target data -end subroutine vmr_2d_to_1d -! ---------------------------------------------------------------------------------- program rte_rrtmgp_clouds_aerosols use mo_rte_kind, only: wp, i8, wl use mo_optical_props, only: ty_optical_props, & @@ -153,17 +115,22 @@ program rte_rrtmgp_clouds_aerosols ! nUserArgs = command_argument_count() if (nUserArgs < 5) call stop_on_err("Need to supply ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") + call get_command_argument(1, char_input) read(char_input, '(i8)') ncol if(ncol <= 0) call stop_on_err("Specify positive ncol.") + call get_command_argument(2, char_input) read(char_input, '(i8)') nlay if(nlay <= 0) call stop_on_err("Specify positive nlay.") + call get_command_argument(3, char_input) read(char_input, '(i8)') nloops if(nloops <= 0) call stop_on_err("Specify positive nreps (number of times to do ncol examples.") + call get_command_argument(4,input_file) call get_command_argument(5,k_dist_file) + if (nUserArgs >= 6) then call get_command_argument(6,cloud_optics_file) do_clouds = .true. @@ -440,6 +407,44 @@ program rte_rrtmgp_clouds_aerosols !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if contains + ! ---------------------------------------------------------------------------------- + subroutine stop_on_err(error_msg) + use iso_fortran_env, only : error_unit + character(len=*), intent(in) :: error_msg + + if(error_msg /= "") then + write (error_unit,*) trim(error_msg) + write (error_unit,*) "rrtmgp_allsky stopping" + error stop 1 + end if + end subroutine stop_on_err + ! ---------------------------------------------------------------------------------- + subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) + use mo_gas_concentrations, only: ty_gas_concs + use mo_rte_kind, only: wp + + type(ty_gas_concs), intent(in) :: gas_concs_garand + type(ty_gas_concs), intent(inout) :: gas_concs + character(len=*), intent(in) :: name + integer, intent(in) :: sz1, sz2 + + real(wp) :: tmp(sz1, sz2), tmp_col(sz2) + + !$acc data create( tmp, tmp_col) + !$omp target data map(alloc:tmp, tmp_col) + print *, "vmr_2d_to_1d " // name + print *, " sizes", sz1, sz2 + call stop_on_err(gas_concs_garand%get_vmr(name, tmp)) + !$acc kernels + !$omp target + tmp_col(:) = tmp(1, :) + !$acc end kernels + !$omp end target + + call stop_on_err(gas_concs%set_vmr (name, tmp_col)) + !$acc end data + !$omp end target data + end subroutine vmr_2d_to_1d ! -------------------------------------------------------------------------------------- ! subroutine compute_clouds From 80272b93cad35f61385ea3d725d1adf163a19401 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 14 Jun 2023 16:57:21 -0400 Subject: [PATCH 11/38] Add RCE profiles to allsky example (unused so far) --- examples/all-sky/rrtmgp_allsky.F90 | 70 ++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 99a699992..b38e20faa 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -407,6 +407,76 @@ program rte_rrtmgp_clouds_aerosols !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if contains + ! ---------------------------------------------------------------------------------- + subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) + real(wp), intent(in ) :: SST + integer, intent(in ) :: ncol, nlay + real(wp), dimension(ncol, nlay ), intent(out) :: p_lay, t_lay, q, o3 + real(wp), dimension(ncol, nlay+1), intent(out) :: p_lev, t_lev + + real(wp) :: z_lay(nlay), z_lev(nlay+1) + real(wp) :: q_tmp, p_hpa + integer :: icol, ilay, i + + real(wp), parameter :: z_trop = 15000._wp, z_top = 70.e3_wp + ! Ozone profile - maybe only a single profile? + real(wp), parameter :: g1 = 3.6478_wp, g2 = 0.83209_wp, g3 = 11.3515_wp, o3_min = 1e-13_wp + ! According to CvH RRTMGP in Single Precision will fail with lower ozone concentrations + + z_lev(:) = [0._wp, 2._wp* z_trop /nlay * [(i, i=1, nlay/2)], & + z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] + z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) + + !$acc data copyin(z_lev, lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target data map(to:z_lev, lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + + !$acc parallel loop collapse(2) + !$omp target teams distribute parallel do simd collapse(2) + do ilay = 1, nlay + do icol = 1, ncol + call thermo(z_lay(ilay), SST, p_lay(icol,ilay), t_lay(icol,ilay), q(icol,ilay)) + p_hpa = p_lay(icol,ilay) / 100._wp + o3(icol, ilay) = max(o3_min, & + g1 * p_hpa**g2 * exp(-p_hpa/g3) * 1.e-6_wp) + call thermo(z_lev(ilay), SST, p_lev(icol,ilay), t_lev(icol,ilay), q_tmp) + end do + end do + + !$acc parallel loop + !$omp target teams distribute parallel do simd + do icol = 1, ncol + call thermo(z_lev(nlay+1), SST, p_lev(icol, nlay+1), t_lev(icol, nlay+1), q_tmp) + end do + !$acc end data + !$omp target end data + end subroutine compute_profiles + ! --------------------------------------- + elemental function Tv(T, q) + real(wp), intent(in) :: T, q + real(wp) :: Tv + + Tv = (1. + 0.608*q) * T + end function Tv + ! --------------------------------------- + elemental subroutine thermo(z, SST, p, T, q) + real(wp), intent(in) :: z, SST ! Height, SST + real(wp), intent(out) :: p, T, q + + real(wp), parameter :: g = 9.79764, Rd = 287.04, p0 = 101480. ! Surface pressure + real(wp), parameter :: z_q1 = 4.0e3, z_q2 = 7.5e3, z_trop = 15.e3, q_t = 1.e-8 + real(wp), parameter :: gamma = 6.7e-3 + + real(wp), parameter :: q_0 = 0.01864 ! for 300 K SST. + + q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) + T = SST - gamma*z / (1. + 0.608*q) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) + if (z > z_trop) then + q = q_t + T = SST - gamma*z_trop/(1. + 0.608*q_0) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + end if + end subroutine thermo ! ---------------------------------------------------------------------------------- subroutine stop_on_err(error_msg) use iso_fortran_env, only : error_unit From 84a3e07064e4176003795ef5f17e4958137615cc Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 14 Jun 2023 21:31:09 -0400 Subject: [PATCH 12/38] All-sky calculations use RCEMIP-ish profiles - no data written so nothing checked --- examples/all-sky/Makefile | 7 +- examples/all-sky/rrtmgp_allsky.F90 | 125 ++++++++++++++++++--------- gas-optics/mo_gas_concentrations.F90 | 3 - 3 files changed, 86 insertions(+), 49 deletions(-) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index 119ed2b44..79220e2af 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -50,11 +50,12 @@ mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merr tests: cp ${RRTMGP_DATA}/examples/all-sky/inputs/garand-atmos-1.nc rrtmgp-allsky.nc - $(RUN_CMD) ./rrtmgp_allsky 128 54 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc - $(RUN_CMD) ./rrtmgp_allsky 128 54 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc check: - python ./compare-to-reference.py + echo "No checking of all sky results at this time" + # python ./compare-to-reference.py clean: -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index b38e20faa..69738ddce 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -22,14 +22,13 @@ program rte_rrtmgp_clouds_aerosols ! Variables ! ---------------------------------------------------------------------------------- ! Arrays: dimensions (col, lay) - real(wp), dimension(:,:), allocatable :: p_lay, t_lay, p_lev - real(wp), dimension(:,:), allocatable :: col_dry + real(wp), dimension(:,:), allocatable :: p_lay, t_lay, p_lev, t_lev ! t_lev is only needed for LW + real(wp), dimension(:,:), allocatable :: q, o3, col_dry real(wp), dimension(:,:), allocatable :: temp_array ! ! Longwave only ! - real(wp), dimension(:,:), allocatable :: t_lev real(wp), dimension(:), allocatable :: t_sfc real(wp), dimension(:,:), allocatable :: emis_sfc ! First dimension is band ! @@ -38,6 +37,11 @@ program rte_rrtmgp_clouds_aerosols real(wp), dimension(:), allocatable :: mu0 real(wp), dimension(:,:), allocatable :: sfc_alb_dir, sfc_alb_dif ! First dimension is band ! + ! Gas concentrations + ! + character(len=3), dimension(8), parameter :: & + gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] + ! ! Source functions ! ! Longwave @@ -91,11 +95,13 @@ program rte_rrtmgp_clouds_aerosols character(len=8) :: char_input integer :: nUserArgs, nloops, ncol, nlay - logical :: do_clouds = .false., use_luts = .true., write_fluxes = .true. + logical :: write_fluxes = .false. + logical :: do_clouds = .false., use_luts = .true. logical :: do_aerosols = .false. + logical :: do_rce = .true. integer, parameter :: ngas = 8 - character(len=3), dimension(ngas) & - :: gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] +! character(len=3), dimension(ngas) & +! :: gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] character(len=256) :: input_file, k_dist_file, cloud_optics_file, aerosol_optics_file ! @@ -140,40 +146,57 @@ program rte_rrtmgp_clouds_aerosols do_aerosols = .true. end if if (nUserArgs > 7) print *, "Ignoring command line arguments beyond the first seven..." - - ! - ! Read temperature, pressure, gas concentrations. - ! Arrays are allocated as they are read - ! - call read_atmos(input_file, & - p_lay, t_lay, p_lev, t_lev, & - gas_concs_garand, col_dry) - deallocate(col_dry) - nlay = size(p_lay, 2) - print *, "In main, Garand atmos sizes : ", size(p_lay, 1), size(p_lay, 2) - - ! For clouds we'll use the first column, repeated over and over - call stop_on_err(gas_concs%init(gas_names)) - do igas = 1, ngas - call vmr_2d_to_1d(gas_concs, gas_concs_garand, gas_names(igas), size(p_lay, 1), nlay) - end do - ! If we trusted in Fortran allocate-on-assign we could skip the temp_array here - allocate(temp_array(ncol, nlay)) - temp_array = spread(p_lay(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, p_lay) - allocate(temp_array(ncol, nlay)) - temp_array = spread(t_lay(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, t_lay) - allocate(temp_array(ncol, nlay+1)) - temp_array = spread(p_lev(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, p_lev) - allocate(temp_array(ncol, nlay+1)) - temp_array = spread(t_lev(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, t_lev) - - ! This puts pressure and temperature arrays on the GPU - !$acc enter data copyin(p_lay, p_lev, t_lay, t_lev) - !$omp target enter data map(to:p_lay, p_lev, t_lay, t_lev) + ! ----------------------------------------------------------------------------------- + if(do_rce) then + allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) + allocate(q (ncol, nlay), o3(ncol, nlay)) + !$acc data create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) + + call stop_on_err(gas_concs%init(gas_names)) + call stop_on_err(gas_concs%set_vmr("h2o", q )) + call stop_on_err(gas_concs%set_vmr("o3", o3)) + call stop_on_err(gas_concs%set_vmr("co2", 348.e-6_wp)) + call stop_on_err(gas_concs%set_vmr("ch4", 1650.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2o", 306.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2", 0.7808_wp)) + call stop_on_err(gas_concs%set_vmr("o2", 0.2095_wp)) + call stop_on_err(gas_concs%set_vmr("co", 0._wp)) + else + ! + ! Read temperature, pressure, gas concentrations. + ! Arrays are allocated as they are read + ! + call read_atmos(input_file, & + p_lay, t_lay, p_lev, t_lev, & + gas_concs_garand, col_dry) + deallocate(col_dry) + nlay = size(p_lay, 2) + + ! For clouds we'll use the first column, repeated over and over + call stop_on_err(gas_concs%init(gas_names)) + do igas = 1, ngas + call vmr_2d_to_1d(gas_concs, gas_concs_garand, gas_names(igas), size(p_lay, 1), nlay) + end do + ! If we trusted in Fortran allocate-on-assign we could skip the temp_array here + allocate(temp_array(ncol, nlay)) + temp_array = spread(p_lay(1,:), dim = 1, ncopies=ncol) + call move_alloc(temp_array, p_lay) + allocate(temp_array(ncol, nlay)) + temp_array = spread(t_lay(1,:), dim = 1, ncopies=ncol) + call move_alloc(temp_array, t_lay) + allocate(temp_array(ncol, nlay+1)) + temp_array = spread(p_lev(1,:), dim = 1, ncopies=ncol) + call move_alloc(temp_array, p_lev) + allocate(temp_array(ncol, nlay+1)) + temp_array = spread(t_lev(1,:), dim = 1, ncopies=ncol) + call move_alloc(temp_array, t_lev) + + ! This puts pressure and temperature arrays on the GPU + !$acc enter data copyin(p_lay, p_lev, t_lay, t_lev) + !$omp target enter data map(to:p_lay, p_lev, t_lay, t_lev) + end if ! ---------------------------------------------------------------------------- ! load data into classes call load_and_init(k_dist, k_dist_file, gas_concs) @@ -371,7 +394,15 @@ program rte_rrtmgp_clouds_aerosols end do ! call system_clock(finish_all, clock_rate) - ! + + ! Release GPU memory for p_lay, t_lay, p_lev, t_lev, q, o3) + if(do_rce) then + !$acc end data + !$omp target end data + else + !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) + !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) + end if if(do_clouds) then !$acc exit data delete( lwp, iwp, rel, rei) !$omp target exit data map(release:lwp, iwp, rel, rei) @@ -379,8 +410,6 @@ program rte_rrtmgp_clouds_aerosols ! ! TK - release aerosol memory ! - !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) - !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) #if defined(_OPENACC) || defined(_OPENMP) avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) @@ -409,6 +438,13 @@ program rte_rrtmgp_clouds_aerosols contains ! ---------------------------------------------------------------------------------- subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) + ! + ! Construct profiles of pressure, temperature, humidity, and ozone + ! more or less following the RCEMIP protocol for a surface temperature of 300K + ! more or less follows a Python implementation by Chiel van Heerwardeen + ! Extensions for future - variable SST and T profile, variable RH, lapse rate in stratosphere + ! will all access absorption coefficient data more realistically + ! real(wp), intent(in ) :: SST integer, intent(in ) :: ncol, nlay real(wp), dimension(ncol, nlay ), intent(out) :: p_lay, t_lay, q, o3 @@ -423,6 +459,9 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) real(wp), parameter :: g1 = 3.6478_wp, g2 = 0.83209_wp, g3 = 11.3515_wp, o3_min = 1e-13_wp ! According to CvH RRTMGP in Single Precision will fail with lower ozone concentrations + ! + ! Split resolution above and below RCE tropopause (15 km or about 125 hPa) + ! z_lev(:) = [0._wp, 2._wp* z_trop /nlay * [(i, i=1, nlay/2)], & z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) diff --git a/gas-optics/mo_gas_concentrations.F90 b/gas-optics/mo_gas_concentrations.F90 index 0ae9965b2..592baf57d 100644 --- a/gas-optics/mo_gas_concentrations.F90 +++ b/gas-optics/mo_gas_concentrations.F90 @@ -336,7 +336,6 @@ function get_vmr_1d(this, gas, array) result(error_msg) if(this%nlay > 0 .and. this%nlay /= size(array)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' - print *, "get_vmr_1d", this%nlay, size(array) end if if(error_msg /= "") return @@ -397,7 +396,6 @@ function get_vmr_2d(this, gas, array) result(error_msg) end if if(this%nlay > 0 .and. this%nlay /= size(array,2)) then error_msg = 'ty_gas_concs%get_vmr; gas ' // trim(gas) // ' array is wrong size (nlay)' - print *, "get_vmr_2d", this%nlay, size(array, 2) end if if(error_msg /= "") return @@ -409,7 +407,6 @@ function get_vmr_2d(this, gas, array) result(error_msg) !$omp target teams distribute parallel do simd do ilay = 1, size(array,2) do icol = 1, size(array,1) - !print *, (size(this%concs)) #ifdef _CRAYFTN array(icol,ilay) = p(icol,ilay) #else From 5c05c913e30fb633d50cde85a11696e5af8e7810 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 16 Jun 2023 22:48:39 -0400 Subject: [PATCH 13/38] Consistent OpenACC memory approach --- examples/all-sky/rrtmgp_allsky.F90 | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 69738ddce..5d8f2640f 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -150,8 +150,8 @@ program rte_rrtmgp_clouds_aerosols if(do_rce) then allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) allocate(q (ncol, nlay), o3(ncol, nlay)) - !$acc data create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + !$acc enter data create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target enter data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) call stop_on_err(gas_concs%init(gas_names)) @@ -396,13 +396,8 @@ program rte_rrtmgp_clouds_aerosols call system_clock(finish_all, clock_rate) ! Release GPU memory for p_lay, t_lay, p_lev, t_lev, q, o3) - if(do_rce) then - !$acc end data - !$omp target end data - else - !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) - !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) - end if + !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) + !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) if(do_clouds) then !$acc exit data delete( lwp, iwp, rel, rei) !$omp target exit data map(release:lwp, iwp, rel, rei) From 0f99ec2caa9a17d87c422416b3f4d4b381b3f752 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 16 Jun 2023 23:20:10 -0400 Subject: [PATCH 14/38] Elemental functions need acc diretives --- examples/all-sky/rrtmgp_allsky.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 5d8f2640f..5138fb090 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -461,8 +461,8 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) - !$acc data copyin(z_lev, lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target data map(to:z_lev, lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + !$acc data copyin(z_lev, z_lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target data map(to:z_lev, z_lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) !$acc parallel loop collapse(2) !$omp target teams distribute parallel do simd collapse(2) @@ -486,13 +486,18 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) end subroutine compute_profiles ! --------------------------------------- elemental function Tv(T, q) - real(wp), intent(in) :: T, q - real(wp) :: Tv + !$acc routine seq + !$omp declare target + + real(wp), intent(in) :: T, q + real(wp) :: Tv - Tv = (1. + 0.608*q) * T + Tv = (1. + 0.608*q) * T end function Tv ! --------------------------------------- elemental subroutine thermo(z, SST, p, T, q) + !$acc routine seq + !$omp declare target real(wp), intent(in) :: z, SST ! Height, SST real(wp), intent(out) :: p, T, q From 3b4bd1344ac696d6b157a2f56f89fa6af2e57080 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 14:24:43 -0400 Subject: [PATCH 15/38] Garand vestiges are removed; fluxes are written but not checked yet --- examples/all-sky/Makefile | 3 - examples/all-sky/rrtmgp_allsky.F90 | 195 ++++++++++++++++------------- 2 files changed, 105 insertions(+), 93 deletions(-) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index 79220e2af..f15b189f1 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -29,7 +29,6 @@ VPATH = ../:$(RRTMGP_ROOT)/rrtmgp-frontend # Needed for cloud_optics and aerosol # Extra sources -- extensions to RRTMGP classes, shared infrastructure, local sources # ADDITIONS = mo_load_coefficients.o mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.o -ADDITIONS += mo_garand_atmos_io.o ADDITIONS += mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.o # @@ -43,13 +42,11 @@ rrtmgp_allsky.o: $(ADDITIONS) rrtmgp_allsky.F90 mo_cloud_optics_rrtmgp.o: mo_cloud_optics_rrtmgp.F90 mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 -mo_garand_atmos_io.o: mo_simple_netcdf.o mo_garand_atmos_io.F90 mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 mo_aerosol_optics_rrtmgp_merra.o: mo_aerosol_optics_rrtmgp_merra.F90 mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.F90 tests: - cp ${RRTMGP_DATA}/examples/all-sky/inputs/garand-atmos-1.nc rrtmgp-allsky.nc $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 5138fb090..1f9fbfe0a 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -15,7 +15,6 @@ program rte_rrtmgp_clouds_aerosols only: load_cld_lutcoeff, load_cld_padecoeff use mo_load_aerosol_coefficients, & only: load_aero_lutcoeff - use mo_garand_atmos_io, only: read_atmos, write_lw_fluxes, write_sw_fluxes use mo_rte_config, only: rte_config_checks implicit none ! ---------------------------------------------------------------------------------- @@ -95,13 +94,10 @@ program rte_rrtmgp_clouds_aerosols character(len=8) :: char_input integer :: nUserArgs, nloops, ncol, nlay - logical :: write_fluxes = .false. + ! logical :: write_fluxes = .false. logical :: do_clouds = .false., use_luts = .true. logical :: do_aerosols = .false. - logical :: do_rce = .true. integer, parameter :: ngas = 8 -! character(len=3), dimension(ngas) & -! :: gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] character(len=256) :: input_file, k_dist_file, cloud_optics_file, aerosol_optics_file ! @@ -147,56 +143,21 @@ program rte_rrtmgp_clouds_aerosols end if if (nUserArgs > 7) print *, "Ignoring command line arguments beyond the first seven..." ! ----------------------------------------------------------------------------------- - if(do_rce) then - allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) - allocate(q (ncol, nlay), o3(ncol, nlay)) - !$acc enter data create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target enter data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) - call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) - - call stop_on_err(gas_concs%init(gas_names)) - call stop_on_err(gas_concs%set_vmr("h2o", q )) - call stop_on_err(gas_concs%set_vmr("o3", o3)) - call stop_on_err(gas_concs%set_vmr("co2", 348.e-6_wp)) - call stop_on_err(gas_concs%set_vmr("ch4", 1650.e-9_wp)) - call stop_on_err(gas_concs%set_vmr("n2o", 306.e-9_wp)) - call stop_on_err(gas_concs%set_vmr("n2", 0.7808_wp)) - call stop_on_err(gas_concs%set_vmr("o2", 0.2095_wp)) - call stop_on_err(gas_concs%set_vmr("co", 0._wp)) - else - ! - ! Read temperature, pressure, gas concentrations. - ! Arrays are allocated as they are read - ! - call read_atmos(input_file, & - p_lay, t_lay, p_lev, t_lev, & - gas_concs_garand, col_dry) - deallocate(col_dry) - nlay = size(p_lay, 2) - - ! For clouds we'll use the first column, repeated over and over - call stop_on_err(gas_concs%init(gas_names)) - do igas = 1, ngas - call vmr_2d_to_1d(gas_concs, gas_concs_garand, gas_names(igas), size(p_lay, 1), nlay) - end do - ! If we trusted in Fortran allocate-on-assign we could skip the temp_array here - allocate(temp_array(ncol, nlay)) - temp_array = spread(p_lay(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, p_lay) - allocate(temp_array(ncol, nlay)) - temp_array = spread(t_lay(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, t_lay) - allocate(temp_array(ncol, nlay+1)) - temp_array = spread(p_lev(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, p_lev) - allocate(temp_array(ncol, nlay+1)) - temp_array = spread(t_lev(1,:), dim = 1, ncopies=ncol) - call move_alloc(temp_array, t_lev) - - ! This puts pressure and temperature arrays on the GPU - !$acc enter data copyin(p_lay, p_lev, t_lay, t_lev) - !$omp target enter data map(to:p_lay, p_lev, t_lay, t_lev) - end if + allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) + allocate(q (ncol, nlay), o3(ncol, nlay)) + !$acc enter data create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target enter data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) + + call stop_on_err(gas_concs%init(gas_names)) + call stop_on_err(gas_concs%set_vmr("h2o", q )) + call stop_on_err(gas_concs%set_vmr("o3", o3)) + call stop_on_err(gas_concs%set_vmr("co2", 348.e-6_wp)) + call stop_on_err(gas_concs%set_vmr("ch4", 1650.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2o", 306.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2", 0.7808_wp)) + call stop_on_err(gas_concs%set_vmr("o2", 0.2095_wp)) + call stop_on_err(gas_concs%set_vmr("co", 0._wp)) ! ---------------------------------------------------------------------------- ! load data into classes call load_and_init(k_dist, k_dist_file, gas_concs) @@ -417,17 +378,13 @@ program rte_rrtmgp_clouds_aerosols print *, " - per column(ms):", (finish_all-start_all) / real(ncol*nloops) / (1.0e-3*clock_rate) #endif + if(.true.) call write_fluxes + if(is_lw) then - !$acc exit data copyout(flux_up, flux_dn) - !$omp target exit data map(from:flux_up, flux_dn) - if(write_fluxes) call write_lw_fluxes(input_file, flux_up, flux_dn) - !$acc exit data delete(t_sfc, emis_sfc) + !$acc exit data delete( t_sfc, emis_sfc) !$omp target exit data map(release:t_sfc, emis_sfc) else - !$acc exit data copyout(flux_up, flux_dn, flux_dir) - !$omp target exit data map(from:flux_up, flux_dn, flux_dir) - if(write_fluxes) call write_sw_fluxes(input_file, flux_up, flux_dn, flux_dir) - !$acc exit data delete(sfc_alb_dir, sfc_alb_dif, mu0) + !$acc exit data delete( sfc_alb_dir, sfc_alb_dif, mu0) !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if contains @@ -527,33 +484,6 @@ subroutine stop_on_err(error_msg) error stop 1 end if end subroutine stop_on_err - ! ---------------------------------------------------------------------------------- - subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) - use mo_gas_concentrations, only: ty_gas_concs - use mo_rte_kind, only: wp - - type(ty_gas_concs), intent(in) :: gas_concs_garand - type(ty_gas_concs), intent(inout) :: gas_concs - character(len=*), intent(in) :: name - integer, intent(in) :: sz1, sz2 - - real(wp) :: tmp(sz1, sz2), tmp_col(sz2) - - !$acc data create( tmp, tmp_col) - !$omp target data map(alloc:tmp, tmp_col) - print *, "vmr_2d_to_1d " // name - print *, " sizes", sz1, sz2 - call stop_on_err(gas_concs_garand%get_vmr(name, tmp)) - !$acc kernels - !$omp target - tmp_col(:) = tmp(1, :) - !$acc end kernels - !$omp end target - - call stop_on_err(gas_concs%set_vmr (name, tmp_col)) - !$acc end data - !$omp end target data - end subroutine vmr_2d_to_1d ! -------------------------------------------------------------------------------------- ! subroutine compute_clouds @@ -731,4 +661,89 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) enddo enddo end subroutine get_relhum + !-------------------------------------------------------------------------------------------------------------------- + subroutine write_fluxes + use netcdf + use mo_simple_netcdf + integer :: ncid, i, col_dim, lay_dim, lev_dim, varid + real(wp) :: vmr(ncol, nlay) + ! + ! Write fluxes - make this optional? + ! + + ! + ! Define dimensions + ! + if(nf90_create("rrmtgp_allsky_" // merge("lw", "sw", is_lw) // ".nc", NF90_CLOBBER, ncid) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't create file rrmtgp_allsky_" // merge("lw", "sw", is_lw) // ".nc") + + if(nf90_def_dim(ncid, "col", ncol, col_dim) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't define col dimension") + if(nf90_def_dim(ncid, "lay", nlay, lay_dim) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't define lay dimension") + if(nf90_def_dim(ncid, "lev", nlay+1, lev_dim) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't define lev dimension") + + ! + ! Define variables + ! + ! State + ! + if(nf90_def_var(ncid, "p_lev", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define p_lev variable") + if(nf90_def_var(ncid, "t_lev", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define t_lev variable") + if(nf90_def_var(ncid, "p_lay", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define p_lay variable") + if(nf90_def_var(ncid, "t_lay", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define t_lay variable") + if(nf90_def_var(ncid, "h2o", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define h2o variable") + if(nf90_def_var(ncid, "o3", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define o3 variable") + ! All the gases except h2o, o3 - write as attributes? + + ! + ! Fluxes - definitions + ! + if(nf90_def_var(ncid, "flux_up", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define flux_up variable") + if(nf90_def_var(ncid, "flux_dn", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define flux_dm variable") + if(.not. is_lw) then + if(nf90_def_var(ncid, "flux_dir", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define flux_dir variable") + end if + if(nf90_enddef(ncid) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't end redefinition??") + + ! + ! Write variables + ! + ! State - writing + !$acc exit data copyout( p_lev, t_lev, p_lay, t_lay) + !$omp target exit data map(from:p_lev, t_lev, p_lay, t_lay) + call stop_on_err(write_field(ncid, "p_lev", p_lev)) + call stop_on_err(write_field(ncid, "t_lev", t_lev)) + call stop_on_err(write_field(ncid, "p_lay", p_lay)) + call stop_on_err(write_field(ncid, "t_lay", t_lay)) + call stop_on_err(gas_concs%get_vmr("h2o", vmr)) + call stop_on_err(write_field(ncid, "h2o", vmr)) + call stop_on_err(gas_concs%get_vmr("o3", vmr)) + call stop_on_err(write_field(ncid, "o3", vmr)) + + ! Fluxes - writing + !$acc exit data copyout( flux_up, flux_dn) + !$omp target exit data map(from:flux_up, flux_dn) + call stop_on_err(write_field(ncid, "flux_up", flux_up)) + call stop_on_err(write_field(ncid, "flux_dn", flux_dn)) + if(.not. is_lw) then + !$acc exit data copyout( flux_dir) + !$omp target exit data map(from:flux_dir) + call stop_on_err(write_field(ncid, "flux_dir", flux_dir)) + end if + + ! Close netCDF + if(nf90_close(ncid) /= NF90_NOERR) call stop_on_err("rrtmgp_allsky: error closing file??") + end subroutine write_fluxes end program rte_rrtmgp_clouds_aerosols From 0c2879129fc602efc0b9d3168950d6b739f376ac Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 15:26:14 -0400 Subject: [PATCH 16/38] Remove elemental subroutine; update Makefile --- .github/workflows/containerized-ci.yml | 1 - examples/all-sky/Makefile | 11 +-- examples/all-sky/rrtmgp_allsky.F90 | 77 ++++++++++--------- examples/{all-sky => }/mo_garand_atmos_io.F90 | 0 4 files changed, 47 insertions(+), 42 deletions(-) rename examples/{all-sky => }/mo_garand_atmos_io.F90 (100%) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index b9e70ccdd..c9c3ce9de 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -53,7 +53,6 @@ jobs: FC: ${{ matrix.fortran-compiler }} FCFLAGS: ${{ matrix.fcflags }} # Make variables: - NCHOME: /dummy NFHOME: /opt/netcdf-fortran RRTMGP_ROOT: ${{ github.workspace }} RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index f15b189f1..228017e70 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -28,7 +28,8 @@ VPATH = ../:$(RRTMGP_ROOT)/rrtmgp-frontend # Needed for cloud_optics and aerosol # # Extra sources -- extensions to RRTMGP classes, shared infrastructure, local sources # -ADDITIONS = mo_load_coefficients.o mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.o +ADDITIONS = mo_load_coefficients.o mo_simple_netcdf.o +ADDITIONS += mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.o ADDITIONS += mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.o # @@ -40,10 +41,10 @@ rrtmgp_allsky: $(ADDITIONS) rrtmgp_allsky.o rrtmgp_allsky.o: $(ADDITIONS) rrtmgp_allsky.F90 -mo_cloud_optics_rrtmgp.o: mo_cloud_optics_rrtmgp.F90 -mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 -mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 +mo_cloud_optics_rrtmgp.o: mo_cloud_optics_rrtmgp.F90 mo_aerosol_optics_rrtmgp_merra.o: mo_aerosol_optics_rrtmgp_merra.F90 +mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 +mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.F90 tests: @@ -55,4 +56,4 @@ check: # python ./compare-to-reference.py clean: - -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod + -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod *.nc diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 1f9fbfe0a..a7efc62fe 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -389,7 +389,7 @@ program rte_rrtmgp_clouds_aerosols end if contains ! ---------------------------------------------------------------------------------- - subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) + subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, o3) ! ! Construct profiles of pressure, temperature, humidity, and ozone ! more or less following the RCEMIP protocol for a surface temperature of 300K @@ -399,11 +399,12 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) ! real(wp), intent(in ) :: SST integer, intent(in ) :: ncol, nlay - real(wp), dimension(ncol, nlay ), intent(out) :: p_lay, t_lay, q, o3 + real(wp), dimension(ncol, nlay ), intent(out) :: p_lay, t_lay, q_lay, o3 real(wp), dimension(ncol, nlay+1), intent(out) :: p_lev, t_lev real(wp) :: z_lay(nlay), z_lev(nlay+1) - real(wp) :: q_tmp, p_hpa + real(wp) :: z, q, T, p + real(wp) :: p_hpa integer :: icol, ilay, i real(wp), parameter :: z_trop = 15000._wp, z_top = 70.e3_wp @@ -411,6 +412,12 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) real(wp), parameter :: g1 = 3.6478_wp, g2 = 0.83209_wp, g3 = 11.3515_wp, o3_min = 1e-13_wp ! According to CvH RRTMGP in Single Precision will fail with lower ozone concentrations + real(wp), parameter :: g = 9.79764, Rd = 287.04, p0 = 101480. ! Surface pressure + real(wp), parameter :: z_q1 = 4.0e3, z_q2 = 7.5e3, q_t = 1.e-8 + real(wp), parameter :: gamma = 6.7e-3 + + real(wp), parameter :: q_0 = 0.01864 ! for 300 K SST. + ! ! Split resolution above and below RCE tropopause (15 km or about 125 hPa) ! @@ -418,28 +425,48 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) - !$acc data copyin(z_lev, z_lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target data map(to:z_lev, z_lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + !$acc enter data copyin(z_lev, z_lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target enter data map(to:z_lev, z_lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) !$acc parallel loop collapse(2) !$omp target teams distribute parallel do simd collapse(2) do ilay = 1, nlay do icol = 1, ncol - call thermo(z_lay(ilay), SST, p_lay(icol,ilay), t_lay(icol,ilay), q(icol,ilay)) + z = z_lay(ilay) + q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) + T = SST - gamma*z / (1. + 0.608*q) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) + if (z > z_trop) then + q = q_t + T = SST - gamma*z_trop/(1. + 0.608*q_0) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + end if + p_lay(icol,ilay) = p + t_lay(icol,ilay) = T + q_lay(icol,ilay) = q p_hpa = p_lay(icol,ilay) / 100._wp o3(icol, ilay) = max(o3_min, & g1 * p_hpa**g2 * exp(-p_hpa/g3) * 1.e-6_wp) - call thermo(z_lev(ilay), SST, p_lev(icol,ilay), t_lev(icol,ilay), q_tmp) end do end do - !$acc parallel loop - !$omp target teams distribute parallel do simd - do icol = 1, ncol - call thermo(z_lev(nlay+1), SST, p_lev(icol, nlay+1), t_lev(icol, nlay+1), q_tmp) - end do - !$acc end data - !$omp target end data + !$acc parallel loop collapse(2) + !$omp target teams distribute parallel do simd collapse(2) + do ilay = 1, nlay+1 + do icol = 1, ncol + z = z_lev(ilay) + q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) + T = SST - gamma*z / (1. + 0.608*q) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) + if (z > z_trop) then + q = q_t + T = SST - gamma*z_trop/(1. + 0.608*q_0) + p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + end if + p_lev(icol,ilay) = p + t_lev(icol,ilay) = T + end do + end do end subroutine compute_profiles ! --------------------------------------- elemental function Tv(T, q) @@ -451,28 +478,6 @@ elemental function Tv(T, q) Tv = (1. + 0.608*q) * T end function Tv - ! --------------------------------------- - elemental subroutine thermo(z, SST, p, T, q) - !$acc routine seq - !$omp declare target - real(wp), intent(in) :: z, SST ! Height, SST - real(wp), intent(out) :: p, T, q - - real(wp), parameter :: g = 9.79764, Rd = 287.04, p0 = 101480. ! Surface pressure - real(wp), parameter :: z_q1 = 4.0e3, z_q2 = 7.5e3, z_trop = 15.e3, q_t = 1.e-8 - real(wp), parameter :: gamma = 6.7e-3 - - real(wp), parameter :: q_0 = 0.01864 ! for 300 K SST. - - q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) - T = SST - gamma*z / (1. + 0.608*q) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) - if (z > z_trop) then - q = q_t - T = SST - gamma*z_trop/(1. + 0.608*q_0) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) - end if - end subroutine thermo ! ---------------------------------------------------------------------------------- subroutine stop_on_err(error_msg) use iso_fortran_env, only : error_unit diff --git a/examples/all-sky/mo_garand_atmos_io.F90 b/examples/mo_garand_atmos_io.F90 similarity index 100% rename from examples/all-sky/mo_garand_atmos_io.F90 rename to examples/mo_garand_atmos_io.F90 From 788808af3140e8b9ade0c60a69bf153afd1a9639 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 15:40:47 -0400 Subject: [PATCH 17/38] Remove elemental function --- examples/all-sky/rrtmgp_allsky.F90 | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index a7efc62fe..7eeb96ff6 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -404,7 +404,7 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, real(wp) :: z_lay(nlay), z_lev(nlay+1) real(wp) :: z, q, T, p - real(wp) :: p_hpa + real(wp) :: Tv, Tv0, p_hpa integer :: icol, ilay, i real(wp), parameter :: z_trop = 15000._wp, z_top = 70.e3_wp @@ -428,6 +428,10 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, !$acc enter data copyin(z_lev, z_lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) !$omp target enter data map(to:z_lev, z_lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + ! + ! The two loops are the same, except applied to layers and levels + ! but nvfortran doesn't seems to support elemental procedures in OpenACC loops + ! !$acc parallel loop collapse(2) !$omp target teams distribute parallel do simd collapse(2) do ilay = 1, nlay @@ -435,11 +439,13 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, z = z_lay(ilay) q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) T = SST - gamma*z / (1. + 0.608*q) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) + Tv = (1. + 0.608*q ) * T + Tv0 = (1. + 0.608*q_0) * SST + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) end if p_lay(icol,ilay) = p t_lay(icol,ilay) = T @@ -457,27 +463,19 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, z = z_lev(ilay) q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) T = SST - gamma*z / (1. + 0.608*q) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) + Tv = (1. + 0.608*q ) * T + Tv0 = (1. + 0.608*q_0) * SST + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) - p = p0 * (Tv(T, q)/Tv(SST, q_0))**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) end if p_lev(icol,ilay) = p t_lev(icol,ilay) = T end do end do end subroutine compute_profiles - ! --------------------------------------- - elemental function Tv(T, q) - !$acc routine seq - !$omp declare target - - real(wp), intent(in) :: T, q - real(wp) :: Tv - - Tv = (1. + 0.608*q) * T - end function Tv ! ---------------------------------------------------------------------------------- subroutine stop_on_err(error_msg) use iso_fortran_env, only : error_unit From bfffbd33c8af973cba878145b4469f905f2394c3 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 15:47:58 -0400 Subject: [PATCH 18/38] All the instances if you please --- examples/all-sky/rrtmgp_allsky.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 7eeb96ff6..2558785d1 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -445,7 +445,7 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) - p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) end if p_lay(icol,ilay) = p t_lay(icol,ilay) = T @@ -469,7 +469,7 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) - p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv(T, q))) ) + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) end if p_lev(icol,ilay) = p t_lev(icol,ilay) = T From 982fd9152e7cc1e6e17fb624afe770dea8bb06c4 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 16:00:46 -0400 Subject: [PATCH 19/38] Recompute Tv in stratosphere --- examples/all-sky/rrtmgp_allsky.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 2558785d1..f10d1b0d4 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -445,6 +445,7 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) + Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) end if p_lay(icol,ilay) = p @@ -469,6 +470,7 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) + Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) end if p_lev(icol,ilay) = p From 1f8313ea795a58c72ae7dd0ea4a6293971243632 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 20 Jun 2023 23:04:08 -0400 Subject: [PATCH 20/38] Create temps, pressures on the device only once. --- examples/all-sky/rrtmgp_allsky.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index f10d1b0d4..fe3f10a23 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -116,7 +116,7 @@ program rte_rrtmgp_clouds_aerosols ! nUserArgs = command_argument_count() - if (nUserArgs < 5) call stop_on_err("Need to supply ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") + if (nUserArgs < 5) call stop_on_err("Usage: rrtmgp_allsky ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") call get_command_argument(1, char_input) read(char_input, '(i8)') ncol @@ -425,8 +425,8 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) - !$acc enter data copyin(z_lev, z_lay) create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target enter data map(to:z_lev, z_lay) map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + !$acc data copyin(z_lev, z_lay) + !$omp target data map(to:z_lev, z_lay) ! ! The two loops are the same, except applied to layers and levels @@ -477,6 +477,8 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, t_lev(icol,ilay) = T end do end do + !$acc end data + !$omp target end data end subroutine compute_profiles ! ---------------------------------------------------------------------------------- subroutine stop_on_err(error_msg) From 62cbf43c04fc840eeff67548c192c38a33934e32 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 21 Jun 2023 12:40:58 -0400 Subject: [PATCH 21/38] Checking of all-sky reference results --- .github/workflows/containerized-ci.yml | 1 + .github/workflows/continuous-integration.yml | 1 + .github/workflows/self-hosted-ci.yml | 1 + examples/all-sky/Makefile | 8 +++---- examples/all-sky/compare-to-reference.py | 10 ++++----- examples/all-sky/rrtmgp_allsky.F90 | 23 +++++++++++--------- 6 files changed, 25 insertions(+), 19 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index c9c3ce9de..8b9fdc028 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -75,6 +75,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-timing # # Cache RFMIP files # diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 928bbb09d..714899782 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -47,6 +47,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-timing # # Synchronize the package index # diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 0b0e1c90c..f07dfd6f9 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -62,6 +62,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-timing # # Finalize build environment # diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index 228017e70..8cc1b4885 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -43,8 +43,8 @@ rrtmgp_allsky.o: $(ADDITIONS) rrtmgp_allsky.F90 mo_cloud_optics_rrtmgp.o: mo_cloud_optics_rrtmgp.F90 mo_aerosol_optics_rrtmgp_merra.o: mo_aerosol_optics_rrtmgp_merra.F90 -mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 -mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 +mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 +mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.F90 tests: @@ -52,8 +52,8 @@ tests: $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc check: - echo "No checking of all sky results at this time" - # python ./compare-to-reference.py + python ./compare-to-reference.py --allsky_file rrtmgp-allsky-lw.nc + python ./compare-to-reference.py --allsky_file rrtmgp-allsky-sw.nc clean: -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod *.nc diff --git a/examples/all-sky/compare-to-reference.py b/examples/all-sky/compare-to-reference.py index e57a42489..54923a6f9 100644 --- a/examples/all-sky/compare-to-reference.py +++ b/examples/all-sky/compare-to-reference.py @@ -16,7 +16,7 @@ parser = argparse.ArgumentParser( description="Compares all-sky example output to file in reference " "directory") - parser.add_argument("--allsky_file", type=str, default="rrtmgp-allsky.nc", + parser.add_argument("--allsky_file", type=str, default="rrtmgp-allsky-lw.nc", dest="file", help="Name of file inputs and outputs for " "all-sky problem (same for test and reference)") @@ -35,13 +35,13 @@ ref = xr.open_dataset(ref_file) failed = False - for v in ['lw_flux_up', 'lw_flux_dn', 'sw_flux_up', 'sw_flux_dn', - 'sw_flux_dir']: + for v in tst.variables: + if np.any(np.isnan(ref.variables[v].values)): + raise Exception(v + ": some ref values are missing. Now that is strange.") if np.all(np.isnan(tst.variables[v].values)): raise Exception("All test values are missing. Were the tests run?") if np.any(np.isnan(tst.variables[v].values)): - raise Exception( - "Some test values are missing. Now that is strange.") + raise Exception(v + ":Some test values are missing. Now that is strange.") # express as (tst-ref).variables[v].values when replacing reference file # to have same number of columns diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index fe3f10a23..1de5737c1 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -357,15 +357,10 @@ program rte_rrtmgp_clouds_aerosols call system_clock(finish_all, clock_rate) ! Release GPU memory for p_lay, t_lay, p_lev, t_lev, q, o3) - !$acc exit data delete(p_lay, p_lev, t_lay, t_lev) + !$acc exit data delete( p_lay, p_lev, t_lay, t_lev) !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) - if(do_clouds) then - !$acc exit data delete( lwp, iwp, rel, rei) - !$omp target exit data map(release:lwp, iwp, rel, rei) - end if - ! - ! TK - release aerosol memory - ! + + ! Memory for clouds and aerosols is released in write-fluxes #if defined(_OPENACC) || defined(_OPENMP) avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) @@ -681,8 +676,8 @@ subroutine write_fluxes ! ! Define dimensions ! - if(nf90_create("rrmtgp_allsky_" // merge("lw", "sw", is_lw) // ".nc", NF90_CLOBBER, ncid) /= NF90_NOERR) & - call stop_on_err("rrtmgp_allsky: can't create file rrmtgp_allsky_" // merge("lw", "sw", is_lw) // ".nc") + if(nf90_create("rrtmgp-allsky-" // merge("lw", "sw", is_lw) // ".nc", NF90_CLOBBER, ncid) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't create file rrmtgp-allsky-" // merge("lw", "sw", is_lw) // ".nc") if(nf90_def_dim(ncid, "col", ncol, col_dim) /= NF90_NOERR) & call stop_on_err("rrtmgp_allsky: can't define col dimension") @@ -739,6 +734,14 @@ subroutine write_fluxes call stop_on_err(gas_concs%get_vmr("o3", vmr)) call stop_on_err(write_field(ncid, "o3", vmr)) + if(do_clouds) then + !$acc exit data delete( lwp, iwp, rel, rei) + !$omp target exit data map(release:lwp, iwp, rel, rei) + end if + + if(do_aerosols) then + end if + ! Fluxes - writing !$acc exit data copyout( flux_up, flux_dn) !$omp target exit data map(from:flux_up, flux_dn) From 4ce5523b9152b9d185c8773b523e3034334f32fa Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 21 Jun 2023 13:55:01 -0400 Subject: [PATCH 22/38] Include cloud and aerosol state in all-sky outputs; add no-aerosols case --- examples/all-sky/Makefile | 12 ++++- examples/all-sky/rrtmgp_allsky.F90 | 78 ++++++++++++++++++++---------- 2 files changed, 63 insertions(+), 27 deletions(-) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index 8cc1b4885..4cdb3207b 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -48,12 +48,20 @@ mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.F90 tests: - $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc - $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw.nc \ + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw.nc \ + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw-no-aerosols.nc \ + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc + $(RUN_CMD) ./rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw-no-aerosols.nc \ + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc check: python ./compare-to-reference.py --allsky_file rrtmgp-allsky-lw.nc python ./compare-to-reference.py --allsky_file rrtmgp-allsky-sw.nc + python ./compare-to-reference.py --allsky_file rrtmgp-allsky-lw-no-aerosols.nc + python ./compare-to-reference.py --allsky_file rrtmgp-allsky-sw-no-aerosols.nc clean: -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod *.nc diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 1de5737c1..c6340617a 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -99,7 +99,7 @@ program rte_rrtmgp_clouds_aerosols logical :: do_aerosols = .false. integer, parameter :: ngas = 8 - character(len=256) :: input_file, k_dist_file, cloud_optics_file, aerosol_optics_file + character(len=256) :: output_file, k_dist_file, cloud_optics_file, aerosol_optics_file ! ! Timing variables ! @@ -116,7 +116,7 @@ program rte_rrtmgp_clouds_aerosols ! nUserArgs = command_argument_count() - if (nUserArgs < 5) call stop_on_err("Usage: rrtmgp_allsky ncol nlay nreps input_file gas-optics [cloud-optics [aerosol-optics]]") + if (nUserArgs < 5) call stop_on_err("Usage: rrtmgp_allsky ncol nlay nreps output_file gas-optics [cloud-optics [aerosol-optics]]") call get_command_argument(1, char_input) read(char_input, '(i8)') ncol @@ -130,7 +130,7 @@ program rte_rrtmgp_clouds_aerosols read(char_input, '(i8)') nloops if(nloops <= 0) call stop_on_err("Specify positive nreps (number of times to do ncol examples.") - call get_command_argument(4,input_file) + call get_command_argument(4,output_file) call get_command_argument(5,k_dist_file) if (nUserArgs >= 6) then @@ -666,7 +666,7 @@ end subroutine get_relhum !-------------------------------------------------------------------------------------------------------------------- subroutine write_fluxes use netcdf - use mo_simple_netcdf + use mo_simple_netcdf, only: write_field integer :: ncid, i, col_dim, lay_dim, lev_dim, varid real(wp) :: vmr(ncol, nlay) ! @@ -676,8 +676,8 @@ subroutine write_fluxes ! ! Define dimensions ! - if(nf90_create("rrtmgp-allsky-" // merge("lw", "sw", is_lw) // ".nc", NF90_CLOBBER, ncid) /= NF90_NOERR) & - call stop_on_err("rrtmgp_allsky: can't create file rrmtgp-allsky-" // merge("lw", "sw", is_lw) // ".nc") + if(nf90_create(trim(output_file), NF90_CLOBBER, ncid) /= NF90_NOERR) & + call stop_on_err("rrtmgp_allsky: can't create file " // trim(output_file)) if(nf90_def_dim(ncid, "col", ncol, col_dim) /= NF90_NOERR) & call stop_on_err("rrtmgp_allsky: can't define col dimension") @@ -691,30 +691,34 @@ subroutine write_fluxes ! ! State ! - if(nf90_def_var(ncid, "p_lev", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define p_lev variable") - if(nf90_def_var(ncid, "t_lev", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define t_lev variable") - if(nf90_def_var(ncid, "p_lay", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define p_lay variable") - if(nf90_def_var(ncid, "t_lay", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define t_lay variable") - if(nf90_def_var(ncid, "h2o", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define h2o variable") - if(nf90_def_var(ncid, "o3", NF90_DOUBLE, [col_dim, lay_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define o3 variable") - ! All the gases except h2o, o3 - write as attributes? + call create_var("p_lev", ncid, [col_dim, lev_dim]) + call create_var("t_lev", ncid, [col_dim, lev_dim]) + call create_var("p_lay", ncid, [col_dim, lay_dim]) + call create_var("t_lay", ncid, [col_dim, lay_dim]) + call create_var("h2o", ncid, [col_dim, lay_dim]) + call create_var("o3", ncid, [col_dim, lay_dim]) + ! All the gases except h2o, o3 - write as attributes? Or not bother? + + if(do_clouds) then + call create_var("lwp", ncid, [col_dim, lay_dim]) + call create_var("iwp", ncid, [col_dim, lay_dim]) + call create_var("rel", ncid, [col_dim, lay_dim]) + call create_var("rei", ncid, [col_dim, lay_dim]) + end if + if(do_aerosols) then + if(nf90_def_var(ncid, "aero_type", NF90_SHORT, [col_dim, lay_dim], varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define variable aero_type") + call create_var("aero_size", ncid, [col_dim, lay_dim]) + call create_var("aero_mass", ncid, [col_dim, lay_dim]) + end if ! ! Fluxes - definitions ! - if(nf90_def_var(ncid, "flux_up", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define flux_up variable") - if(nf90_def_var(ncid, "flux_dn", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define flux_dm variable") + call create_var("flux_up", ncid, [col_dim, lev_dim]) + call create_var("flux_dn", ncid, [col_dim, lev_dim]) if(.not. is_lw) then - if(nf90_def_var(ncid, "flux_dir", NF90_DOUBLE, [col_dim, lev_dim], varid) /= NF90_NOERR) & - call stop_on_err("create_var: can't define flux_dir variable") + call create_var("flux_dir", ncid, [col_dim, lev_dim]) end if if(nf90_enddef(ncid) /= NF90_NOERR) & call stop_on_err("rrtmgp_allsky: can't end redefinition??") @@ -735,11 +739,22 @@ subroutine write_fluxes call stop_on_err(write_field(ncid, "o3", vmr)) if(do_clouds) then + !$acc exit data copyout( lwp, iwp, rel, rei) + !$omp target exit data map(from:lwp, iwp, rel, rei) + call stop_on_err(write_field(ncid, "lwp", lwp)) + call stop_on_err(write_field(ncid, "iwp", iwp)) + call stop_on_err(write_field(ncid, "rel", rel)) + call stop_on_err(write_field(ncid, "rei", rei)) !$acc exit data delete( lwp, iwp, rel, rei) !$omp target exit data map(release:lwp, iwp, rel, rei) end if if(do_aerosols) then + !$acc exit data copyout( aero_size, aero_mass, aero_type) + !$omp target exit data map(from:aero_size, aero_mass, aero_type) + call stop_on_err(write_field(ncid, "aero_size", aero_size)) + call stop_on_err(write_field(ncid, "aero_mass", aero_mass)) + call stop_on_err(write_field(ncid, "aero_type", aero_type)) end if ! Fluxes - writing @@ -756,4 +771,17 @@ subroutine write_fluxes ! Close netCDF if(nf90_close(ncid) /= NF90_NOERR) call stop_on_err("rrtmgp_allsky: error closing file??") end subroutine write_fluxes + ! --------------------------------------------------------- + subroutine create_var(name, ncid, dim_ids) + use netcdf + character(len=*), intent(in) :: name + integer, intent(in) :: ncid + integer, dimension(:), intent(in) :: dim_ids + + integer :: varid + + if(nf90_def_var(ncid, trim(name), NF90_DOUBLE, dim_ids, varid) /= NF90_NOERR) & + call stop_on_err("create_var: can't define " // trim(name) // " variable") + end subroutine create_var + ! --------------------------------------------------------- end program rte_rrtmgp_clouds_aerosols From 2caba63103a2faa1b54882e39768a2e1f765db9a Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 22 Jun 2023 21:27:49 -0400 Subject: [PATCH 23/38] Accelerator data management in all-sky example --- examples/all-sky/rrtmgp_allsky.F90 | 86 +++++++++++++++++------------- 1 file changed, 48 insertions(+), 38 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index c6340617a..e338044e5 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -1,4 +1,4 @@ -program rte_rrtmgp_clouds_aerosols +program rte_rrtmgp_allsky use mo_rte_kind, only: wp, i8, wl use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str @@ -145,8 +145,8 @@ program rte_rrtmgp_clouds_aerosols ! ----------------------------------------------------------------------------------- allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) allocate(q (ncol, nlay), o3(ncol, nlay)) - !$acc enter data create( p_lay, t_lay, p_lev, t_lev, q, o3) - !$omp target enter data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) + !$acc data create( p_lay, t_lay, p_lev, t_lev, q, o3) + !$omp target data map(alloc:p_lay, t_lay, p_lev, t_lev, q, o3) call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) call stop_on_err(gas_concs%init(gas_names)) @@ -214,7 +214,7 @@ program rte_rrtmgp_clouds_aerosols !$acc enter data copyin(atmos) create(atmos%tau, atmos%ssa, atmos%g) !$omp target enter data map(alloc:atmos%tau, atmos%ssa, atmos%g) class default - call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") + call stop_on_err("rte_rrtmgp_allsky: Don't recognize the kind of optical properties ") end select ! ---------------------------------------------------------------------------- @@ -227,8 +227,8 @@ program rte_rrtmgp_clouds_aerosols !!$omp end parallel ! allocate(sfc_alb_dir(nbnd, ncol), sfc_alb_dif(nbnd, ncol), mu0(ncol)) - !$acc enter data create(sfc_alb_dir, sfc_alb_dif, mu0) - !$omp target enter data map(alloc:sfc_alb_dir, sfc_alb_dif, mu0) + !$acc enter data create( sfc_alb_dir, sfc_alb_dif, mu0) + !$omp target enter data map(alloc:sfc_alb_dir, sfc_alb_dif, mu0) ! Ocean-ish values for no particular reason !$acc kernels !$omp target @@ -244,8 +244,8 @@ program rte_rrtmgp_clouds_aerosols !!$omp end parallel allocate(t_sfc(ncol), emis_sfc(nbnd, ncol)) - !$acc enter data create(t_sfc, emis_sfc) - !$omp target enter data map(alloc:t_sfc, emis_sfc) + !$acc enter data create (t_sfc, emis_sfc) + !$omp target enter data map(alloc:t_sfc, emis_sfc) ! Surface temperature !$acc kernels !$omp target @@ -262,8 +262,8 @@ program rte_rrtmgp_clouds_aerosols allocate(flux_up(ncol,nlay+1), flux_dn(ncol,nlay+1)) !!$omp end parallel - !$acc enter data create(flux_up, flux_dn) - !$omp target enter data map(alloc:flux_up, flux_dn) + !$acc data create( flux_up, flux_dn) + !$omp target data map(alloc:flux_up, flux_dn) if(is_sw) then allocate(flux_dir(ncol,nlay+1)) !$acc enter data create(flux_dir) @@ -305,6 +305,9 @@ program rte_rrtmgp_clouds_aerosols fluxes%flux_up => flux_up(:,:) fluxes%flux_dn => flux_dn(:,:) if(is_lw) then + ! + ! Should we allocate these once, rather than once per loop? They're big. + ! !$acc data create( lw_sources, lw_sources%lay_source, lw_sources%lev_source_inc) & !$acc create( lw_sources%lev_source_dec, lw_sources%sfc_source, lw_sources%sfc_source_Jac) !$omp target data map(alloc: lw_sources%lay_source, lw_sources%lev_source_inc) & @@ -356,12 +359,6 @@ program rte_rrtmgp_clouds_aerosols ! call system_clock(finish_all, clock_rate) - ! Release GPU memory for p_lay, t_lay, p_lev, t_lev, q, o3) - !$acc exit data delete( p_lay, p_lev, t_lay, t_lev) - !$omp target exit data map(release:p_lay, p_lev, t_lay, t_lev) - - ! Memory for clouds and aerosols is released in write-fluxes - #if defined(_OPENACC) || defined(_OPENMP) avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) @@ -382,6 +379,15 @@ program rte_rrtmgp_clouds_aerosols !$acc exit data delete( sfc_alb_dir, sfc_alb_dif, mu0) !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if + + ! Clouds and aerosols used enter data - we won't release explicitly + + ! fluxes - but not flux_dir, which used enter data + !$acc end data + !$omp target end data + ! p_lay etc + !$acc end data + !$omp target end data contains ! ---------------------------------------------------------------------------------- subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, o3) @@ -511,7 +517,7 @@ subroutine compute_clouds !$acc enter data copyin(clouds) create(clouds%tau, clouds%ssa, clouds%g) !$omp target enter data map(alloc:clouds%tau, clouds%ssa, clouds%g) class default - call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") + call stop_on_err("rte_rrtmgp_allsky: Don't recognize the kind of optical properties ") end select ! ! Cloud physical properties @@ -565,28 +571,31 @@ subroutine compute_aerosols select type(aerosols) class is (ty_optical_props_1scl) call stop_on_err(aerosols%alloc_1scl(ncol, nlay)) - !$acc enter data copyin(aerosols) create(aerosols%tau) - !$omp target enter data map(alloc:aerosols%tau) + !$acc enter data copyin(aerosols) create(aerosols%tau) + !$omp target enter data map (alloc:aerosols%tau) class is (ty_optical_props_2str) call stop_on_err(aerosols%alloc_2str(ncol, nlay)) - !$acc enter data copyin(aerosols) create(aerosols%tau, aerosols%ssa, aerosols%g) - !$omp target enter data map(alloc:aerosols%tau, aerosols%ssa, aerosols%g) + !$acc enter data copyin(aerosols) create(aerosols%tau, aerosols%ssa, aerosols%g) + !$omp target enter data map (alloc:aerosols%tau, aerosols%ssa, aerosols%g) class default - call stop_on_err("rte_rrtmgp_clouds_aerosols: Don't recognize the kind of optical properties ") + call stop_on_err("rte_rrtmgp_allsky: Don't recognize the kind of optical properties ") end select ! ! Derive relative humidity from profile + ! relhum is on the CPU at this point ! call stop_on_err(gas_concs%get_vmr("h2o",vmr_h2o)) allocate(relhum(ncol, nlay)) + !$acc enter data create( relhum) + !$omp target enter data map(alloc:relhum) call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! ! Aerosol properties ! allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & - aero_mass(ncol,nlay), aero_mask(ncol,nlay)) - !$acc enter data create( aero_mask, aero_type, aero_size, aero_mass) - !$omp target enter data map(alloc:aero_mask, aero_type, aero_size, aero_mass) + aero_mass(ncol,nlay)) + !$acc enter data create( aero_type, aero_size, aero_mass) + !$omp target enter data map(alloc:aero_type, aero_size, aero_mass) ! Restrict sulfate aerosols to lower stratosphere (> 50 hPa = 50*100 Pa; < 100 hPa = 100*100 Pa) ! and dust aerosols to the lower troposphere (> 700 hPa; < 900 hPa), and @@ -650,6 +659,8 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! ------------------- ! Derive layer virtual temperature + !$acc parallel loop collapse(2) copyin(p_lay, vmr_h2o, t_lay) copyout( relhum) + !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay, vmr_h2o, t_lay) map(from:relhum) do i = 1, ncol do k = 1, nlay ! Convert h2o vmr to mmr @@ -727,44 +738,43 @@ subroutine write_fluxes ! Write variables ! ! State - writing - !$acc exit data copyout( p_lev, t_lev, p_lay, t_lay) - !$omp target exit data map(from:p_lev, t_lev, p_lay, t_lay) + !$acc update host ( p_lev, t_lev, p_lay, t_lay) + !$omp target update map(from:p_lev, t_lev, p_lay, t_lay) call stop_on_err(write_field(ncid, "p_lev", p_lev)) call stop_on_err(write_field(ncid, "t_lev", t_lev)) call stop_on_err(write_field(ncid, "p_lay", p_lay)) call stop_on_err(write_field(ncid, "t_lay", t_lay)) + ! Array vmr is on the host, not the device, but is copied-out call stop_on_err(gas_concs%get_vmr("h2o", vmr)) call stop_on_err(write_field(ncid, "h2o", vmr)) call stop_on_err(gas_concs%get_vmr("o3", vmr)) call stop_on_err(write_field(ncid, "o3", vmr)) if(do_clouds) then - !$acc exit data copyout( lwp, iwp, rel, rei) - !$omp target exit data map(from:lwp, iwp, rel, rei) + !$acc update host( lwp, iwp, rel, rei) + !$omp target update map(from:lwp, iwp, rel, rei) call stop_on_err(write_field(ncid, "lwp", lwp)) call stop_on_err(write_field(ncid, "iwp", iwp)) call stop_on_err(write_field(ncid, "rel", rel)) call stop_on_err(write_field(ncid, "rei", rei)) - !$acc exit data delete( lwp, iwp, rel, rei) - !$omp target exit data map(release:lwp, iwp, rel, rei) end if if(do_aerosols) then - !$acc exit data copyout( aero_size, aero_mass, aero_type) - !$omp target exit data map(from:aero_size, aero_mass, aero_type) + !$acc update host( aero_size, aero_mass, aero_type) + !$omp target update map(from:aero_size, aero_mass, aero_type) call stop_on_err(write_field(ncid, "aero_size", aero_size)) call stop_on_err(write_field(ncid, "aero_mass", aero_mass)) call stop_on_err(write_field(ncid, "aero_type", aero_type)) end if ! Fluxes - writing - !$acc exit data copyout( flux_up, flux_dn) - !$omp target exit data map(from:flux_up, flux_dn) + !$acc update host( flux_up, flux_dn) + !$omp target update map(from:flux_up, flux_dn) call stop_on_err(write_field(ncid, "flux_up", flux_up)) call stop_on_err(write_field(ncid, "flux_dn", flux_dn)) if(.not. is_lw) then - !$acc exit data copyout( flux_dir) - !$omp target exit data map(from:flux_dir) + !$acc update host( flux_dir) + !$omp target update map(from:flux_dir) call stop_on_err(write_field(ncid, "flux_dir", flux_dir)) end if @@ -784,4 +794,4 @@ subroutine create_var(name, ncid, dim_ids) call stop_on_err("create_var: can't define " // trim(name) // " variable") end subroutine create_var ! --------------------------------------------------------- -end program rte_rrtmgp_clouds_aerosols +end program rte_rrtmgp_allsky From 384beab97a08b4e6be103fa5f85999b3055d7e13 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 23 Jun 2023 14:51:58 -0400 Subject: [PATCH 24/38] omp directive syntax --- examples/all-sky/rrtmgp_allsky.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index e338044e5..6459032cd 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -383,11 +383,11 @@ program rte_rrtmgp_allsky ! Clouds and aerosols used enter data - we won't release explicitly ! fluxes - but not flux_dir, which used enter data - !$acc end data - !$omp target end data + !$acc end data + !$omp end target data ! p_lay etc - !$acc end data - !$omp target end data + !$acc end data + !$omp end target data contains ! ---------------------------------------------------------------------------------- subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, o3) @@ -478,8 +478,8 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, t_lev(icol,ilay) = T end do end do - !$acc end data - !$omp target end data + !$acc end data + !$omp end target data end subroutine compute_profiles ! ---------------------------------------------------------------------------------- subroutine stop_on_err(error_msg) From 808121bb2a9cfdbe249c331dbe2ac4dce189560e Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 23 Jun 2023 15:58:52 -0400 Subject: [PATCH 25/38] Syntax, relhum on GPU for compute_aerosols --- examples/all-sky/rrtmgp_allsky.F90 | 51 +++++++++++++++--------------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 6459032cd..bbc69eeae 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -557,9 +557,9 @@ end subroutine compute_clouds ! -------------------------------------------------------------------------------------- ! subroutine compute_aerosols - real(wp), dimension(ncol,nlay) :: vmr_h2o ! h2o vmr - logical :: is_sulfate, is_dust, is_even_column - ! + real(wp), dimension(ncol,nlay) :: vmr_h2o ! h2o vmr + logical :: is_sulfate, is_dust, is_even_column + ! ! Variable and memory allocation ! if(is_sw) then @@ -582,20 +582,21 @@ subroutine compute_aerosols end select ! ! Derive relative humidity from profile - ! relhum is on the CPU at this point + ! Keep vmr_h2o on the GPU ! + !$acc data create( vmr_h2o) + !$omp target data map(alloc:vmr_h2o) call stop_on_err(gas_concs%get_vmr("h2o",vmr_h2o)) - allocate(relhum(ncol, nlay)) - !$acc enter data create( relhum) - !$omp target enter data map(alloc:relhum) - call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! ! Aerosol properties ! allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & - aero_mass(ncol,nlay)) - !$acc enter data create( aero_type, aero_size, aero_mass) - !$omp target enter data map(alloc:aero_type, aero_size, aero_mass) + aero_mass(ncol,nlay), relhum (ncol,nlay)) + !$acc enter data create( aero_type, aero_size, aero_mass, relhum) + !$omp target enter data map(alloc:aero_type, aero_size, aero_mass, relhum) + call get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) + !$acc end data + !$omp end target data ! Restrict sulfate aerosols to lower stratosphere (> 50 hPa = 50*100 Pa; < 100 hPa = 100*100 Pa) ! and dust aerosols to the lower troposphere (> 700 hPa; < 900 hPa), and @@ -647,14 +648,14 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! Local variables integer :: i, k - real(wp) :: mmr_h2o ! water mass mixing ratio - real(wp) :: q_lay ! water specific humidity + real(wp) :: mmr_h2o ! water mass mixing ratio + real(wp) :: q_lay ! water specific humidity real(wp) :: q_lay_min, q_tmp, es_tmp real(wp) :: mwd, t_ref, rh ! Set constants - mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights - t_ref = 273.16_wp ! reference temperature (K) + mwd = m_h2o/m_dry ! ratio of water to dry air molecular weights + t_ref = 273.16_wp ! reference temperature (K) q_lay_min = 1.e-7_wp ! minimum water mass mixing ratio ! ------------------- @@ -738,8 +739,8 @@ subroutine write_fluxes ! Write variables ! ! State - writing - !$acc update host ( p_lev, t_lev, p_lay, t_lay) - !$omp target update map(from:p_lev, t_lev, p_lay, t_lay) + !$acc update host(p_lev, t_lev, p_lay, t_lay) + !$omp target update from(p_lev, t_lev, p_lay, t_lay) call stop_on_err(write_field(ncid, "p_lev", p_lev)) call stop_on_err(write_field(ncid, "t_lev", t_lev)) call stop_on_err(write_field(ncid, "p_lay", p_lay)) @@ -751,8 +752,8 @@ subroutine write_fluxes call stop_on_err(write_field(ncid, "o3", vmr)) if(do_clouds) then - !$acc update host( lwp, iwp, rel, rei) - !$omp target update map(from:lwp, iwp, rel, rei) + !$acc update host(lwp, iwp, rel, rei) + !$omp target update from(lwp, iwp, rel, rei) call stop_on_err(write_field(ncid, "lwp", lwp)) call stop_on_err(write_field(ncid, "iwp", iwp)) call stop_on_err(write_field(ncid, "rel", rel)) @@ -760,21 +761,21 @@ subroutine write_fluxes end if if(do_aerosols) then - !$acc update host( aero_size, aero_mass, aero_type) - !$omp target update map(from:aero_size, aero_mass, aero_type) + !$acc update host(aero_size, aero_mass, aero_type) + !$omp target update from(aero_size, aero_mass, aero_type) call stop_on_err(write_field(ncid, "aero_size", aero_size)) call stop_on_err(write_field(ncid, "aero_mass", aero_mass)) call stop_on_err(write_field(ncid, "aero_type", aero_type)) end if ! Fluxes - writing - !$acc update host( flux_up, flux_dn) - !$omp target update map(from:flux_up, flux_dn) + !$acc update host(flux_up, flux_dn) + !$omp target update from(flux_up, flux_dn) call stop_on_err(write_field(ncid, "flux_up", flux_up)) call stop_on_err(write_field(ncid, "flux_dn", flux_dn)) if(.not. is_lw) then - !$acc update host( flux_dir) - !$omp target update map(from:flux_dir) + !$acc update host(flux_dir) + !$omp target update from(flux_dir) call stop_on_err(write_field(ncid, "flux_dir", flux_dir)) end if From a1c0d7d7d9507a5242f1033dc31b00affb42178e Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 23 Jun 2023 17:06:45 -0400 Subject: [PATCH 26/38] Revised timing output from all-sky example --- .github/workflows/self-hosted-ci.yml | 2 -- examples/all-sky/rrtmgp_allsky.F90 | 23 ++++++++++++++--------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index f07dfd6f9..e7cf0bbab 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -43,8 +43,6 @@ jobs: FC: ftn FCFLAGS: ${{ matrix.fcflags }} # Make variables: - NCHOME: /dummy - NFHOME: /dummy RRTMGP_ROOT: ${{ github.workspace }} RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data RTE_KERNELS: ${{ matrix.rte-kernels }} diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index bbc69eeae..c107de774 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -1,4 +1,6 @@ program rte_rrtmgp_allsky + use, intrinsic :: iso_fortran_env, & + only: output_unit use mo_rte_kind, only: wp, i8, wl use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str @@ -104,7 +106,7 @@ program rte_rrtmgp_allsky ! Timing variables ! integer(kind=i8) :: start, finish, start_all, finish_all, clock_rate - real(wp) :: avg + real(wp) :: avg, mint integer(kind=i8), allocatable :: elapsed(:) ! NAR OpenMP CPU directives in compatible with OpenMP GPU directives !!$omp threadprivate( lw_sources, toa_flux, flux_up, flux_dn, flux_dir ) @@ -360,15 +362,18 @@ program rte_rrtmgp_allsky call system_clock(finish_all, clock_rate) #if defined(_OPENACC) || defined(_OPENMP) - avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) - - print *, "Execution times - min(s) :", minval(elapsed) / real(clock_rate) - print *, " - avg(s) :", avg / real(clock_rate) - print *, " - per column(ms):", avg / real(ncol) / (1.0e-3*clock_rate) + avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) + mint = minval(elapsed) #else - print *, "Execution times - total(s) :", (finish_all-start_all) / real(clock_rate) - print *, " - per column(ms):", (finish_all-start_all) / real(ncol*nloops) / (1.0e-3*clock_rate) -#endif + avg = (finish_all-start_all) + mint = avg +#endif + ! What to print? + ! ncol, nlay, ngpt; are clouds used, are aerosols used; time per column, total, min; + print *, " ncol nlay ngpt clouds aerosols time_per_col_ms nloops time_total_s time_min_s" + write(output_unit, '(3(i6, 1x), 6x, 2(i1, 8x), 1x, f7.3, 1x, i6, 2x, 2(4x,f7.3))') & + ncol, nlay, ngpt, merge(1,0,do_clouds), merge(1,0,do_aerosols), & + avg/(real(ncol*nloops) * (1.0e-3*clock_rate)), nloops, avg / real(clock_rate), mint / real(clock_rate) if(.true.) call write_fluxes From 2c1bae84a8e28648e996eed61673e3fcf757b270 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 26 Jun 2023 22:31:43 -0400 Subject: [PATCH 27/38] Conditional evaluation when computing profiles --- examples/all-sky/rrtmgp_allsky.F90 | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index c107de774..b51008982 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -423,7 +423,8 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, real(wp), parameter :: gamma = 6.7e-3 real(wp), parameter :: q_0 = 0.01864 ! for 300 K SST. - + ! ------------------- + Tv0 = (1. + 0.608*q_0) * SST ! ! Split resolution above and below RCE tropopause (15 km or about 125 hPa) ! @@ -443,16 +444,16 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, do ilay = 1, nlay do icol = 1, ncol z = z_lay(ilay) - q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) - T = SST - gamma*z / (1. + 0.608*q) - Tv = (1. + 0.608*q ) * T - Tv0 = (1. + 0.608*q_0) * SST - p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) + else + q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) + T = SST - gamma*z / (1. + 0.608*q) + Tv = (1. + 0.608*q ) * T + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) end if p_lay(icol,ilay) = p t_lay(icol,ilay) = T @@ -468,16 +469,16 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, do ilay = 1, nlay+1 do icol = 1, ncol z = z_lev(ilay) - q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) - T = SST - gamma*z / (1. + 0.608*q) - Tv = (1. + 0.608*q ) * T - Tv0 = (1. + 0.608*q_0) * SST - p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) + else + q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) + T = SST - gamma*z / (1. + 0.608*q) + Tv = (1. + 0.608*q ) * T + p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) end if p_lev(icol,ilay) = p t_lev(icol,ilay) = T From 6446ab4f8eb2e25b57d18a6a712a21c9c3972875 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 27 Jun 2023 15:53:28 -0400 Subject: [PATCH 28/38] Add Python script to loop over ncol, nlay while specifying clouds, aerosols --- examples/all-sky/make_problem_size_loop.py | 63 ++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 examples/all-sky/make_problem_size_loop.py diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py new file mode 100644 index 000000000..8aa1cf4ce --- /dev/null +++ b/examples/all-sky/make_problem_size_loop.py @@ -0,0 +1,63 @@ +#! /usr/bin/env python +# +# This script... +# +import argparse +import csv +import os + +# template: exe ncol nlay nloops output_file k_dist clouds aeorsols + +# specify: kdist, optional clouds, aerosols. specify nloops +# Toggle clouds and aerosols? +# Loop over sets of ncol, nlay, +# output name + + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description="Description here ") + # Argument parseing described at + # https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse + parser.add_argument("-x", "--executable", + type=str, + default="./rrtmgp_allsky", + help="Path to exectuable") + parser.add_argument("-c", "--ncol", + type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + default="2,4,8,16", + help="Number of columns (multiple e.g. 2,4,8,16)") + parser.add_argument("-l", "--nlay", + type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + default="32, 64, 96", + help="Number of layers (multiple e.g. 32,64.96)") + parser.add_argument("-i", "--nloops", + type=int, default=1, + help="Number of loops (same for all ncol)") + parser.add_argument("-o", "--output_file", + type=str, + default="rrtmgp-allsky-output.nc", + help="Path to output file") + parser.add_argument("-k", "--k_distribution", + type=str, + required = True, + help="Path to gas optics file [required]") + parser.add_argument("-cl", "--cloud_optics", + type=str, default="", + help="Path to cloud optics file") + parser.add_argument("-a", "--aerosol_optics", + type=str, default="", + help="Path to aerosol optics file") + args = parser.parse_args() + + # Can't supply aerosols without clouds + if(args.cloud_optics == "" and args.aerosol_optics != ""): + raise AssertionError("Need to supply cloud optics if providing aerosol optics") + + print (args.output_file) + # Every combo of ncol, nlay + for l in args.nlay: + for i in args.ncol: + print(f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + \ + f"{args.output_file} {args.k_distribution}" + \ + f"{args.cloud_optics} {args.aerosol_optics}") From 2f512a2cc016f1f3c839c0eac26f6ece73ca4fbe Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 27 Jun 2023 16:03:05 -0400 Subject: [PATCH 29/38] Extra print --- examples/all-sky/make_problem_size_loop.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py index 8aa1cf4ce..1f45d02da 100644 --- a/examples/all-sky/make_problem_size_loop.py +++ b/examples/all-sky/make_problem_size_loop.py @@ -54,7 +54,6 @@ if(args.cloud_optics == "" and args.aerosol_optics != ""): raise AssertionError("Need to supply cloud optics if providing aerosol optics") - print (args.output_file) # Every combo of ncol, nlay for l in args.nlay: for i in args.ncol: From 0d96bdb07d07f2bc32e2812c42530b471babad1b Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 27 Jun 2023 16:45:34 -0400 Subject: [PATCH 30/38] Explicit memory release in all-sky example --- examples/all-sky/rrtmgp_allsky.F90 | 44 ++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 6 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index b51008982..a3a97f460 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -330,8 +330,8 @@ program rte_rrtmgp_allsky !$omp end target data else - !$acc enter data create( toa_flux) - !$omp target enter data map(alloc:toa_flux) + !$acc data create( toa_flux) + !$omp target data map(alloc:toa_flux) fluxes%flux_dn_dir => flux_dir(:,:) call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & @@ -351,10 +351,9 @@ program rte_rrtmgp_allsky mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - !$acc exit data delete( toa_flux) - !$omp target exit data map(release:toa_flux) + !$acc end data + !$omp end target data end if - !print *, "******************************************************************" call system_clock(finish, clock_rate) elapsed(iloop) = finish - start end do @@ -377,6 +376,10 @@ program rte_rrtmgp_allsky if(.true.) call write_fluxes + call stop_on_err(fluxes%finalize()) + + ! Memory for bounday conditions on the GPU was allocated with unstructured data dataments + ! (acc enter data). Deallocate it expliicity if(is_lw) then !$acc exit data delete( t_sfc, emis_sfc) !$omp target exit data map(release:t_sfc, emis_sfc) @@ -385,7 +388,36 @@ program rte_rrtmgp_allsky !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if - ! Clouds and aerosols used enter data - we won't release explicitly + ! Clouds and aerosols also used enter data + if(do_clouds) then + !$acc exit data delete( cloud_mask, lwp, iwp, rel, rei) + !$omp target exit data map(release:cloud_mask, lwp, iwp, rel, rei) + select type(clouds) + class is (ty_optical_props_1scl) + !$acc exit data delete (clouds%tau, clouds) + !$omp target exit data map(release:clouds%tau) + class is (ty_optical_props_2str) + !$acc exit data delete (clouds%tau, clouds%ssa, clouds%g, clouds) + !$omp target exit data map(release:clouds%tau, clouds%ssa, clouds%g) + end select + end if + if(do_aerosols) then + !$acc exit data delete( aero_type, aero_size, aero_mass, relhum) + !$omp target exit data map(release:aero_type, aero_size, aero_mass, relhum) + select type(aerosols) + class is (ty_optical_props_1scl) + !$acc exit data delete (aerosols%tau, aerosols) + !$omp target exit data map(release:aerosols%tau) + class is (ty_optical_props_2str) + !$acc exit data delete (aerosols%tau, aerosols%ssa, aerosols%g, aerosols) + !$omp target exit data map(release:aerosols%tau, aerosols%ssa, aerosols%g) + end select + end if + + if(.not. is_lw) then + !$acc exit data delete( flux_dir) + !$omp target exit data map(release:flux_dir) + end if ! fluxes - but not flux_dir, which used enter data !$acc end data From 0d482e0799505085231a76afa73aa2fb975e284e Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 27 Jun 2023 16:47:29 -0400 Subject: [PATCH 31/38] Misplaced finalize() --- examples/all-sky/rrtmgp_allsky.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index a3a97f460..a10d1cd27 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -376,8 +376,6 @@ program rte_rrtmgp_allsky if(.true.) call write_fluxes - call stop_on_err(fluxes%finalize()) - ! Memory for bounday conditions on the GPU was allocated with unstructured data dataments ! (acc enter data). Deallocate it expliicity if(is_lw) then From 1b0bb64b876cd7bd97316f64271173695549ae25 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 27 Jun 2023 17:16:05 -0400 Subject: [PATCH 32/38] Revised timing info in all-sky example --- examples/all-sky/rrtmgp_allsky.F90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index a10d1cd27..ef2617520 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -360,19 +360,15 @@ program rte_rrtmgp_allsky ! call system_clock(finish_all, clock_rate) -#if defined(_OPENACC) || defined(_OPENMP) avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) mint = minval(elapsed) -#else - avg = (finish_all-start_all) - mint = avg -#endif + ! What to print? ! ncol, nlay, ngpt; are clouds used, are aerosols used; time per column, total, min; print *, " ncol nlay ngpt clouds aerosols time_per_col_ms nloops time_total_s time_min_s" write(output_unit, '(3(i6, 1x), 6x, 2(i1, 8x), 1x, f7.3, 1x, i6, 2x, 2(4x,f7.3))') & ncol, nlay, ngpt, merge(1,0,do_clouds), merge(1,0,do_aerosols), & - avg/(real(ncol*nloops) * (1.0e-3*clock_rate)), nloops, avg / real(clock_rate), mint / real(clock_rate) + avg/(real(ncol) * (1.0e-3*clock_rate)), nloops, sum(elapsed) / real(clock_rate), mint / real(clock_rate) if(.true.) call write_fluxes From 039f6c84eb0a8f5fc2533ce2f50f9ae1a547e7a2 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 28 Jun 2023 13:21:53 -0400 Subject: [PATCH 33/38] Explicit memory release; no checks in any loop iteration beyond the first --- examples/all-sky/rrtmgp_allsky.F90 | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index ef2617520..164b3608b 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -287,7 +287,7 @@ program rte_rrtmgp_allsky !!$omp parallel do firstprivate(fluxes) do iloop = 1, nloops ! Omit the checks starting with the second iteration - if (iloop == 2) call rte_config_checks(logical(.false., wl)) + if (iloop > 1) call rte_config_checks(logical(.false., wl)) call system_clock(start) ! @@ -372,8 +372,10 @@ program rte_rrtmgp_allsky if(.true.) call write_fluxes + ! ! Memory for bounday conditions on the GPU was allocated with unstructured data dataments ! (acc enter data). Deallocate it expliicity + ! if(is_lw) then !$acc exit data delete( t_sfc, emis_sfc) !$omp target exit data map(release:t_sfc, emis_sfc) @@ -381,8 +383,10 @@ program rte_rrtmgp_allsky !$acc exit data delete( sfc_alb_dir, sfc_alb_dif, mu0) !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if - + + ! ! Clouds and aerosols also used enter data + ! if(do_clouds) then !$acc exit data delete( cloud_mask, lwp, iwp, rel, rei) !$omp target exit data map(release:cloud_mask, lwp, iwp, rel, rei) @@ -394,6 +398,11 @@ program rte_rrtmgp_allsky !$acc exit data delete (clouds%tau, clouds%ssa, clouds%g, clouds) !$omp target exit data map(release:clouds%tau, clouds%ssa, clouds%g) end select + ! + ! Explicit finalization of cloud optical properties - not really necessary since memory + ! will be freed when the program ends, but useful for testing + ! + call clouds%finalize end if if(do_aerosols) then !$acc exit data delete( aero_type, aero_size, aero_mass, relhum) @@ -406,8 +415,17 @@ program rte_rrtmgp_allsky !$acc exit data delete (aerosols%tau, aerosols%ssa, aerosols%g, aerosols) !$omp target exit data map(release:aerosols%tau, aerosols%ssa, aerosols%g) end select + ! + ! Explicit finalization of aerosol optical properties - not really necessary since memory + ! will be freed when the program ends, but useful for testing + ! + call aerosols%finalize end if - + ! + ! k-distribution + ! + call k_dist%finalize + if(.not. is_lw) then !$acc exit data delete( flux_dir) !$omp target exit data map(release:flux_dir) From 905a1fe54f0dec77d097e7961f3f04417f272adb Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 28 Jun 2023 13:32:28 -0400 Subject: [PATCH 34/38] Update run-allsky-example.py --- examples/all-sky/run-allsky-example.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/examples/all-sky/run-allsky-example.py b/examples/all-sky/run-allsky-example.py index 3d9dcac03..2b7f604a4 100644 --- a/examples/all-sky/run-allsky-example.py +++ b/examples/all-sky/run-allsky-example.py @@ -21,8 +21,6 @@ # In the local directory all_sky_exe_name = os.path.join(all_sky_dir, "rrtmgp_allsky") -input_file = os.path.join(os.environ["RRTMGP_DATA"], "examples", "all-sky", "inputs", "garand-atmos-1.nc") -atmos_file = os.path.join(all_sky_dir, "rrtmgp-allsky.nc") if __name__ == '__main__': parser = argparse.ArgumentParser( @@ -30,28 +28,27 @@ parser.add_argument("--run_command", type=str, default="", help="Prefix ('jsrun' etc.) for running commands. " "Use quote marks to enclose multi-part commands.") - parser.add_argument("--ncol", type=int, default=128, - help="Number of cloudy columns to compute " + parser.add_argument("--ncol", type=int, default=24, + help="Number of columns to compute") + parser.add_argument("--nlay", type=int, default=72, + help="Number of layers " "(every one will have the same clouds)") parser.add_argument("--nloops", type=int, default=1, help="Number of times to compute 'nloops' " "cloudy columns") args = parser.parse_args() - ncol_str = '{0:5d}'.format(args.ncol) + ncol_str = '{0:5d}'.format(args.ncol) + nlay_str = '{0:5d}'.format(args.nlay) nloops_str = '{0:5d}'.format(args.nloops) if args.run_command: print("using the run command") all_sky_exe_name = args.run_command + " " + all_sky_exe_name - os.chdir(all_sky_dir) # Remove cloudy-sky fluxes from the file containing the atmospheric profiles - shutil.copyfile(input_file, atmos_file) subprocess.run( - [all_sky_exe_name, atmos_file, lw_gas_coeffs_file, lw_clouds_coeff_file, - ncol_str, nloops_str]) + [all_sky_exe_name, ncol_str, nlay_str, nloops_str, "rrtmgp-allsky-lw-no-aerosols.nc", lw_gas_coeffs_file, lw_clouds_coeff_file]) subprocess.run( - [all_sky_exe_name, atmos_file, sw_gas_coeffs_file, sw_clouds_coeff_file, - ncol_str, nloops_str]) + [all_sky_exe_name, ncol_str, nlay_str, nloops_str, "rrtmgp-allsky-sw-no-aerosols.nc", sw_gas_coeffs_file, sw_clouds_coeff_file]) # end main From adddafa7b7e7a7432dbb3ee5883b0a211af2dde7 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 3 Jul 2023 16:07:04 -0400 Subject: [PATCH 35/38] OpenACC/OpenMP memory management --- .../mo_aerosol_optics_rrtmgp_merra.F90 | 54 ++++++++++--------- rte-frontend/mo_rte_util_array_validation.F90 | 24 ++++----- 2 files changed, 42 insertions(+), 36 deletions(-) diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index b2d4ce997..db1249c23 100644 --- a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 +++ b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 @@ -54,7 +54,7 @@ module mo_aerosol_optics_rrtmgp_merra ! index identifiers for aerosol optical property tables integer, parameter, private :: ext = 1 ! extinction integer, parameter, private :: ssa = 2 ! single scattering albedo - integer, parameter, private :: g = 3 ! asymmetry parameter + integer, parameter, private :: g = 3 ! asymmetry parameter private ! ----------------------------------------------------------------------------------- @@ -66,7 +66,7 @@ module mo_aerosol_optics_rrtmgp_merra ! Table upper and lower aerosol size (radius) bin limits (microns) real(wp),dimension(:,:), allocatable :: merra_aero_bin_lims ! Dimensions (pair,nbin) ! Table relative humidity values - real(wp),dimension(:), allocatable :: aero_rh(:) + real(wp),dimension(:), allocatable :: aero_rh(:) ! ! The aerosol tables themselves. ! extinction (m2/kg) @@ -83,10 +83,9 @@ module mo_aerosol_optics_rrtmgp_merra ! ! ----- contains - generic, public :: load => load_lut - procedure, public :: finalize - procedure, public :: aerosol_optics - + generic, public :: load => load_lut + procedure, public :: finalize + procedure, public :: aerosol_optics ! Internal procedures procedure, private :: load_lut end type ty_aerosol_optics_rrtmgp_merra @@ -194,10 +193,6 @@ function load_lut(this, band_lims_wvn, & this%aero_ocar_rh_tbl = aero_ocar_rh_tbl !$acc end kernels !$omp end target - ! - - - end function load_lut !-------------------------------------------------------------------------------------------------------------------- ! @@ -209,16 +204,13 @@ subroutine finalize(this) ! Lookup table aerosol optics interpolation arrays if(allocated(this%merra_aero_bin_lims)) then - deallocate(this%merra_aero_bin_lims, this%aero_rh) + !$acc exit data delete( this%merra_aero_bin_lims, this%aero_rh) + !$omp target exit data map(release:this%merra_aero_bin_lims, this%aero_rh) & end if ! Lookup table aerosol optics coefficients if(allocated(this%aero_dust_tbl)) then - - deallocate(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl, & - this%aero_bcar_tbl, this%aero_bcar_rh_tbl, & - this%aero_ocar_tbl, this%aero_ocar_rh_tbl) !$acc exit data delete(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & !$acc delete(this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & !$acc delete(this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & @@ -226,6 +218,9 @@ subroutine finalize(this) !$omp target exit data map(release:this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & !$omp map(release:this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & !$omp map(release:this%aero_ocar_tbl, this%aero_ocar_rh_tbl) + deallocate(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl, & + this%aero_bcar_tbl, this%aero_bcar_rh_tbl, & + this%aero_ocar_tbl, this%aero_ocar_rh_tbl) end if end subroutine finalize @@ -271,6 +266,8 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & integer :: icol, ilay, ibnd, ibin ! scalars for total tau, tau*ssa real(wp) :: tau, taussa + ! Scalars to work around OpenACC/OMP issues + real(wp) :: minSize, maxSize ! ---------------------------------------- ! @@ -290,6 +287,8 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & nrh = size(this%aero_rh,1) nval = size(this%aero_dust_tbl,1) nbnd = size(this%aero_dust_tbl,3) + minSize = this%merra_aero_bin_lims(1,1) + maxSize = this%merra_aero_bin_lims(2,nbin) ! ! Array sizes @@ -322,14 +321,14 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & if(error_msg /= "") return end if - !$acc data copyin(aero_type, aero_size, aero_mass, relhum) & - !$acc create(atau, ataussa, ataussag, aeromsk) - !$omp target data map(to:aero_type, aero_size, aero_mass, relhum) & - !$omp map(alloc:atau, ataussa, ataussag, aeromsk) + !$acc data copyin(aero_type, aero_size, aero_mass, relhum) + !$omp target data map(to:aero_type, aero_size, aero_mass, relhum) ! ! Aerosol mask; don't need aerosol optics if there's no aerosol ! - !$acc parallel loop gang vector default(none) collapse(2) + !$acc data create(aeromsk) + !$omp target data map(alloc:aeromsk) + !$acc parallel loop default(none) collapse(2) !$omp target teams distribute parallel do simd collapse(2) do ilay = 1, nlay do icol = 1, ncol @@ -341,13 +340,18 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & ! Aerosol size, relative humidity ! if(check_values) then - if(any_vals_outside(aero_size, aeromsk, & - this%merra_aero_bin_lims(1,1), this%merra_aero_bin_lims(2,nbin))) & + if(any_vals_outside(aero_size, aeromsk, minSize, maxSize) & error_msg = 'aerosol optics: requested aerosol size is out of bounds' - if(any_vals_outside(relhum, 0._wp, 1._wp)) & + if(any_vals_outside(relhum, aeromsk, 0._wp, 1._wp)) & error_msg = 'aerosol optics: relative humidity fraction is out of bounds' end if + ! Release aerosol mask + !$acc end data + !$omp end target data + if(error_msg == "") then + !$acc data create(atau, ataussa, ataussag) + !$omp target data map(alloc:atau, ataussa, ataussag) ! ! ! ---------------------------------------- @@ -411,9 +415,11 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & type is (ty_optical_props_nstr) error_msg = "aerosol optics: n-stream calculations not yet supported" end select + !$acc end data + !$omp end target data end if !$acc end data - !$omp end target data + !$omp end target data end function aerosol_optics !-------------------------------------------------------------------------------------------------------------------- ! diff --git a/rte-frontend/mo_rte_util_array_validation.F90 b/rte-frontend/mo_rte_util_array_validation.F90 index 6acd99b7e..79584efe3 100644 --- a/rte-frontend/mo_rte_util_array_validation.F90 +++ b/rte-frontend/mo_rte_util_array_validation.F90 @@ -123,8 +123,8 @@ logical function any_vals_less_than_1D_masked(array, mask, check_value) real(wp) :: minValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue) minValue = minval(array, mask=mask) !$acc end kernels !$omp end target @@ -140,8 +140,8 @@ logical function any_vals_less_than_2D_masked(array, mask, check_value) real(wp) :: minValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue) minValue = minval(array, mask=mask) !$acc end kernels !$omp end target @@ -157,8 +157,8 @@ logical function any_vals_less_than_3D_masked(array, mask, check_value) real(wp) :: minValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue) minValue = minval(array, mask=mask) !$acc end kernels !$omp end target @@ -249,8 +249,8 @@ logical function any_vals_outside_1D_masked(array, mask, checkMin, checkMax) real(wp) :: minValue, maxValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue, maxValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue, maxValue) minValue = minval(array, mask=mask) maxValue = maxval(array, mask=mask) !$acc end kernels @@ -266,8 +266,8 @@ logical function any_vals_outside_2D_masked(array, mask, checkMin, checkMax) real(wp) :: minValue, maxValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue, maxValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue, maxValue) minValue = minval(array, mask=mask) maxValue = maxval(array, mask=mask) !$acc end kernels @@ -283,8 +283,8 @@ logical function any_vals_outside_3D_masked(array, mask, checkMin, checkMax) real(wp) :: minValue, maxValue - !$acc kernels copyin(array) - !$omp target map(to:array, mask) map(from:minValue, maxValue) + !$acc kernels copyin(array, mask) + !$omp target map(to: array, mask) map(from:minValue, maxValue) minValue = minval(array, mask=mask) maxValue = maxval(array, mask=mask) !$acc end kernels From 453857b5d48bfffab91b1780782e873b2837a387 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 3 Jul 2023 16:09:08 -0400 Subject: [PATCH 36/38] Syntax.... --- rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index db1249c23..ce66114e8 100644 --- a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 +++ b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 @@ -340,7 +340,7 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & ! Aerosol size, relative humidity ! if(check_values) then - if(any_vals_outside(aero_size, aeromsk, minSize, maxSize) & + if(any_vals_outside(aero_size, aeromsk, minSize, maxSize)) & error_msg = 'aerosol optics: requested aerosol size is out of bounds' if(any_vals_outside(relhum, aeromsk, 0._wp, 1._wp)) & error_msg = 'aerosol optics: relative humidity fraction is out of bounds' From fb3e92a9205827014fa61ef8df5eea5b716b8799 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 3 Jul 2023 16:18:41 -0400 Subject: [PATCH 37/38] More robust updating of derived types in aerosol optics --- rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index ce66114e8..60ef85ee4 100644 --- a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 +++ b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 @@ -287,6 +287,9 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & nrh = size(this%aero_rh,1) nval = size(this%aero_dust_tbl,1) nbnd = size(this%aero_dust_tbl,3) + + !$acc update host(this%merra_aero_bin_lims) + !$omp target update from(this%merra_aero_bin_lims) minSize = this%merra_aero_bin_lims(1,1) maxSize = this%merra_aero_bin_lims(2,nbin) From 95c051a31c088a1645690c4068df5df8403d4336 Mon Sep 17 00:00:00 2001 From: Dmitry Alexeev Date: Tue, 4 Jul 2023 08:55:39 -0700 Subject: [PATCH 38/38] reorder aerosol LUTs for much better GPU perf --- .../mo_aerosol_optics_rrtmgp_merra.F90 | 80 +++++++++---------- 1 file changed, 39 insertions(+), 41 deletions(-) diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index 60ef85ee4..ae020f94a 100644 --- a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 +++ b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 @@ -161,38 +161,36 @@ function load_lut(this, band_lims_wvn, & allocate(this%aero_rh(nrh)) ! Allocate LUT coefficients allocate(this%aero_dust_tbl(nval, nbin, nband), & - this%aero_salt_tbl(nval, nrh, nbin, nband), & - this%aero_sulf_tbl(nval, nrh, nband), & + this%aero_salt_tbl(nrh, nval, nbin, nband), & + this%aero_sulf_tbl(nrh, nval, nband), & this%aero_bcar_tbl(nval, nband), & - this%aero_bcar_rh_tbl(nval, nrh, nband), & + this%aero_bcar_rh_tbl(nrh, nval, nband), & this%aero_ocar_tbl(nval, nband), & - this%aero_ocar_rh_tbl(nval, nrh, nband)) - !$acc enter data create(this) & - !$acc create(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & - !$acc create(this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & - !$acc create(this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & - !$acc create(this%merra_aero_bin_lims, this%aero_rh) - !$omp target enter data & - !$omp map(alloc:this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & - !$omp map(alloc:this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & - !$omp map(alloc:this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & - !$omp map(alloc:this%merra_aero_bin_lims, this%aero_rh) - - - ! Load LUT coefficients - !$acc kernels - !$omp target + this%aero_ocar_rh_tbl(nrh, nval, nband)) + + ! Copy LUT coefficients this%merra_aero_bin_lims = merra_aero_bin_lims this%aero_rh = aero_rh this%aero_dust_tbl = aero_dust_tbl - this%aero_salt_tbl = aero_salt_tbl - this%aero_sulf_tbl = aero_sulf_tbl this%aero_bcar_tbl = aero_bcar_tbl - this%aero_bcar_rh_tbl = aero_bcar_rh_tbl this%aero_ocar_tbl = aero_ocar_tbl - this%aero_ocar_rh_tbl = aero_ocar_rh_tbl - !$acc end kernels - !$omp end target + + this%aero_salt_tbl = reshape( aero_salt_tbl, shape=(/nrh, nval, nbin, nband/), order=(/2,1,3,4/) ) + this%aero_sulf_tbl = reshape( aero_sulf_tbl, shape=(/nrh, nval, nband/), order=(/2,1,3/) ) + this%aero_bcar_rh_tbl = reshape( aero_bcar_rh_tbl, shape=(/nrh, nval, nband/), order=(/2,1,3/) ) + this%aero_ocar_rh_tbl = reshape( aero_ocar_rh_tbl, shape=(/nrh, nval, nband/), order=(/2,1,3/) ) + + !$acc enter data create(this) & + !$acc copyin(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & + !$acc copyin(this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & + !$acc copyin(this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & + !$acc copyin(this%merra_aero_bin_lims, this%aero_rh) + !$omp target enter data & + !$omp map(to:this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & + !$omp map(to:this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & + !$omp map(to:this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & + !$omp map(to:this%merra_aero_bin_lims, this%aero_rh) + end function load_lut !-------------------------------------------------------------------------------------------------------------------- ! @@ -451,11 +449,11 @@ subroutine compute_all_from_table(ncol, nlay, npair, nval, nrh, nbin, nbnd, real(wp), dimension(nrh), intent(in) :: aero_rh real(wp), dimension(nval, nbin,nbnd), intent(in) :: aero_dust_tbl - real(wp), dimension(nval,nrh,nbin,nbnd), intent(in) :: aero_salt_tbl - real(wp), dimension(nval,nrh, nbnd), intent(in) :: aero_sulf_tbl - real(wp), dimension(nval,nrh, nbnd), intent(in) :: aero_bcar_rh_tbl + real(wp), dimension(nrh,nval,nbin,nbnd), intent(in) :: aero_salt_tbl + real(wp), dimension(nrh,nval, nbnd), intent(in) :: aero_sulf_tbl + real(wp), dimension(nrh,nval, nbnd), intent(in) :: aero_bcar_rh_tbl real(wp), dimension(nval, nbnd), intent(in) :: aero_bcar_tbl - real(wp), dimension(nval,nrh, nbnd), intent(in) :: aero_ocar_rh_tbl + real(wp), dimension(nrh,nval, nbnd), intent(in) :: aero_ocar_rh_tbl real(wp), dimension(nval, nbnd), intent(in) :: aero_ocar_tbl real(wp), dimension(ncol,nlay,nbnd), intent(out) :: tau, taussa, taussag @@ -507,28 +505,28 @@ subroutine compute_all_from_table(ncol, nlay, npair, nval, nrh, nbin, nbnd, ! sea-salt case(merra_aero_salt) tau (icol,ilay,ibnd) = mass (icol,ilay) * & - linear_interp_aero_table(aero_salt_tbl(ext,:,ibin,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_salt_tbl(:,ext,ibin,ibnd),irh1,irh2,rdrh) taussa (icol,ilay,ibnd) = tau (icol,ilay,ibnd) * & - linear_interp_aero_table(aero_salt_tbl(ssa,:,ibin,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_salt_tbl(:,ssa,ibin,ibnd),irh1,irh2,rdrh) taussag(icol,ilay,ibnd) = taussa(icol,ilay,ibnd) * & - linear_interp_aero_table(aero_salt_tbl(g, :,ibin,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_salt_tbl(:,g, ibin,ibnd),irh1,irh2,rdrh) ! sulfate case(merra_aero_sulf) tau (icol,ilay,ibnd) = mass (icol,ilay) * & - linear_interp_aero_table(aero_sulf_tbl(ext,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_sulf_tbl(:,ext,ibnd),irh1,irh2,rdrh) taussa (icol,ilay,ibnd) = tau (icol,ilay,ibnd) * & - linear_interp_aero_table(aero_sulf_tbl(ssa,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_sulf_tbl(:,ssa,ibnd),irh1,irh2,rdrh) taussag(icol,ilay,ibnd) = taussa(icol,ilay,ibnd) * & - linear_interp_aero_table(aero_sulf_tbl(g, :,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_sulf_tbl(:,g, ibnd),irh1,irh2,rdrh) ! black carbon - hydrophilic case(merra_aero_bcar_rh) tau (icol,ilay,ibnd) = mass (icol,ilay) * & - linear_interp_aero_table(aero_bcar_rh_tbl(ext,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_bcar_rh_tbl(:,ext,ibnd),irh1,irh2,rdrh) taussa (icol,ilay,ibnd) = tau (icol,ilay,ibnd) * & - linear_interp_aero_table(aero_bcar_rh_tbl(ssa,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_bcar_rh_tbl(:,ssa,ibnd),irh1,irh2,rdrh) taussag(icol,ilay,ibnd) = taussa(icol,ilay,ibnd) * & - linear_interp_aero_table(aero_bcar_rh_tbl(g, :,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_bcar_rh_tbl(:,g, ibnd),irh1,irh2,rdrh) ! black carbon - hydrophobic case(merra_aero_bcar) tau (icol,ilay,ibnd) = mass (icol,ilay) * aero_bcar_tbl(ext,ibnd) @@ -537,11 +535,11 @@ subroutine compute_all_from_table(ncol, nlay, npair, nval, nrh, nbin, nbnd, ! organic carbon - hydrophilic case(merra_aero_ocar_rh) tau (icol,ilay,ibnd) = mass (icol,ilay) * & - linear_interp_aero_table(aero_ocar_rh_tbl(ext,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_ocar_rh_tbl(:,ext,ibnd),irh1,irh2,rdrh) taussa (icol,ilay,ibnd) = tau (icol,ilay,ibnd) * & - linear_interp_aero_table(aero_ocar_rh_tbl(ssa,:,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_ocar_rh_tbl(:,ssa,ibnd),irh1,irh2,rdrh) taussag(icol,ilay,ibnd) = taussa(icol,ilay,ibnd) * & - linear_interp_aero_table(aero_ocar_rh_tbl(g, :,ibnd),irh1,irh2,rdrh) + linear_interp_aero_table(aero_ocar_rh_tbl(:,g, ibnd),irh1,irh2,rdrh) ! organic carbon - hydrophobic case(merra_aero_ocar) tau (icol,ilay,ibnd) = mass (icol,ilay) * aero_ocar_tbl(ext,ibnd)