From 36d79691feddf32b63b569f157494ae7847c548a Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 4 Jul 2023 16:42:56 -0400 Subject: [PATCH] All-sky example useful for timing (#220) Revises the all-sky example to be useful for timing benchmarks. Removes use of NCHOME environment variable, removes (optional) use of GPTL timing library in RFMIP examples. Revises aerosol optics for GPUs. --- .github/workflows/containerized-ci.yml | 2 +- .github/workflows/continuous-integration.yml | 1 + .github/workflows/self-hosted-ci.yml | 3 +- examples/all-sky/Makefile | 29 +- examples/all-sky/compare-to-reference.py | 10 +- examples/all-sky/make_problem_size_loop.py | 62 ++ examples/all-sky/rrtmgp_allsky.F90 | 927 +++++++++++------- examples/all-sky/run-allsky-example.py | 19 +- examples/{all-sky => }/mo_garand_atmos_io.F90 | 4 + examples/rfmip-clear-sky/Makefile | 16 - examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 | 41 - examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 | 38 - gas-optics/mo_gas_concentrations.F90 | 1 - .../mo_aerosol_optics_rrtmgp_merra.F90 | 131 +-- rte-frontend/mo_rte_util_array_validation.F90 | 24 +- 15 files changed, 772 insertions(+), 536 deletions(-) create mode 100644 examples/all-sky/make_problem_size_loop.py rename examples/{all-sky => }/mo_garand_atmos_io.F90 (96%) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index b9e70ccdd..8b9fdc028 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 @@ -76,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..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 }} @@ -62,6 +60,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 21e129a0d..4cdb3207b 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -28,8 +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_garand_atmos_io.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 # @@ -41,20 +41,27 @@ 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_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_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: - 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 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 + 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 + -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/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py new file mode 100644 index 000000000..1f45d02da --- /dev/null +++ b/examples/all-sky/make_problem_size_loop.py @@ -0,0 +1,62 @@ +#! /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") + + # 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}") diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index 65b983003..164b3608b 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -1,84 +1,6 @@ -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) - 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 -! -------------------------------------------------------------------------------------- -! -! 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 +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 @@ -95,21 +17,19 @@ 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 ! ---------------------------------------------------------------------------------- ! 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 ! @@ -118,6 +38,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 @@ -133,27 +58,16 @@ 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 ! Relative humidity (fraction) logical, dimension(:,:), allocatable :: aero_mask ! Aerosol mask - real(wp), allocatable, dimension(:,:) :: vmr_h2o - ! h2o vmr ! ! Output variables @@ -177,23 +91,22 @@ 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 - logical :: use_luts = .true., write_fluxes = .true. + integer :: nUserArgs, nloops, ncol, nlay + ! logical :: write_fluxes = .false. + logical :: do_clouds = .false., use_luts = .true. + logical :: do_aerosols = .false. 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 + character(len=256) :: output_file, k_dist_file, cloud_optics_file, aerosol_optics_file ! ! 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 ) @@ -201,88 +114,76 @@ program rte_rrtmgp_clouds_aerosols ! Code ! ---------------------------------------------------------------------------------- ! - ! Parse command line for any file names, block size - ! - 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 + ! Parse command line: rrtmgp_allsky ncol nlay nreps kdist [clouds [aerosols]] ! - ! 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) + nUserArgs = command_argument_count() + 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 + 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,output_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..." + ! ----------------------------------------------------------------------------------- + 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) - ! 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) + 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) 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) + end if ! ---------------------------------------------------------------------------- ! @@ -298,17 +199,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())) - ! Clouds 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. ! @@ -323,32 +216,9 @@ 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 ") - 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 ") + call stop_on_err("rte_rrtmgp_allsky: Don't recognize the kind of optical properties ") end select + ! ---------------------------------------------------------------------------- ! Boundary conditions depending on whether the k-distribution being supplied ! is LW or SW @@ -359,8 +229,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 @@ -376,8 +246,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 @@ -394,91 +264,16 @@ 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) !$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 ! ---------------------------------------------------------------------------- ! @@ -492,26 +287,29 @@ program rte_rrtmgp_clouds_aerosols !!$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) ! ! 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 ! 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) & @@ -522,8 +320,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, & @@ -532,8 +330,8 @@ program rte_rrtmgp_clouds_aerosols !$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, & @@ -541,53 +339,510 @@ 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, & 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 ! 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) - !$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)) + avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) + mint = minval(elapsed) + + ! 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) * (1.0e-3*clock_rate)), nloops, sum(elapsed) / real(clock_rate), mint / real(clock_rate) - 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) -#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 + 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 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 - !$acc enter data create(lwp, iwp, rel, rei) - !$omp target enter data map(alloc:lwp, iwp, rel, rei) -end program rte_rrtmgp_clouds_aerosols + + ! + ! 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 + ! + ! 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) + !$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 + ! + ! 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) + end if + + ! fluxes - but not flux_dir, which used enter data + !$acc end data + !$omp end target data + ! p_lay etc + !$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) + ! + ! 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_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) :: z, q, T, p + real(wp) :: Tv, Tv0, 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 + + 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. + ! ------------------- + Tv0 = (1. + 0.608*q_0) * SST + ! + ! 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)) + + !$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 + ! 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 + do icol = 1, ncol + z = z_lay(ilay) + 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 + 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) + end do + end do + + !$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) + 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 + end do + end do + !$acc end data + !$omp end target data + end subroutine compute_profiles + ! ---------------------------------------------------------------------------------- + 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 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_allsky: 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 + 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 + 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_allsky: Don't recognize the kind of optical properties ") + end select + ! + ! Derive relative humidity from profile + ! 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)) + ! + ! Aerosol properties + ! + allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & + 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 + ! 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 + 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 + 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 + 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 + !$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 + 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 + !-------------------------------------------------------------------------------------------------------------------- + subroutine write_fluxes + use netcdf + use mo_simple_netcdf, only: write_field + 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(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") + 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 + ! + 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 + ! + 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 + 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??") + + ! + ! Write variables + ! + ! State - writing + !$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)) + 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 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)) + call stop_on_err(write_field(ncid, "rei", rei)) + end if + + if(do_aerosols) then + !$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 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 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 + ! --------------------------------------------------------- + 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_allsky 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 diff --git a/examples/all-sky/mo_garand_atmos_io.F90 b/examples/mo_garand_atmos_io.F90 similarity index 96% rename from examples/all-sky/mo_garand_atmos_io.F90 rename to examples/mo_garand_atmos_io.F90 index 5b730c9ab..dde3bfbf0 100644 --- a/examples/all-sky/mo_garand_atmos_io.F90 +++ b/examples/mo_garand_atmos_io.F90 @@ -70,11 +70,15 @@ 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) 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)) 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 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) diff --git a/gas-optics/mo_gas_concentrations.F90 b/gas-optics/mo_gas_concentrations.F90 index 26de09998..592baf57d 100644 --- a/gas-optics/mo_gas_concentrations.F90 +++ b/gas-optics/mo_gas_concentrations.F90 @@ -407,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 diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index b2d4ce997..ae020f94a 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 @@ -162,41 +161,35 @@ 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 !-------------------------------------------------------------------------------------------------------------------- @@ -209,16 +202,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 +216,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 +264,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 ! ---------------------------------------- ! @@ -291,6 +286,11 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & 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) + ! ! Array sizes ! @@ -322,14 +322,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 +341,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 +416,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 !-------------------------------------------------------------------------------------------------------------------- ! @@ -442,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 @@ -498,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) @@ -528,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) 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