diff --git a/.github/workflows/check-api.yml b/.github/workflows/check-api.yml new file mode 100644 index 000000000..b3681e666 --- /dev/null +++ b/.github/workflows/check-api.yml @@ -0,0 +1,34 @@ +name: Check Fortran API +on: + push: + branches: + - main + - develop + pull_request: + branches-ignore: + - documentation + workflow_dispatch: + + +jobs: + API: + runs-on: ubuntu-22.04 + env: + # Core variables: + FC: gfortran-12 + FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan -g -DRTE_USE_CBOOL" + RRTMGP_ROOT: ${{ github.workspace }} + RTE_KERNELS: extern + steps: + # + # Check out repository under $GITHUB_WORKSPACE + # + - name: Check out code + uses: actions/checkout@v4 + # + # Build libraries + # + - name: Build libraries + run: | + $FC --version + make -j4 libs \ No newline at end of file diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 3cb19a757..138994fd4 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -86,14 +86,14 @@ jobs: # # Build libraries, examples and tests # - - name: Build libraries, examples and tests + - name: Build libraries run: | $FC --version make -j4 libs # # Run examples and tests # - - name: Run examples and tests + - name: Build and run examples and tests run: make -j4 tests # # Compare the results diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 21fe2116e..5ef3cd98b 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -42,7 +42,7 @@ jobs: env: # Core variables: FC: ftn - FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel}} + FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel }} # Make variables: RRTMGP_ROOT: ${{ github.workspace }} RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data @@ -93,21 +93,21 @@ jobs: # # Build libraries, examples and tests # - - name: Build libraries, examples and tests + - name: Build libraries run: | $FC --version make -j8 libs # # Run examples and tests (expect success) # - - name: Run examples and tests (expect success) + - name: Build and run examples and tests (expect success) id: run-success if: matrix.config-name != 'cce-gpu-openmp' run: make -j8 tests # # Run examples and tests (expect failure) # - - name: Run examples and tests (expect failure) + - name: Build and run examples and tests (expect failure) if: steps.run-success.outcome == 'skipped' run: | make -j8 tests && { diff --git a/Makefile b/Makefile index ec6eb75a3..96c1e2edd 100644 --- a/Makefile +++ b/Makefile @@ -1,31 +1,28 @@ # # Top-level Makefile # -.PHONY: libs tests check docs -all: libs tests check docs +.PHONY: libs tests check +all: libs tests check libs: - $(MAKE) -C build - $(MAKE) -C tests - $(MAKE) -C examples/all-sky - $(MAKE) -C examples/rfmip-clear-sky + $(MAKE) -C build $@ tests: + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ check: + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ docs: @cd doc; ./build_documentation.sh clean: $(MAKE) -C build $@ + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ rm -rf public diff --git a/build/Makefile b/build/Makefile index 800bbed3e..efd0341d2 100644 --- a/build/Makefile +++ b/build/Makefile @@ -12,6 +12,7 @@ RRTMGP_KERNEL_DIR = ../rrtmgp-kernels all: librte.a librrtmgp.a \ librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a separate-libs: librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a +libs: all COMPILE = $(FC) $(FCFLAGS) $(FCINCLUDE) -c %.o: %.F90 @@ -22,23 +23,38 @@ include $(RRTMGP_DIR)/Make.depends include $(RTE_KERNEL_DIR)/Make.depends include $(RRTMGP_KERNEL_DIR)/Make.depends -VPATH = $(RTE_DIR):$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) # # If using OpenACC/OpenMP files in *-kernels/accel take precendence # ifeq ($(RTE_KERNELS), accel) - VPATH = $(RTE_DIR):$(RTE_KERNEL_DIR)/accel:$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR)/accel:$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) + VPATH = $(RTE_KERNEL_DIR)/accel:$(RRTMGP_KERNEL_DIR)/accel endif +# +# If using external libraries just compile the interfaces +# +ifeq ($(RTE_KERNELS), extern) + VPATH = $(RTE_KERNEL_DIR)/api:$(RRTMGP_KERNEL_DIR)/api +endif +VPATH += $(RTE_DIR):$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) +# +# Complete library - kernels plus Fortran front end +# librte.a: $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE) ar -rvs librte.a $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE) - +# +# Library with just the kernels... +# librtekernels.a: $(RTE_FORTRAN_KERNELS) ar -rvs librtekernels.a $(RTE_FORTRAN_KERNELS) - +# +# ... and just the Fortran front-end +# librtef.a: $(RTE_FORTRAN_INTERFACE) ar -rvs librtef.a $(RTE_FORTRAN_INTERFACE) - +# +# As with RTE, libraries with Fortran front-end and kernels, separate and combined +# librrtmgp.a: $(RRTMGP_FORTRAN_KERNELS) $(RRTMGP_FORTRAN_INTERFACE) ar -rvs librrtmgp.a $(RRTMGP_FORTRAN_KERNELS) $(RRTMGP_FORTRAN_INTERFACE) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index a3da6beb1..d404181c1 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -46,7 +47,7 @@ mo_load_coefficients.o: mo_simple_netcdf.o 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: +tests: rrtmgp_allsky $(RUN_CMD) bash all_tests.sh check: diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index 0d90cea7a..bad1bc80f 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -5,11 +6,10 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build # # RRTMGP library, module files # -# LDFLAGS += -L$(RRTMGP_BUILD) -# LIBS += -lrrtmgp -lrte +LDFLAGS += -L$(RRTMGP_BUILD) +LIBS += -lrrtmgp -lrte FCINCLUDE += -I$(RRTMGP_BUILD) -# # netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. #FCINCLUDE += -I$(NFHOME)/include @@ -38,11 +38,11 @@ ADDITIONS = mo_simple_netcdf.o mo_rfmip_io.o mo_load_coefficients.o all: rrtmgp_rfmip_lw rrtmgp_rfmip_sw -rrtmgp_rfmip_lw: rrtmgp_rfmip_lw.o $(ADDITIONS) $(RRTMGP_BUILD)/librrtmgp.a $(RRTMGP_BUILD)/librte.a +rrtmgp_rfmip_lw: rrtmgp_rfmip_lw.o $(ADDITIONS) rrtmgp_rfmip_lw.o: rrtmgp_rfmip_lw.F90 $(ADDITIONS) -rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) $(RRTMGP_BUILD)/librrtmgp.a $(RRTMGP_BUILD)/librte.a +rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) rrtmgp_rfmip_sw.o: rrtmgp_rfmip_sw.F90 $(ADDITIONS) @@ -50,7 +50,8 @@ mo_rfmip_io.o: mo_rfmip_io.F90 mo_simple_netcdf.o mo_load_coefficients.o: mo_load_coefficients.F90 mo_simple_netcdf.o -tests: multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc \ +tests: rrtmgp_rfmip_lw rrtmgp_rfmip_sw \ + multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc \ rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc \ rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc $(RUN_CMD) ./rrtmgp_rfmip_lw 8 multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc diff --git a/extensions/mo_fluxes_byband.F90 b/extensions/mo_fluxes_byband.F90 index 587ab42bf..7bf334b66 100644 --- a/extensions/mo_fluxes_byband.F90 +++ b/extensions/mo_fluxes_byband.F90 @@ -156,7 +156,7 @@ end function are_desired_byband ! ! Spectral reduction over all points ! - subroutine sum_byband(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux, byband_flux) bind (C) + subroutine sum_byband(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux, byband_flux) bind(C, name="rte_sum_byband") integer, intent(in ) :: ncol, nlev, ngpt, nbnd integer, dimension(2, nbnd), intent(in ) :: band_lims real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux @@ -181,7 +181,7 @@ end subroutine sum_byband ! ! Net flux: Spectral reduction over all points ! - subroutine net_byband_full(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux_dn, spectral_flux_up, byband_flux_net) bind (C) + subroutine net_byband_full(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux_dn, spectral_flux_up, byband_flux_net) bind(C, name="rte_net_byband_full") integer, intent(in ) :: ncol, nlev, ngpt, nbnd integer, dimension(2, nbnd), intent(in ) :: band_lims real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux_dn, spectral_flux_up diff --git a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 new file mode 100644 index 000000000..9272bb16e --- /dev/null +++ b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 @@ -0,0 +1,245 @@ +module mo_gas_optics_rrtmgp_kernels + use mo_rte_kind, only : wp, wl + use mo_rte_util_array,only : zero_array + implicit none + private + public :: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine interpolation( & + ncol,nlay,ngas,nflav,neta, npres, ntemp, & + flavor, & + press_ref_log, temp_ref,press_ref_log_delta, & + temp_ref_min,temp_ref_delta,press_ref_trop_log, & + vmr_ref, & + play,tlay,col_gas, & + jtemp,fmajor,fminor,col_mix,tropo,jeta,jpress) bind(C, name="rrtmgp_interpolation") + use mo_rte_kind, only : wp, wl + ! input dimensions + integer, intent(in) :: ncol,nlay + !! physical domain size + integer, intent(in) :: ngas,nflav,neta,npres,ntemp + !! k-distribution table dimensions + integer, dimension(2,nflav), intent(in) :: flavor + !! index into vmr_ref of major gases for each flavor + real(wp), dimension(npres), intent(in) :: press_ref_log + !! log of pressure dimension in RRTMGP tables + real(wp), dimension(ntemp), intent(in) :: temp_ref + !! temperature dimension in RRTMGP tables + real(wp), intent(in) :: press_ref_log_delta, & + temp_ref_min, temp_ref_delta, & + press_ref_trop_log + !! constants related to RRTMGP tables + real(wp), dimension(2,0:ngas,ntemp), intent(in) :: vmr_ref + !! reference volume mixing ratios used in compute "binary species parameter" eta + + ! inputs from profile or parent function + real(wp), dimension(ncol,nlay), intent(in) :: play, tlay + !! input pressure (Pa?) and temperature (K) + real(wp), dimension(ncol,nlay,0:ngas), intent(in) :: col_gas + !! input column gas amount - molecules/cm^2 + ! outputs + integer, dimension(ncol,nlay), intent(out) :: jtemp, jpress + !! temperature and pressure interpolation indexes + logical(wl), dimension(ncol,nlay), intent(out) :: tropo + !! use lower (or upper) atmosphere tables + integer, dimension(2, ncol,nlay,nflav), intent(out) :: jeta + !! Index for binary species interpolation +#if !defined(__INTEL_LLVM_COMPILER) && __INTEL_COMPILER >= 2021 + ! A performance-hitting workaround for the vectorization problem reported in + ! https://github.com/earth-system-radiation/rte-rrtmgp/issues/159 + ! The known affected compilers are Intel Fortran Compiler Classic + ! 2021.4, 2021.5 and 2022.1. We do not limit the workaround to these + ! versions because it is not clear when the compiler bug will be fixed, see + ! https://community.intel.com/t5/Intel-Fortran-Compiler/Compiler-vectorization-bug/m-p/1362591. + ! We, however, limit the workaround to the Classic versions only since the + ! problem is not confirmed for the Intel Fortran Compiler oneAPI (a.k.a + ! 'ifx'), which does not mean there is none though. + real(wp), dimension(:, :, :, :), intent(out) :: col_mix +#else + real(wp), dimension(2, ncol,nlay,nflav), intent(out) :: col_mix + !! combination of major species's column amounts (first index is strat/trop) +#endif + real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(out) :: fmajor + !! Interpolation weights in pressure, eta, strat/trop + real(wp), dimension(2,2, ncol,nlay,nflav), intent(out) :: fminor + !! Interpolation fraction in eta, strat/trop + end subroutine interpolation + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_tau_absorption( & + ncol,nlay,nbnd,ngpt, & ! dimensions + ngas,nflav,neta,npres,ntemp, & + nminorlower, nminorklower, & ! number of minor contributors, total num absorption coeffs + nminorupper, nminorkupper, & + idx_h2o, & + gpoint_flavor, & + band_lims_gpt, & + kmajor, & + kminor_lower, & + kminor_upper, & + minor_limits_gpt_lower, & + minor_limits_gpt_upper, & + minor_scales_with_density_lower, & + minor_scales_with_density_upper, & + scale_by_complement_lower, & + scale_by_complement_upper, & + idx_minor_lower, & + idx_minor_upper, & + idx_minor_scaling_lower, & + idx_minor_scaling_upper, & + kminor_start_lower, & + kminor_start_upper, & + tropo, & + col_mix,fmajor,fminor, & + play,tlay,col_gas, & + jeta,jtemp,jpress, & + tau) bind(C, name="rrtmgp_compute_tau_absorption") + ! --------------------- + use mo_rte_kind, only : wp, wl + ! input dimensions + integer, intent(in) :: ncol,nlay,nbnd,ngpt !! array sizes + integer, intent(in) :: ngas,nflav,neta,npres,ntemp !! tables sizes + integer, intent(in) :: nminorlower, nminorklower,nminorupper, nminorkupper + !! table sizes + integer, intent(in) :: idx_h2o !! index of water vapor in col_gas + ! --------------------- + ! inputs from object + integer, dimension(2,ngpt), intent(in) :: gpoint_flavor + !! major gas flavor (pair) by upper/lower, g-point + integer, dimension(2,nbnd), intent(in) :: band_lims_gpt + !! beginning and ending g-point for each band + real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor + !! absorption coefficient table - major gases + real(wp), dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower + !! absorption coefficient table - minor gases, lower atmosphere + real(wp), dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper + !! absorption coefficient table - minor gases, upper atmosphere + integer, dimension(2,nminorlower), intent(in) :: minor_limits_gpt_lower + !! beginning and ending g-point for each minor gas + integer, dimension(2,nminorupper), intent(in) :: minor_limits_gpt_upper + logical(wl), dimension( nminorlower), intent(in) :: minor_scales_with_density_lower + !! generic treatment of minor gases - scales with density (e.g. continuum, collision-induced absorption)? + logical(wl), dimension( nminorupper), intent(in) :: minor_scales_with_density_upper + logical(wl), dimension( nminorlower), intent(in) :: scale_by_complement_lower + !! generic treatment of minor gases - scale by density (e.g. self-continuum) or complement? + logical(wl), dimension( nminorupper), intent(in) :: scale_by_complement_upper + integer, dimension( nminorlower), intent(in) :: idx_minor_lower + !! index of each minor gas in col_gas + integer, dimension( nminorupper), intent(in) :: idx_minor_upper + integer, dimension( nminorlower), intent(in) :: idx_minor_scaling_lower + !! for this minor gas, index of the "scaling gas" in col_gas + integer, dimension( nminorupper), intent(in) :: idx_minor_scaling_upper + integer, dimension( nminorlower), intent(in) :: kminor_start_lower + !! starting g-point index in minor gas absorption table + integer, dimension( nminorupper), intent(in) :: kminor_start_upper + logical(wl), dimension(ncol,nlay), intent(in) :: tropo + !! use upper- or lower-atmospheric tables? + ! --------------------- + ! inputs from profile or parent function + real(wp), dimension(2, ncol,nlay,nflav ), intent(in) :: col_mix + !! combination of major species's column amounts - computed in interpolation() + real(wp), dimension(2,2,2,ncol,nlay,nflav ), intent(in) :: fmajor + !! interpolation weights for major gases - computed in interpolation() + real(wp), dimension(2,2, ncol,nlay,nflav ), intent(in) :: fminor + !! interpolation weights for minor gases - computed in interpolation() + real(wp), dimension( ncol,nlay ), intent(in) :: play, tlay + !! input temperature and pressure + real(wp), dimension( ncol,nlay,0:ngas), intent(in) :: col_gas + !! input column gas amount (molecules/cm^2) + integer, dimension(2, ncol,nlay,nflav ), intent(in) :: jeta + !! interpolation indexes in eta - computed in interpolation() + integer, dimension( ncol,nlay ), intent(in) :: jtemp + !! interpolation indexes in temperature - computed in interpolation() + integer, dimension( ncol,nlay ), intent(in) :: jpress + !! interpolation indexes in pressure - computed in interpolation() + ! --------------------- + ! output - optical depth + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth + end subroutine compute_tau_absorption + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, & + ngas,nflav,neta,npres,ntemp, & + gpoint_flavor,band_lims_gpt, & + krayl, & + idx_h2o, col_dry,col_gas, & + fminor,jeta,tropo,jtemp, & + tau_rayleigh) bind(C, name="rrtmgp_compute_tau_rayleigh") + use mo_rte_kind, only : wp, wl + integer, intent(in ) :: ncol,nlay,nbnd,ngpt + !! input dimensions + integer, intent(in ) :: ngas,nflav,neta,npres,ntemp + !! table dimensions + integer, dimension(2,ngpt), intent(in ) :: gpoint_flavor + !! major gas flavor (pair) by upper/lower, g-point + integer, dimension(2,nbnd), intent(in ) :: band_lims_gpt + !! start and end g-point for each band + real(wp), dimension(ntemp,neta,ngpt,2), intent(in ) :: krayl + !! Rayleigh scattering coefficients + integer, intent(in ) :: idx_h2o + !! index of water vapor in col_gas + real(wp), dimension(ncol,nlay), intent(in ) :: col_dry + !! column amount of dry air + real(wp), dimension(ncol,nlay,0:ngas), intent(in ) :: col_gas + !! input column gas amount (molecules/cm^2) + real(wp), dimension(2,2,ncol,nlay,nflav), intent(in ) :: fminor + !! interpolation weights for major gases - computed in interpolation() + integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta + !! interpolation indexes in eta - computed in interpolation() + logical(wl), dimension(ncol,nlay), intent(in ) :: tropo + !! use upper- or lower-atmospheric tables? + integer, dimension(ncol,nlay), intent(in ) :: jtemp + !! interpolation indexes in temperature - computed in interpolation() + ! outputs + real(wp), dimension(ncol,nlay,ngpt), intent(out) :: tau_rayleigh + !! Rayleigh optical depth + end subroutine compute_tau_rayleigh + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_Planck_source( & + ncol, nlay, nbnd, ngpt, & + nflav, neta, npres, ntemp, nPlanckTemp,& + tlay, tlev, tsfc, sfc_lay, & + fmajor, jeta, tropo, jtemp, jpress, & + gpoint_bands, band_lims_gpt, & + pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, & + sfc_src, lay_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") + use mo_rte_kind, only : wp, wl + integer, intent(in) :: ncol, nlay, nbnd, ngpt + !! input dimensions + integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp + !! table dimensions + real(wp), dimension(ncol,nlay ), intent(in) :: tlay !! temperature at layer centers (K) + real(wp), dimension(ncol,nlay+1), intent(in) :: tlev !! temperature at interfaces (K) + real(wp), dimension(ncol ), intent(in) :: tsfc !! surface temperture + integer, intent(in) :: sfc_lay !! index into surface layer + ! Interpolation variables + real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor + !! interpolation weights for major gases - computed in interpolation() + integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta + !! interpolation indexes in eta - computed in interpolation() + logical(wl), dimension( ncol,nlay), intent(in) :: tropo + !! use upper- or lower-atmospheric tables? + integer, dimension( ncol,nlay), intent(in) :: jtemp, jpress + !! interpolation indexes in temperature and pressure - computed in interpolation() + ! Table-specific + integer, dimension(ngpt), intent(in) :: gpoint_bands !! band to which each g-point belongs + integer, dimension(2, nbnd), intent(in) :: band_lims_gpt !! start and end g-point for each band + real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: pfracin !! Fraction of the Planck function in each g-point + real(wp), intent(in) :: temp_ref_min, totplnk_delta !! interpolation constants + real(wp), dimension(nPlanckTemp,nbnd), intent(in) :: totplnk !! Total Planck function by band at each temperature + integer, dimension(2,ngpt), intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point + + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emssion from the surface + real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src !! Planck emssion from layer centers + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src !! Planck emission at layer boundaries + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac + !! Jacobian (derivative) of the surface Planck source with respect to surface temperature + end subroutine compute_Planck_source + end interface + ! ------------------------------------------------------------------------------------------------------------------ +end module mo_gas_optics_rrtmgp_kernels diff --git a/rrtmgp-kernels/api/rrtmgp_kernels.h b/rrtmgp-kernels/api/rrtmgp_kernels.h new file mode 100644 index 000000000..92c36957c --- /dev/null +++ b/rrtmgp-kernels/api/rrtmgp_kernels.h @@ -0,0 +1,123 @@ +/* This code is part of RRTMGP + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + All right reserved. + +Use and duplication is permitted under the terms of the + BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause + +This header files defines the C bindings for the kernels used in RRTMGP + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ +#pragma once +#include "rte_types.h" + +extern "C" +{ + void rrtmgp_interpolation( + const int& ncol, const int& nlay, + const int& ngas, const int& nflav, const int& neta, + const int& npres, const int& ntemp, + const int* flavor, // (2,nflav) + const Float* press_ref_log, // (npres) + const Float* temp_ref, // (ntemp) + const Float& press_ref_log_delta, + const Float& temp_ref_min, + const Float& temp_ref_delta, + const Float& press_ref_trop_log, + const Float* vmr_ref, //(2,ngas+1,ntemp) + const Float* play, // (ncol,nlay) + const Float* tlay, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + int* jtemp, // [out] (ncol*nlay) + Float* fmajor, // [out] (2,2,2,ncol,nlay,nflav) + Float* fminor, // [out[ (2,2, ncol,nlay,nflav)) + Float* col_mix, // [out] (2, ncol,nlay,nflav) + Bool* tropo, // [out] size (ncol*nlay) + int* jeta,// [out] size (2*ncol*nlay*nflav) + int* jpress // [out] size (ncol*nlay) + ); + + void rrtmgp_compute_tau_absorption( + const int& ncol, const int& nlay, const int& nband, const int& ngpt, + const int& ngas, const int& nflav, const int& neta, + const int& npres, const int& ntemp, + const int& nminorlower, const int& nminorklower, + const int& nminorupper, const int& nminorkupper, + const int& idx_h2o, + const int* gpoint_flavor, // (2,ngpt) + const int* band_lims_gpt, // (2,nbnd) + const Float* kmajor, // (ntemp,neta,npres+1,ngpt) + const Float* kminor_lower, // (ntemp,neta,nminorklower) + const Float* kminor_upper, // (ntemp,neta,nminorkupper) + const int* minor_limits_gpt_lower, // (2,nminorlower) + const int* minor_limits_gpt_upper, // (2,nminorupper) + const Bool* minor_scales_with_density_lower, // ( nminorlower) + const Bool* minor_scales_with_density_upper,// ( nminorupper) + const Bool* scale_by_complement_lower,// ( nminorlower) + const Bool* scale_by_complement_upper,// ( nminorupper) + const int* idx_minor_lower, // ( nminorlower) + const int* idx_minor_upper, // ( nminorupper) + const int* idx_minor_scaling_lower,// ( nminorlower) + const int* idx_minor_scaling_upper,// ( nminorupper) + const int* kminor_start_lower, // ( nminorlower) + const int* kminor_start_upper,// ( nminorupper) + const Bool* tropo, // (ncol,nlay) + const Float* col_mix, // (2, ncol,nlay,nflav ) + const Float* fmajor, // (2,2,2,ncol,nlay,nflav ) + const Float* fminor, // (2,2, ncol,nlay,nflav ) + const Float* play, // (ncol,nlay) + const Float* tlay, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + const int* jeta, // (2, ncol,nlay,nflav ) + const int* jtemp, // (ncol,nlay) + const int* jpress, // (ncol,nlay) + Float* tau // [inout] (ncol,nlay.ngpt) + ); + + void rrtmgp_compute_tau_rayleigh( + const int& ncol, const int& nlay, const int& nband, const int& ngpt, + const int& ngas, const int& nflav, const int& neta, const int& npres, const int& ntemp, + const int* gpoint_flavor, // (2,ngpt) + const int* band_lims_gpt, // (2,nbnd) + const Float* krayl, // (ntemp,neta,ngpt,2) + const int& idx_h2o, + const Float* col_dry, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + const Float* fminor, // (2,2,ncol,nlay,nflav) + const int* jeta, // (2, ncol,nlay,nflav) + const Bool* tropo, // (ncol,nlay) + const int* jtemp, // (ncol,nlay) + Float* tau_rayleigh // [inout] (ncol,nlay.ngpt) + ); + + void rrtmgp_compute_Planck_source( + const int& ncol, const int& nlay, const int& nbnd, const int& ngpt, + const int& nflav, const int& neta, const int& npres, const int& ntemp, + const int& nPlanckTemp, + const Float* tlay, // (ncol,nlay ) + const Float* tlev, // (ncol,nlay+1) + const Float* tsfc, //(ncol ) + const int& sfc_lay, + const Float* fmajor, // (2,2,2,ncol,nlay,nflav) + const int* jeta, // (2, ncol,nlay,nflav) + const Bool* tropo, // ( ncol,nlay) + const int* jtemp, // ( ncol,nlay) + const int* jpress, // ( ncol,nlay) + const int* gpoint_bands, // (ngpt) + const int* band_lims_gpt, // (2, nbnd) + const Float* pfracin, // (ntemp,neta,npres+1,ngpt) + const Float& temp_ref_min, const Float& totplnk_delta, + const Float* totplnk, // (nPlanckTemp,nbnd) + const int* gpoint_flavor, // (2,ngpt) + Float* sfc_src, // [out] (ncol, ngpt) + Float* lay_src, // [out] (ncol,nlay, ngpt) + Float* lev_src, // [out] (ncol,nlay+1,ngpt) + Float* sfc_src_jac // [out] (ncol, ngpt) + ); +} diff --git a/rte-kernels/api/mo_fluxes_broadband_kernels.F90 b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 new file mode 100644 index 000000000..d64c1fa29 --- /dev/null +++ b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 @@ -0,0 +1,74 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, Trustees of Columbia University. All right reserved. +! +! Use and duplication is permitted under the terms of the +! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause +! ------------------------------------------------------------------------------------------------- +!> +!> ## Kernels for computing broadband fluxes +!> +! ------------------------------------------------------------------------------------------------- +module mo_fluxes_broadband_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp + implicit none + private + public :: sum_broadband, net_broadband + + ! ---------------------------------------------------------------------------- + !> + !> Spectral reduction over all points + !> + interface + subroutine sum_broadband(ncol, nlev, ngpt, spectral_flux, broadband_flux) bind(C, name="rte_sum_broadband") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev, ngpt + !! Array sizes + real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux + !! Spectrally-resolved flux + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux + !! Sum of spectrally-resolved flux over `ngpt` + end subroutine sum_broadband + end interface + ! ---------------------------------------------------------------------------- + !> + !> Spectral reduction over all points for net flux + !> Overloaded - which routine is called depends on arguments + !> + interface net_broadband + ! ---------------------------------------------------------------------------- + !> + !> Net flux from g-point fluxes up and down + !> + subroutine net_broadband_full(ncol, nlev, ngpt, spectral_flux_dn, spectral_flux_up, broadband_flux_net) & + bind(C, name="rte_net_broadband_full") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev, ngpt + !! Array sizes + real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux_dn, spectral_flux_up + !! Spectrally-resolved flux up and down + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux_net + !! Net (down minus up) summed over `ngpt` + end subroutine net_broadband_full + ! ---------------------------------------------------------------------------- + !> + !> Net flux when bradband flux up and down are already available + !> + subroutine net_broadband_precalc(ncol, nlev, flux_dn, flux_up, broadband_flux_net) & + bind(C, name="rte_net_broadband_precalc") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev + !! Array sizes + real(wp), dimension(ncol, nlev), intent(in ) :: flux_dn, flux_up + !! Broadband downward and upward fluxes + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux_net + !! Net (down minus up) + end subroutine net_broadband_precalc + end interface net_broadband + ! ---------------------------------------------------------------------------- +end module mo_fluxes_broadband_kernels diff --git a/rte-kernels/api/mo_optical_props_kernels.F90 b/rte-kernels/api/mo_optical_props_kernels.F90 new file mode 100644 index 000000000..bb7e5116f --- /dev/null +++ b/rte-kernels/api/mo_optical_props_kernels.F90 @@ -0,0 +1,377 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, Trustees of Columbia University. All right reserved. +! +! Use and duplication is permitted under the terms of the +! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause +! ------------------------------------------------------------------------------------------------- +! +!> ## Kernels for arrays of optical properties: +!> - delta-scaling +!> - adding two sets of properties +!> - extracting subsets along the column dimension +! +! ------------------------------------------------------------------------------------------------- + +module mo_optical_props_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp, wl + implicit none + + public + + ! ------------------------------------------------------------------------------------------------- + ! + ! Delta-scaling is provided only for two-stream properties at present + ! + interface delta_scale_2str_kernel + ! ------------------------------------------------------------------------------------------------- + !> Delta-scale two-stream optical properties given user-provided value of \(f\) (forward scattering) + ! + pure subroutine delta_scale_2str_f_k(ncol, nlay, ngpt, tau, ssa, g, f) & + bind(C, name="rte_delta_scale_2str_f_k") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Array sizes + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g + !! Optical depth, single-scattering albedo, asymmetry parameter + real(wp), dimension(ncol, nlay, ngpt), intent(in ) :: f + !! User-provided forward-scattering fraction + end subroutine delta_scale_2str_f_k + ! --------------------------------- + !> Delta-scale assuming forward-scatternig fraction is the square of the asymmetry parameter + !> i.e. \(f = g^2\) + ! + pure subroutine delta_scale_2str_k(ncol, nlay, ngpt, tau, ssa, g) & + bind(C, name="rte_delta_scale_2str_k") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Array sizes + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g + !! Optical depth, single-scattering albedo, asymmetry parameter + end subroutine delta_scale_2str_k + end interface delta_scale_2str_kernel + ! ------------------------------------------------------------------------------------------------- + ! + ! Addition of optical properties: the first set are incremented by the second set. + ! + ! There are three possible representations of optical properties (scalar = optical depth only; + ! two-stream = tau, single-scattering albedo, and asymmetry factor g, and + ! n-stream = tau, ssa, and phase function moments p.) Thus we need nine routines, three for + ! each choice of representation on the left hand side times three representations of the + ! optical properties to be added. + ! + ! There are two sets of these nine routines. In the first the two sets of optical + ! properties are defined at the same spectral resolution. There is also a set of routines + ! to add properties defined at lower spectral resolution to a set defined at higher spectral + ! resolution (adding properties defined by band to those defined by g-point) + ! + ! ------------------------------------------------------------------------------------------------- + !> increase one absorption optical depth by a second value + interface + pure subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, & + tau1, & + tau2) bind(C, name="rte_increment_1scalar_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_1scalar_by_1scalar + end interface + ! --------------------------------- + !> increase absorption optical depth with extinction optical depth (2-stream form) + interface + pure subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2) bind(C, name="rte_increment_1scalar_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + end subroutine increment_1scalar_by_2stream + end interface + ! --------------------------------- + !> increase absorption optical depth with extinction optical depth (n-stream form) + interface + pure subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2) bind(C, name="rte_increment_1scalar_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + end subroutine increment_1scalar_by_nstream + end interface + ! --------------------------------- + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with absorption optical depth + interface + pure subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2) bind(C, name="rte_increment_2stream_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_2stream_by_1scalar + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with a second set + interface + pure subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, & + tau1, ssa1, g1, & + tau2, ssa2, g2) bind(C, name="rte_increment_2stream_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original + end subroutine increment_2stream_by_2stream + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with _n_-stream + interface + pure subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, & + tau1, ssa1, g1, & + tau2, ssa2, p2) bind(C, name="rte_increment_2stream_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom2 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + real(wp), dimension(nmom2, & + ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added + end subroutine increment_2stream_by_nstream + end interface + ! --------------------------------- + ! --------------------------------- + !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with absorption optical depth + interface + pure subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2) bind(C, name="rte_increment_nstream_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_nstream_by_1scalar + end interface + ! --------------------------------- + !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with two-stream values + interface + pure subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, nmom1, & + tau1, ssa1, p1, & + tau2, ssa2, g2) bind(C, name="rte_increment_nstream_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original + end subroutine increment_nstream_by_2stream + end interface + ! --------------------------------- + !> increment one set of _n_-stream optical properties with another set + interface + pure subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, & + tau1, ssa1, p1, & + tau2, ssa2, p2) bind(C, name="rte_increment_nstream_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + real(wp), dimension(nmom2, & + ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added + end subroutine increment_nstream_by_nstream + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Incrementing when the second set of optical properties is defined at lower spectral resolution + ! (e.g. by band instead of by gpoint) + ! + ! ------------------------------------------------------------------------------------------------- + !> increase one absorption optical depth defined on g-points by a second value defined on bands + interface + pure subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increase absorption optical depth defined on g-points with extinction optical depth (2-stream form) defined on bands + interface + pure subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_2stream_bybnd + end interface + ! --------------------------------- + !> increase absorption optical depth defined on g-points with extinction optical depth (n-stream form) defined on bands + interface + pure subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_nstream_bybnd + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) defined on g-points with absorption optical depth defined on bands + interface + pure subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increment 2-stream optical properties defined on g-points with another set defined on bands + interface + pure subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, g1, & + tau2, ssa2, g2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_2stream_bybnd + end interface + ! --------------------------------- + !> increment 2-stream optical properties defined on g-points with _n_-stream properties set defined on bands + interface + pure subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, nmom2, & + tau1, ssa1, g1, & + tau2, ssa2, p2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom2, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + real(wp), dimension(nmom2, & + ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_nstream_bybnd + end interface + ! --------------------------------- + ! --------------------------------- + !> increment _n_-stream optical properties defined on g-points with absorption optical depth defined on bands + interface + pure subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increment n-stream optical properties defined on g-points with 2-stream properties set defined on bands + interface + pure subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, & + tau1, ssa1, p1, & + tau2, ssa2, g2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_2stream_bybnd + end interface + ! --------------------------------- + !> increment _n_-stream optical properties defined on g-points with a second set defined on bands + interface + pure subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, nmom2, & + tau1, ssa1, p1, & + tau2, ssa2, p2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + real(wp), dimension(nmom2, & + ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_nstream_bybnd + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Subsetting, meaning extracting some portion of the 3D domain + ! + ! ------------------------------------------------------------------------------------------------- + !> + !> Extract a subset from the first dimension (normally columns) of a 3D field. + !> Applicable to most variables e.g. tau, ssa, g + !> + interface extract_subset + pure subroutine extract_subset_dim1_3d(ncol, nlay, ngpt, array_in, colS, colE, array_out) & + bind (C, name="rte_extract_subset_dim1_3d") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(colE-colS+1,& + nlay,ngpt), intent(out) :: array_out !! subset of the input array + end subroutine extract_subset_dim1_3d + ! --------------------------------- + !> Extract a subset from the second dimension (normally columns) of a 4D field. + !> Applicable to phase function moments, where the first dimension is the moment + pure subroutine extract_subset_dim2_4d(nmom, ncol, nlay, ngpt, array_in, colS, colE, array_out) & + bind (C, name="rte_extract_subset_dim2_4d") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: nmom, ncol, nlay, ngpt !! Array sizes + real(wp), dimension(nmom,ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(nmom,colE-colS+1,& + nlay,ngpt), intent(out) :: array_out !! subset of the input array + end subroutine extract_subset_dim2_4d + ! --------------------------------- + ! + !> Extract the absorption optical thickness \(\tau_{abs} = 1 - \omega_0 \tau_{ext}\) + ! + pure subroutine extract_subset_absorption_tau(ncol, nlay, ngpt, tau_in, ssa_in, & + colS, colE, tau_out) & + bind (C, name="rte_extract_subset_absorption_tau") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau_in, ssa_in !! Optical thickness, single scattering albedo + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(colE-colS+1,& + nlay,ngpt), intent(out) :: tau_out !! absorption optical thickness subset + end subroutine extract_subset_absorption_tau + end interface extract_subset +end module mo_optical_props_kernels diff --git a/rte-kernels/api/mo_rte_solver_kernels.F90 b/rte-kernels/api/mo_rte_solver_kernels.F90 new file mode 100644 index 000000000..da8d0f706 --- /dev/null +++ b/rte-kernels/api/mo_rte_solver_kernels.F90 @@ -0,0 +1,202 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, Trustees of Columbia University. All right reserved. +! +! Use and duplication is permitted under the terms of the +! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause +! ------------------------------------------------------------------------------------------------- +! +!>## Numeric calculations for radiative transfer solvers +!> - Emission/absorption (no-scattering) calculations +!> - solver for multi-angle Gaussian quadrature +!> - Extinction-only calculation (direct solar beam) +!> - Two-stream calculations: +!> solvers for LW and SW with different boundary conditions and source functions +! +! ------------------------------------------------------------------------------------------------- +module mo_rte_solver_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp, wl + implicit none + private + + public :: lw_solver_noscat, lw_solver_2stream, & + sw_solver_noscat, sw_solver_2stream + ! ------------------------------------------------------------------------------------------------- + ! + ! Top-level longwave kernels + ! + ! ------------------------------------------------------------------------------------------------- + ! + !> LW transport, no scattering, multi-angle quadrature + !> Users provide a set of weights and quadrature angles + !> Routine sums over single-angle solutions for each sets of angles/weights + ! + ! --------------------------------------------------------------- + interface + subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & + nmus, Ds, weights, & + tau, & + lay_source, lev_source, & + sfc_emis, sfc_src, & + inc_flux, & + flux_up, flux_dn, & + do_broadband, broadband_up, broadband_dn, & + do_Jacobians, sfc_srcJac, flux_upJac, & + do_rescaling, ssa, g) bind(C, name="rte_lw_solver_noscat") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + integer, intent(in ) :: nmus + !! number of quadrature angles + real(wp), dimension (ncol, ngpt, & + nmus), intent(in ) :: Ds + !! quadrature secants + real(wp), dimension(nmus), intent(in ) :: weights + !! quadrature weights + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau + !! Absorption optical thickness [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source + !! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source + !! Planck source at layer edge for radiation [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis + !! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src + !! Surface source function [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux + !! Incident diffuse flux, probably 0 [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), target, & + intent( out) :: flux_up, flux_dn + !! Fluxes [W/m2] + ! + ! Optional variables - arrays aren't referenced if corresponding logical == False + ! + logical(wl), intent(in ) :: do_broadband + real(wp), dimension(ncol,nlay+1 ), target, & + intent( out) :: broadband_up, broadband_dn + !! Spectrally-integrated fluxes [W/m2] + logical(wl), intent(in ) :: do_Jacobians + !! compute Jacobian with respect to surface temeprature? + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac + !! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1 ), target, & + intent( out) :: flux_upJac + !! surface temperature Jacobian of Radiances [W/m2-str / K] + logical(wl), intent(in ) :: do_rescaling + !! Approximate treatment of scattering (10.1175/JAS-D-18-0014.1) + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: ssa, g + !! single-scattering albedo, asymmetry parameter + end subroutine lw_solver_noscat + end interface + ! ------------------------------------------------------------------------------------------------- + ! + !> Longwave two-stream calculation: + !> - combine RRTMGP-specific sources at levels + !> - compute layer reflectance, transmittance + !> - compute total source function at levels using linear-in-tau + !> - transport + ! + ! ------------------------------------------------------------------------------------------------- + interface + subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & + tau, ssa, g, & + lay_source, lev_source, sfc_emis, sfc_src, & + inc_flux, & + flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g + !! Optical thickness, single-scattering albedo, asymmetry parameter [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source + !! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source + !! Planck source at layer edge for radiation [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis + !! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src + !! Surface source function [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux + !! Incident diffuse flux, probably 0 [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up, flux_dn + !! Fluxes [W/m2] + end subroutine lw_solver_2stream + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Top-level shortwave kernels + ! + ! ------------------------------------------------------------------------------------------------- + ! + ! !> Extinction-only shortwave solver i.e. solar direct beam + ! + ! ------------------------------------------------------------------------------------------------- + interface + pure subroutine sw_solver_noscat(ncol, nlay, ngpt, top_at_1, & + tau, mu0, inc_flux_dir, flux_dir) bind(C, name="rte_sw_solver_noscat") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau + !! Absorption optical thickness [] + real(wp), dimension(ncol,nlay ), intent(in ) :: mu0 + !! cosine of solar zenith angle + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir + !! Direct beam incident flux [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dir + end subroutine sw_solver_noscat + end interface + ! ------------------------------------------------------------------------------------------------- + ! + !> Shortwave two-stream calculation: + !> compute layer reflectance, transmittance + !> compute solar source function for diffuse radiation + !> transport + ! + ! ------------------------------------------------------------------------------------------------- + interface + subroutine sw_solver_2stream (ncol, nlay, ngpt, top_at_1, & + tau, ssa, g, mu0, & + sfc_alb_dir, sfc_alb_dif, & + inc_flux_dir, & + flux_up, flux_dn, flux_dir, & + has_dif_bc, inc_flux_dif, & + do_broadband, broadband_up, & + broadband_dn, broadband_dir) bind(C, name="rte_sw_solver_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g + !! Optical thickness, single-scattering albedo, asymmetry parameter [] + real(wp), dimension(ncol,nlay ), intent(in ) :: mu0 + !! cosine of solar zenith angle + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_alb_dir, sfc_alb_dif + !! Spectral surface albedo for direct and diffuse radiation + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir + !! Direct beam incident flux + real(wp), dimension(ncol,nlay+1,ngpt), target, & + intent( out) :: flux_up, flux_dn, flux_dir + !! Fluxes [W/m2] + logical(wl), intent(in ) :: has_dif_bc + !! Is a boundary condition for diffuse flux supplied? + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dif + !! Boundary condition for diffuse flux [W/m2] + logical(wl), intent(in ) :: do_broadband + !! Provide broadband-integrated, not spectrally-resolved, fluxes? + real(wp), dimension(ncol,nlay+1 ), intent( out) :: broadband_up, broadband_dn, broadband_dir + end subroutine sw_solver_2stream + end interface +end module mo_rte_solver_kernels diff --git a/rte-kernels/api/mo_rte_util_array.F90 b/rte-kernels/api/mo_rte_util_array.F90 new file mode 100644 index 000000000..cdae473c4 --- /dev/null +++ b/rte-kernels/api/mo_rte_util_array.F90 @@ -0,0 +1,47 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! All right reserved. +! +! Use and duplication is permitted under the terms of the +! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause +! ------------------------------------------------------------------------------------------------- +module mo_rte_util_array + use mo_rte_kind, only: wp, wl + implicit none + public :: zero_array + + !------------------------------------------------------------------------------------------------- + ! Initializing arrays to 0 + !------------------------------------------------------------------------------------------------- + interface zero_array + subroutine zero_array_1D(ni, array) bind(C, name="zero_array_1D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni + real(wp), dimension(ni), intent(out) :: array + end subroutine zero_array_1D + ! ---------------------------------------------------------- + subroutine zero_array_2D(ni, nj, array) bind(C, name="zero_array_2D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj + real(wp), dimension(ni, nj), intent(out) :: array + end subroutine zero_array_2D + ! ---------------------------------------------------------- + subroutine zero_array_3D(ni, nj, nk, array) bind(C, name="zero_array_3D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj, nk + real(wp), dimension(ni, nj, nk), intent(out) :: array + end subroutine zero_array_3D + ! ---------------------------------------------------------- + subroutine zero_array_4D(ni, nj, nk, nl, array) bind(C, name="zero_array_4D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj, nk, nl + real(wp), dimension(ni, nj, nk, nl), intent(out) :: array + end subroutine zero_array_4D + end interface zero_array +end module mo_rte_util_array diff --git a/rte-kernels/api/rte_kernels.h b/rte-kernels/api/rte_kernels.h new file mode 100644 index 000000000..71f86b81e --- /dev/null +++ b/rte-kernels/api/rte_kernels.h @@ -0,0 +1,389 @@ +/* This code is part of Radiative Transfer for Energetics (RTE) + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + All right reserved. + +Use and duplication is permitted under the terms of the + BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause + +This header files defines the C bindings for the kernels used in RTE + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ + +#include "rte_types.h" + +extern "C" +{ + // + // Shortwave solvers + // + void rte_sw_solver_noscat( + const int& ncol, const int& nlay, const int& ngpt, const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* mu0, // (ncol,nlay) + const Float* inc_flux_dir, // (ncol, ngpt) + Float* flux_dir); // [out] (ncol,nlay+1,ngpt) + + void rte_sw_solver_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g, // (ncol,nlay, ngpt) + const Float* mu0, // (ncol,nlay) + const Float* sfc_alb_dir, // (ncol, ngpt) + const Float* sfc_alb_dif, // (ncol, ngpt) + const Float* inc_flux_dir, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn, // [out] (ncol,nlay+1,ngpt) + Float* flux_dir, // [out] (ncol,nlay+1,ngpt) + const Bool& has_dif_bc, + const Float* inc_flux_dif, // (ncol, ngpt) + const Bool& do_broadband, + Float* broadband_up, // [out] (ncol,nlay+1) + Float* broadband_dn, // [out] (ncol,nlay+1) + Float* broadband_dir); // [out] (ncol,nlay+1) + + void rte_lw_solver_noscat( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const int& nmus, + const Float* secants, // (nmus) + const Float* weights, // (nmus) + const Float* tau, // (ncol,nlay, ngpt) + const Float* lay_source, // (ncol,nlay, ngpt) + const Float* lev_source, // (ncol,nlay+1,ngpt) + const Float* sfc_emis, // (ncol, ngpt) + const Float* sfc_src, // (ncol, ngpt) + const Float* inc_flux, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn, // [out] (ncol,nlay+1,ngpt) + const Bool& do_broadband, + Float* broadband_up, + // [out] (ncol,nlay+1) + Float* broadband_dn, + // [out] (ncol,nlay+1) + const Bool& do_jacobians, + const Float* sfc_src_jac, + // (ncol, ngpt) + Float* flux_up_jac, + // [out] (ncol,nlay+1,ngpt) + const Bool& do_rescaling, + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g); // (ncol,nlay, ngpt) + + + void rte_lw_solver_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g, // (ncol,nlay, ngpt) + const Float* lay_source, // (ncol,nlay, ngpt) + const Float* lev_source, // (ncol,nlay+1,ngpt) + const Float* sfc_emis, // (ncol, ngpt) + const Float* sfc_src, // (ncol, ngpt) + const Float* inc_flux, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn // [out] (ncol,nlay+1,ngpt) + ); + + // + // OPTICAL PROPS - INCREMENT + // + void rte_increment_1scalar_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in); // (ncol,nlay,ngpt) + + + void rte_increment_1scalar_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_1scalar_by_nstream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_2stream_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in); // (ncol,nlay,ngpt) + + void rte_increment_2stream_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* g_in); // (ncol,nlay,ngpt) + + + void rte_increment_2stream_by_nstream( + const int& ncol, const int& nlay, const int& ngpt, const int& nmom, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* p_in); //(nmom,ncol,nlay,ngpt) + + void rte_increment_nstream_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_nstream_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, // [inout] (nmom,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* g_in); // (ncol,nlay,ngpt) + + void rte_increment_nstream_by_nstream( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + const int& nmom2, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, // [inout](nmom1,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* p_in); // (nmom2,ncol,nlay,ngpt) + + // + // OPTICAL PROPS - INCREMENT BYBND + // + void rte_inc_1scalar_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout,// [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_1scalar_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_1scalar_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* g_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* p_in, // (nmom,ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, + // [inout] (nomo,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* g_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + const int& nmom2, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, + // [inout] (nomo,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* p_in, // (nmom,ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + // + // OPTICAL PROPS - DELTA SCALING + // + void rte_delta_scale_2str_k( + const int& ncol, const int& nlay, const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlev,ngpt) + Float* ssa_inout, // [inout] (ncol,nlev,ngpt) + Float* g_inout); // [inout] (ncol,nlev,ngpt) + + void rte_delta_scale_2str_f_k( + const int& ncol, const int& nlay, const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlev,ngpt) + Float* ssa_inout, // [inout] (ncol,nlev,ngpt) + Float* g_inout, // [inout] (ncol,nlev,ngpt) + const Float* f); // (ncol,nlev,ngpt) + + // + // OPTICAL PROPS - SUBSET + // + void rte_extract_subset_dim1_3d( + const int& ncol, const int& nlay, const int& ngpt, + Float* array_in, // (ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* array_out); // [out] (ncol_end-ncol_start+1,nlay,ngpt) + + void rte_extract_subset_dim2_4d( + const int& nmom, const int& ncol, const int& nlay, const int& ngpt, + const Float* array_in, // (nmom,ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* array_out); // [out] (nmom,ncol_end-ncol_start+1,nlay,ngpt) + + void rte_extract_subset_absorption_tau( + const int& ncol, const int& nlay, const int& ngpt, + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* tau_out); // [out] (ncol_end-ncol_start+1,nlay,ngpt) + + // + // Fluxes - reduction + // + void rte_sum_broadband( + const int& ncol, const int& nlev, const int& ngpt, + const Float* gpt_flux, // (ncol,nlev,ngpt) + Float* flux); // [out] (ncol,nlev) + + void rte_net_broadband_full( + const int& ncol, const int& nlev, const int& ngpt, + const Float* gpt_flux_dn, // (ncol,nlev,ngpt) + const Float* gpt_flux_up, // (ncol,nlev,ngpt) + Float* flux_net); // [out] (ncol,nlev) + + void rte_net_broadband_precalc( + const int& ncol, const int& nlev, + const Float* broadband_flux_dn, // (ncol, nlev) + const Float* broadband_flux_up, // (ncol, nlev) + Float* broadband_flux_net);//[out] (ncol, nlev) + + void rte_sum_byband( + const int& ncol, const int& nlev, const int& ngpt, const int& nbnd, + const int* band_lims, // dimension(2, nbnd) + const Float* gpt_flux, // (ncol, nlev, ngpt) + Float* bnd_flux); // [out] (ncol, nlev, ngpt) + + void rte_net_byband_full( + const int& ncol, const int& nlev, const int& ngpt, const int& nbnd, int* band_lims, + const Float* bnd_flux_dn, // (ncol,nlev,nbnd) + const Float* bnd_flux_up, // (ncol,nlev,nbnd) + Float* bnd_flux_net); // [out] (ncol,nlev) + // + // Array utilities + // + void zero_array_1D( + const int& ni, + Float* array); // [out] (ni) + + void zero_array_2D( + const int& ni, const int& nj, + Float* array); // [out] (ni, nj) + + void zero_array_3D( + const int& ni, const int& nj, const int& nk, + Float* array); // [out] (ni, nj, nk) + + void zero_array_4D( + const int& ni, const int& nj, const int& nk, const int& nl, + Float* array); // [out] (ni, nj, nk, nl) + +} diff --git a/rte-kernels/api/rte_types.h b/rte-kernels/api/rte_types.h new file mode 100644 index 000000000..fc645f203 --- /dev/null +++ b/rte-kernels/api/rte_types.h @@ -0,0 +1,34 @@ +/* This code is part of Radiative Transfer for Energetics (RTE) + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + All right reserved. + +Use and duplication is permitted under the terms of the + BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause + +This header files C-compatible Boolean and floating point types (see mo_rte_type.F90 for the Fortran equivalent) + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ +#pragma once + +#ifdef RTE_USE_CBOOL +using Bool = signed char; +#else +using Bool = int; +#endif + +#ifdef RTE_USE_SP +using Float = float; +const Float Float_epsilon = FLT_EPSILON; +#else +using Float = double; +const Float Float_epsilon = DBL_EPSILON; +#endif + + + diff --git a/tests/Makefile b/tests/Makefile index 1b90ed26c..a5060c623 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -5,20 +6,18 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build # # RRTMGP library, module files # -# LDFLAGS += -L$(RRTMGP_BUILD) -# LIBS += -lrrtmgp -lrte +LDFLAGS += -L$(RRTMGP_BUILD) +LIBS += -lrrtmgp -lrte FCINCLUDE += -I$(RRTMGP_BUILD) -# + # netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. #FCINCLUDE += -I$(NFHOME)/include # netcdf C and Fortran libraries have to be in the search path or added via environment variable LDFLAGS e.g. #LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib -LDFLAGS += -L$(RRTMGP_BUILD) -LIBS += -lrte -lrrtmgp LIBS += -lnetcdff -lnetcdf -VPATH = .:$(RRTMGP_ROOT)/examples:$(RRTMGP_ROOT)/examples/rfmip-clear-sky:$(RRTMGP_ROOT)/examples/all-sky +VPATH = $(RRTMGP_ROOT)/examples:$(RRTMGP_ROOT)/examples/rfmip-clear-sky:$(RRTMGP_ROOT)/examples/all-sky VPATH += $(RRTMGP_ROOT)/rrtmgp-frontend:$(RRTMGP_ROOT)/extensions:$(RRTMGP_ROOT)/:$(RRTMGP_ROOT)/extensions/solar_variability # Compilation rules @@ -80,7 +79,7 @@ rte_sw_solver_unit_tests : $(LIB_DEPS) mo_testing_utils.o rte_sw_solver_unit_te .PHONY: tests -tests: +tests: check_variants check_equivalence test_zenith_angle_spherical_correction rte_sw_solver_unit_tests rte_optic_prop_unit_tests rte_lw_solver_unit_tests cp ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc ./test_atmospheres.nc $(RUN_CMD) bash all_tests.sh