Accumulated changes from develop (#37)

RFMIP test cases run completely on GPUs (tested with PGI 16.4). Merged default and OpenACC kernel implementations where possible. Updated DOI in Readme to point to accepted paper. Coefficient files opened read-only (not read-write) Reorganized RFMIP example scripts. Thread-safety by not initializing pointers to NULL(). Closes #35. Specific commits: * Update README.md * Travis CI integration (#24) * Merging Travis CI onto GPU-Hackathon 2019 branch (#25) * Missed one file in last commit * tweaks for CPU compilation * In OpenACC: allocating types as well as data components in ACC copyins, deletes * Ignoring things for Luis... * Missing statement * Aligning array sizes in kernels with arguments * Refined argument intents for some kernels * Workaround for PGI compiler problem with logicals * Input sanitizing gets it own module * Yeah, we'll need the sanitizing module too. * OpenACC-compatible checking for max and min values * OpenACC value checking working; OLCF makefile doesn't use managed memory * Logical kind chosen with pre-preprocessor flags * End results from GPU Hackathon19 (#27) -- Several small bug fixes (argument intent, maximum interpolation indices, thanks to Sebastian Rast) -- Parameterized checking for out-of-range values (works also on GPU) -- Continuous integration with Travis (thanks to Valentin Clement) -- Logical type defaults to Fortran; can be set to use c_bool with -DUSE_CBOOL -- Internal build system can use environmental variables instead of specified files (Makefile.conf etc.) to define compilers, flags, choose kernel directory -- Python scripts to automate running and testing of RFMIP examples -- Update RFMIP examples to use version 1.2 of atmospheres file -- End-to-end RFMIP examples on GPU are broken; fixes pending * Removing unneeded USE statement (thanks to Cheil van Heerwarden). * Remove nullify on declaration for thread safety (#29) Remove nullify statements on declaration of pointers in subroutines to ensure thread safety for mo_gas_optics_rrtmgp. When pointers get assigned in declarations, they implicitly get a save attribute and are assumed static. This is a problem when then occurs in a threaded region, so this code was NOT thread-safe before. Removing the `=> NULL()` does not change the behavior of the code for non-threaded applications, but does ensure thread-safety. * Array size bug fix in compute_bc() * Open coefficients files for read-only access (#32) Open coefficients files for read-only, rather than read-write access because we do NOT want to accidentally write to these files, nor do we want to require users to have write-permissions to load these files. Closes #31, #32. * Moved downloading of reference results for RFMIP from file staging script to comparision * Updating README with DOI for overview paper. * GPU refinement (#34) * Shortwave RFMIP running end-to-end on GPU. Boundary conditions still on CPU. * Upper boundary condition lives on GPU in LW no-scattering calculation. * Moved optical props validation in rte_lw(); simplified data movement in gas optics. Source function still sloshing back and forth between host and device. * Surface emissivity computed on GPU in LW RFMIP example * Moved transposition of surface Planck source onto GPU, clumsily; LW RFMIP cases now running end-to-end on device. * RFMIP boundary conditions on GPU; removing async (may add back later) * Reorder kernels use a single source * Moving array-zeroing routines into mo_util_array * Single-source for array utilities * rte_sw uses array utilities to check validity of boundary conditions * Adding 1D array-zeroing routine * Some SW RFMIP boundary conditions on GPU. * Single-source for fluxes_broadband_kernels * Removing an unneeded OpenACC data transfer * Array value checking uses functions in mo_rte_lw; syntactic cleanup * Correcting mal-formed Makefile * Refined copying of one array in SW examples. * Ben Hillman spots a GPU array being initialized on the CPU. Fixed that. * Updating CSCS compiler and module information as suggesed by Phillipe Marti. Closes #36. * Further updates to Daint modules, library paths from Philippe Marti.
earth-system-radiation · Aug 8, 2019 · ed5b011 · ed5b011
1 parent d7f19ed
commit ed5b011
Show file tree

Hide file tree

Showing 24 changed files with 396 additions and 545 deletions.
diff --git a/README.md b/README.md
@@ -1,6 +1,6 @@
 # RTE+RRTMGP
 
-This is the repository for RTE+RRTMGP, a set of codes for computing radiative fluxes in planetary atmospheres. RTE+RRTMGP is described in a [manuscript revised May 10, 2019](https://doi.org/10.1002/essoar.10500964.1) to [Journal of Advances in Modeling Earth Systems](http://james.agu.org). 
+This is the repository for RTE+RRTMGP, a set of codes for computing radiative fluxes in planetary atmospheres. RTE+RRTMGP is described in a [paper](https://doi.org/10.1029/2019MS001621) in [Journal of Advances in Modeling Earth Systems](http://james.agu.org). 
 
 RRTMGP uses a k-distribution to provide an optical description (absorption and possibly Rayleigh optical depth) of the gaseous atmosphere, along with the relevant source functions, on a pre-determined spectral grid given temperatures, pressures, and gas concentration. The k-distribution currently distributed with this package is applicable to the Earth's atmosphere under present-day, pre-industrial, and 4xCO2 conditions.
 
@@ -11,13 +11,13 @@ Example programs and documenation are evolving - please see examples/ in the rep
 ## Building the libraries.
 
 1. `cd build`
-2. Set environment variables `FC` (the Fortran 2003 compiler) and `FCFLAGS` (compiler flags). Alternately create a Makefile.conf that sets these variables. You could also link to an existing file. 
-3. Set environment variable `RTE_KERNELS` to `openacc` if you want the OpenACC kernels rather than the default. 
+2. Set environment variables `FC` (the Fortran 2003 compiler) and `FCFLAGS` (compiler flags). Alternately create a Makefile.conf that sets these variables. You could also link to an existing file.
+3. Set environment variable `RTE_KERNELS` to `openacc` if you want the OpenACC kernels rather than the default.
 4. `make`
 
 ## Building and running the examples.
 
-1. From the root RTE+RRTMGP directory: `cd examples/rfmip-clear-sky`. 
-2. Set environment variables `NCHOME` and `NFHOME` to the root directories of the Netcdf C and Fortran libraries respectively. Set environment variable `RRTMGP_DIR` to the location of the libraries (`../../build`) in the default layout). 
+1. From the root RTE+RRTMGP directory: `cd examples/rfmip-clear-sky`.
+2. Set environment variables `NCHOME` and `NFHOME` to the root directories of the Netcdf C and Fortran libraries respectively. Set environment variable `RRTMGP_DIR` to the location of the libraries (`../../build`) in the default layout).
 3. `make`
 4. Python scripts are provided to stage the files needed (`stage_files.py`),  run the examples `run-rfmip-examples.py`), and compare to results computed on an example host (`compare-to-reference.py`). The python scripts require modules xarray and netCDF.  
diff --git a/build/Makefile.conf.ifort b/build/Makefile.conf.ifort
@@ -8,8 +8,8 @@ export FC = ifort
 #
 # Optimized
 #
-export FCFLAGS  += -m64 -O3 -traceback -assume realloc_lhs -extend-source 132
-export F77FLAGS += -m64 -O3 -traceback
+export FCFLAGS  += -m64 -O3 -g -traceback -assume realloc_lhs -extend-source 132
+export F77FLAGS += -m64 -O3 -g -traceback
 # can add  -qopt-report-phase=vec
 #
 # Debugging

diff --git a/build/Makefile.conf.pgfortran-cscs b/build/Makefile.conf.pgfortran-cscs
@@ -2,10 +2,11 @@
 
 # Load the following modules to compile with PGI for CPU
 #
-#  module swap PrgEnv-cray PrgEnv-pgi
-#  module swap pgi pgi/18.10.0
-#  module load cray-netcdf cray-hdf5
-#  module unload cray-libsci_acc
+# module load cdt/19.06
+# module swap PrgEnv-cray PrgEnv-pgi
+# module load cray-netcdf cray-hdf5
+# module load craype-accel-nvidia60
+# module unload cray-libsci_acc
 #
 #
 # Fortran compiler command

diff --git a/build/Makefile.conf.pgfortran-cscs-gpu b/build/Makefile.conf.pgfortran-cscs-gpu
@@ -2,10 +2,11 @@
 
 # Load the following modules
 #
-#  module swap PrgEnv-cray PrgEnv-pgi
-#  module load cray-netcdf cray-hdf5
-#  module load craype-accel-nvidia60
-#  module unload cray-libsci_acc
+# module load cdt/19.06
+# module swap PrgEnv-cray PrgEnv-pgi
+# module load cray-netcdf cray-hdf5
+# module load craype-accel-nvidia60
+# module unload cray-libsci_acc
 #
 #
 # Fortran compiler command
@@ -14,7 +15,7 @@ NCHOME =
 
 # Fortran compiler flags
 #FCFLAGS = -g -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk -Mipa=fast,inline -acc -ta=tesla:6.5
-FCFLAGS = -g -ta=tesla:cc60,cuda9.1 -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk                   -Mallocatable=03
+FCFLAGS = -g -ta=tesla:cc60,cuda9.2 -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk -Mallocatable=03
 
 # Fortran .mod files, e.g. -I<include dir> if you have headers in a nonstandard directory <include dir>
 FCINCLUDE =

diff --git a/examples/mo_load_coefficients.F90 b/examples/mo_load_coefficients.F90
@@ -97,7 +97,7 @@ subroutine load_and_init(kdist, filename, available_gases)
     !
     ! How big are the various arrays?
     !
-    if(nf90_open(trim(fileName), NF90_WRITE, ncid) /= NF90_NOERR) &
+    if(nf90_open(trim(fileName), NF90_NOWRITE, ncid) /= NF90_NOERR) &
       call stop_on_err("load_and_init(): can't open file " // trim(fileName))
     ntemps            = get_dim_size(ncid,'temperature')
     npress            = get_dim_size(ncid,'pressure')

diff --git a/examples/rfmip-clear-sky/Makefile.libs.pgfortran-cscs b/examples/rfmip-clear-sky/Makefile.libs.pgfortran-cscs
@@ -1,8 +1,11 @@
-# Load the following modules
+# Load the following modules and set the library path
 #
-#  module swap PrgEnv-cray PrgEnv-pgi
-#  module load cray-netcdf cray-hdf5
-#  module unload cray-libsci_acc
+# module load cdt/19.06
+# module swap PrgEnv-cray PrgEnv-pgi
+# module load cray-netcdf cray-hdf5
+# module load craype-accel-nvidia60
+# module unload cray-libsci_acc
+# export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
 
 export FC = ftn
 export FCFLAGS = -g -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk -Mipa=fast,inline -Mallocatable=03

diff --git a/examples/rfmip-clear-sky/compare-to-reference.py b/examples/rfmip-clear-sky/compare-to-reference.py
@@ -5,12 +5,31 @@
 import os
 import numpy as np
 import xarray as xr
+import urllib.request
 
 ref_dir  = "reference"
 tst_dir = "."
 
-rrtmg_suffix = "_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
+#
+# Download reference results
+#
+print("Downloading reference results")
+ref_dir = "./reference/"
+if not os.path.exists(ref_dir):
+    os.makedirs(ref_dir)
+urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/kbhl3JOSccGtR0m/download", \
+                           os.path.join(ref_dir, "rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
+urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/iFa28GFxRaNGKU1/download", \
+                           os.path.join(ref_dir, "rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
+urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/uCemCHlGxbGK0gJ/download", \
+                           os.path.join(ref_dir, "rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
+urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/l8ZG28j9ttZWD9r/download", \
+                           os.path.join(ref_dir, "rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
 
+#
+# Comparing reference and test results 
+#
+rrtmg_suffix = "_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"
 tst = xr.open_mfdataset(os.path.join(tst_dir, "r??" + rrtmg_suffix))
 ref = xr.open_mfdataset(os.path.join(ref_dir, "r??" + rrtmg_suffix))
 

diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90
@@ -225,10 +225,9 @@ program rrtmgp_rfmip_lw
   ! NOTE: these are causing problems right now, most likely due to a compiler
   ! bug related to the use of Fortran classes on the GPU.
   !
-  !!!$acc enter data create(sfc_emis_spec)
-  !!$acc enter data create(optical_props, optical_props%tau)
-  !!$acc enter data create(source, source%lay_source, source%lev_source_inc, source%lev_source_dec)
-  !!$acc enter data create(source%sfc_source, source%band2gpt, source%gpt2band, source%band_lims_wvn)
+  !$acc enter data create(sfc_emis_spec)
+  !$acc enter data create(optical_props, optical_props%tau)
+  !$acc enter data create(source, source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source)
   ! --------------------------------------------------
 #ifdef USE_TIMING
   !
@@ -251,7 +250,7 @@ program rrtmgp_rfmip_lw
     ! Expand the spectrally-constant surface emissivity to a per-band emissivity for each column
     !   (This is partly to show how to keep work on GPUs using OpenACC)
     !
-    !!$acc parallel loop collapse(2) copyin(sfc_emis)
+    !$acc parallel loop collapse(2) copyin(sfc_emis)
     do icol = 1, block_size
       do ibnd = 1, nbnd
         sfc_emis_spec(ibnd,icol) = sfc_emis(icol,b)
@@ -299,10 +298,10 @@ program rrtmgp_rfmip_lw
   ret = gptlpr(block_size)
   ret = gptlfinalize()
 #endif
-  !!!$acc exit data delete(sfc_emis_spec)
-  !!$acc exit data delete(optical_props%tau, optical_props)
-  !!$acc exit data delete(source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source)
-  !!$acc exit data delete(source%band2gpt, source%gpt2band, source%band_lims_wvn, source)
+  !$acc exit data delete(sfc_emis_spec)
+  !$acc exit data delete(optical_props%tau, optical_props)
+  !$acc exit data delete(source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source)
+  !$acc exit data delete(source)
   ! --------------------------------------------------m
   call unblock_and_write(trim(flxup_file), 'rlu', flux_up)
   call unblock_and_write(trim(flxdn_file), 'rld', flux_dn)

diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90
@@ -46,6 +46,10 @@ program rrtmgp_rfmip_sw
   !
   use mo_rte_kind,           only: wp
   !
+  ! Array utilities
+  !
+  use mo_util_array,         only: zero_array
+  !
   ! Optical properties of the atmosphere as array of values
   !   In the longwave we include only absorption optical depth (_1scl)
   !   Shortwave calculations use optical depth, single-scattering albedo, asymmetry parameter (_2str)
@@ -224,11 +228,9 @@ program rrtmgp_rfmip_sw
            flux_dn(block_size, nlay+1, nblocks))
   allocate(mu0(block_size), sfc_alb_spec(nbnd,block_size))
   call stop_on_err(optical_props%alloc_2str(block_size, nlay, k_dist))
-  ! Handle GPU data. Leave mu0, sfc_alb_spec, toa_flux, and def_tsi on CPU for
-  ! now, and let compiler or CUDA runtime handle data movement because not
-  ! everything is in kernels at the next level down yet.
-  !!$acc enter data create(optical_props, optical_props%tau, optical_props%ssa, optical_props%g)
-  !!!$acc enter data create(mu0,sfc_alb_spec,toa_flux,def_tsi)
+  !$acc enter data create(optical_props, optical_props%tau, optical_props%ssa, optical_props%g)
+  !$acc enter data create (toa_flux, def_tsi)
+  !$acc enter data create (sfc_alb_spec, mu0)
   ! --------------------------------------------------
 #ifdef USE_TIMING
   !
@@ -265,25 +267,27 @@ program rrtmgp_rfmip_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?
-    !  The first two loops could be more expressed more ompactly as def_tsi(1:block_size) = sum(toa_flux, dim=2)
     !
-    !!!$acc parallel loop
-    do icol = 1, block_size
-      def_tsi(icol) = toa_flux(icol, 1)
-    end do
-    !!!$acc parallel loop collapse(2)
+#ifdef _OPENACC
+    call zero_array(block_size, def_tsi)
+    !$acc parallel loop collapse(2) copy(def_tsi) copyin(toa_flux)
     do igpt = 1, ngpt
       do icol = 1, block_size
-        !!$acc atomic update
+        !$acc atomic update
         def_tsi(icol) = def_tsi(icol) + toa_flux(icol, igpt)
       end do
     end do
+#else
+    !
+    ! More compactly...
+    !
+    def_tsi(1:block_size) = sum(toa_flux, dim=2)
+#endif
     !
     ! Normalize incoming solar flux to match RFMIP specification
     !
-    !!!$acc parallel loop collapse(2)
+    !$acc parallel loop collapse(2) copyin(total_solar_irradiance, def_tsi) copy(toa_flux)
     do igpt = 1, ngpt
       do icol = 1, block_size
         toa_flux(icol,igpt) = toa_flux(icol,igpt) * total_solar_irradiance(icol,b)/def_tsi(icol)
@@ -292,7 +296,7 @@ program rrtmgp_rfmip_sw
     !
     ! Expand the spectrally-constant surface albedo to a per-band albedo for each column
     !
-    !!!$acc parallel loop collapse(2)
+    !$acc parallel loop collapse(2) copyin(surface_albedo)
     do icol = 1, block_size
       do ibnd = 1, nbnd
         sfc_alb_spec(ibnd,icol) = surface_albedo(icol,b)
@@ -301,7 +305,7 @@ program rrtmgp_rfmip_sw
     !
     ! Cosine of the solar zenith angle
     !
-    !!!$acc parallel loop
+    !$acc parallel loop copyin(solar_zenith_angle, usecol)
     do icol = 1, block_size
       mu0(icol) = merge(cos(solar_zenith_angle(icol,b)*deg_to_rad), 1._wp, usecol(icol,b))
     end do
@@ -341,8 +345,9 @@ program rrtmgp_rfmip_sw
   ret = gptlpr(block_size)
   ret = gptlfinalize()
 #endif
-  !!$acc exit data delete(optical_props, optical_props%tau, optical_props%ssa, optical_props%g)
-  !!!$acc exit data delete(mu0,sfc_alb_spec,toa_flux,def_tsi)
+  !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g, optical_props)
+  !$acc exit data delete(sfc_alb_spec, mu0)
+  !$acc exit data delete(toa_flux, def_tsi)
   ! --------------------------------------------------
   call unblock_and_write(trim(flxup_file), 'rsu', flux_up)
   call unblock_and_write(trim(flxdn_file), 'rsd', flux_dn)

diff --git a/examples/rfmip-clear-sky/stage_files.py b/examples/rfmip-clear-sky/stage_files.py
@@ -32,21 +32,4 @@
 urllib.request.urlretrieve(conds_url,     conds_file)
 print("Dowloading scripts for generating output templates")
 urllib.request.urlretrieve(templ_scr_url, templ_scr)
-#%run -i generate-output-file-templates.py --source_id RTE-RRTMGP-181204
 subprocess.run(["python3", templ_scr, "--source_id", "RTE-RRTMGP-181204"])
-
-#
-# Reference results
-#
-print("Downloading reference results")
-ref_dir = "./reference/"
-if not os.path.exists(ref_dir):
-    os.makedirs(ref_dir)
-urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/kbhl3JOSccGtR0m/download", \
-                           os.path.join(ref_dir, "rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
-urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/iFa28GFxRaNGKU1/download", \
-                           os.path.join(ref_dir, "rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
-urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/uCemCHlGxbGK0gJ/download", \
-                           os.path.join(ref_dir, "rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
-urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/l8ZG28j9ttZWD9r/download", \
-                           os.path.join(ref_dir, "rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"))
diff --git a/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 b/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90
@@ -17,10 +17,6 @@
 module mo_gas_optics_kernels
   use mo_rte_kind,      only: wp, wl
   implicit none
-
-  interface zero_array
-    module procedure zero_array_3D, zero_array_4D
-  end interface
 contains
   ! --------------------------------------------------------------------------------------
   ! Compute interpolation coefficients
@@ -454,7 +450,8 @@ subroutine gas_optical_depths_minor(ncol,nlay,ngpt,        &
                 tau_minor = 0._wp
                 iflav = gpt_flv(idx_tropo,igpt) ! eta interpolation depends on flavor
                 minor_loc = minor_start + (igpt - gptS) ! add offset to starting point
-                kminor_loc = interpolate2D(fminor(:,:,iflav,icol,ilay), kminor, minor_loc, jeta(:,iflav,icol,ilay), jtemp(icol,ilay))
+                kminor_loc = interpolate2D(fminor(:,:,iflav,icol,ilay), kminor, minor_loc, &
+                                           jeta(:,iflav,icol,ilay), jtemp(icol,ilay))
                 tau_minor = kminor_loc * scaling
                 !$acc atomic update
                 tau(igpt,ilay,icol) = tau(igpt,ilay,icol) + tau_minor
@@ -553,7 +550,7 @@ subroutine compute_Planck_source(                        &
     real(wp) :: planck_function(nbnd,nlay+1,ncol)
     ! -----------------
 
-    !$acc enter data copyin(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,temp_ref_min,totplnk_delta,pfracin,totplnk,gpoint_flavor,one)
+    !$acc enter data copyin(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor,one)
     !$acc enter data create(sfc_src,lay_src,lev_src_inc,lev_src_dec)
     !$acc enter data create(pfrac,planck_function)
 
@@ -636,9 +633,9 @@ subroutine compute_Planck_source(                        &
       end do ! ilay
     end do ! icol
 
-    !$acc exit data delete(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,temp_ref_min,totplnk_delta,pfracin,totplnk,gpoint_flavor,one)
-    !$acc exit data copyout(sfc_src,lay_src,lev_src_inc,lev_src_dec)
+    !$acc exit data delete(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor,one)
     !$acc exit data delete(pfrac,planck_function)
+    !$acc exit data copyout(sfc_src,lay_src,lev_src_inc,lev_src_dec)
 
   end subroutine compute_Planck_source
   ! ----------------------------------------------------------
@@ -784,42 +781,4 @@ subroutine combine_and_reorder_nstr(ncol, nlay, ngpt, nmom, tau_abs, tau_rayleig
     end do
   end subroutine combine_and_reorder_nstr
   ! ----------------------------------------------------------
-  subroutine zero_array_3D(ni, nj, nk, array) bind(C, name="zero_array_3D")
-    integer, intent(in) :: ni, nj, nk
-    real(wp), dimension(ni, nj, nk), intent(out) :: array
-    ! -----------------------
-    integer :: i,j,k
-    ! -----------------------
-    !$acc parallel loop collapse(3) &
-    !$acc&     copyout(array(:ni,:nj,:nk))
-    do k = 1, nk
-      do j = 1, nj
-        do i = 1, ni
-          array(i,j,k) = 0.0_wp
-        end do
-      end do
-    end do
-
-  end subroutine zero_array_3D
-  ! ----------------------------------------------------------
-  subroutine zero_array_4D(ni, nj, nk, nl, array) bind(C, name="zero_array_4D")
-    integer, intent(in) :: ni, nj, nk, nl
-    real(wp), dimension(ni, nj, nk, nl), intent(out) :: array
-    ! -----------------------
-    integer :: i,j,k,l
-    ! -----------------------
-    !$acc parallel loop collapse(4) &
-    !$acc&     copyout(array(:ni,:nj,:nk,:nl))
-    do l = 1, nl
-      do k = 1, nk
-        do j = 1, nj
-          do i = 1, ni
-            array(i,j,k,l) = 0.0_wp
-          end do
-        end do
-      end do
-    end do
-
-  end subroutine zero_array_4D
-  ! ----------------------------------------------------------
 end module mo_gas_optics_kernels