WIP: compiles and runs tests; clear-sky pasess but all-sky not

examples/all-sky/rrtmgp_allsky.F90

@@ -89,7 +89,7 @@ program rte_rrtmgp_allsky
   !
   ! Inputs to RRTMGP
   !
-  logical :: top_at_1, is_sw, is_lw
+  logical :: is_sw, is_lw
 
   integer  :: nbnd, ngpt
   integer  :: icol, ilay, ibnd, iloop, igas
@@ -191,7 +191,6 @@ program rte_rrtmgp_allsky
   !
   nbnd = k_dist%get_nband()
   ngpt = k_dist%get_ngpt()
-  top_at_1 = p_lay(1, 1) < p_lay(1, nlay)
 
   ! ----------------------------------------------------------------------------
   ! LW calculations neglect scattering; SW calculations use the 2-stream approximation
@@ -251,7 +250,7 @@ program rte_rrtmgp_allsky
     ! Surface temperature
     !$acc kernels
     !$omp target
-    t_sfc = t_lev(1, merge(nlay+1, 1, top_at_1))
+    t_sfc = t_lev(1, merge(nlay+1, 1, atmos%top_is_at_1()))
     emis_sfc = 0.98_wp
     !$acc end kernels
     !$omp end target
@@ -322,9 +321,9 @@ program rte_rrtmgp_allsky
                                          tlev = t_lev))
       if(do_clouds)   call stop_on_err(clouds%increment(atmos))
       if(do_aerosols) call stop_on_err(aerosols%increment(atmos))
-      call stop_on_err(rte_lw(atmos, top_at_1, &
-                              lw_sources,      &
-                              emis_sfc,        &
+      call stop_on_err(rte_lw(atmos,      &
+                              lw_sources, &
+                              emis_sfc,   &
                               fluxes))
       !$acc        end data
       !$omp end target data
@@ -347,8 +346,7 @@ program rte_rrtmgp_allsky
         call stop_on_err(aerosols%delta_scale())
         call stop_on_err(aerosols%increment(atmos))
       end if 
-      call stop_on_err(rte_sw(atmos, top_at_1, &
-                              mu0,   toa_flux, &
+      call stop_on_err(rte_sw(atmos, mu0,   toa_flux, &
                               sfc_alb_dir, sfc_alb_dif, &
                               fluxes))
       !$acc        end data   

examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90

@@ -90,7 +90,6 @@ program rrtmgp_rfmip_lw
   character(len=132) :: kdist_file = 'coefficients_lw.nc'
   character(len=132) :: flxdn_file, flxup_file
   integer            :: nargs, ncol, nlay, nbnd, nexp, nblocks, block_size, forcing_index, physics_index, n_quad_angles = 1
-  logical            :: top_at_1
   integer            :: b, icol, ibnd
   character(len=4)   :: block_size_char, forcing_index_char = '1', physics_index_char = '1'
 
@@ -168,10 +167,6 @@ program rrtmgp_rfmip_lw
   ! Allocation on assignment within reading routines
   !
   call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev)
-  !
-  ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain?
-  !
-  top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1)
 
   !
   ! Read the gas concentrations and surface properties
@@ -193,7 +188,8 @@ program rrtmgp_rfmip_lw
   !   is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep.
   !   This introduces an error but shows input sanitizing.
   !
-  if(top_at_1) then
+  ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain?
+  if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then
     p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min())
   else
     p_lev(:,nlay+1,:) &
@@ -256,7 +252,6 @@ program rrtmgp_rfmip_lw
     !    via ty_fluxes_broadband
     !
     call stop_on_err(rte_lw(optical_props,   &
-                            top_at_1,        &
                             source,          &
                             sfc_emis_spec,   &
                             fluxes, n_gauss_angles = n_quad_angles))

examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90

@@ -90,7 +90,6 @@ program rrtmgp_rfmip_sw
   character(len=132) :: kdist_file = 'coefficients_sw.nc'
   character(len=132) :: flxdn_file, flxup_file
   integer            :: nargs, ncol, nlay, nbnd, ngpt, nexp, nblocks, block_size, forcing_index
-  logical            :: top_at_1
   integer            :: b, icol, ibnd, igpt
   character(len=4)   :: block_size_char, forcing_index_char = '1'
 
@@ -166,10 +165,6 @@ program rrtmgp_rfmip_sw
   ! Allocation on assignment within reading routines
   !
   call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev)
-  !
-  ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain?
-  !
-  top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1)
 
   !
   ! Read the gas concentrations and surface properties
@@ -193,7 +188,10 @@ program rrtmgp_rfmip_sw
   !   is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep.
   !   This introduces an error but shows input sanitizing.
   !
-  if(top_at_1) then
+  !
+  ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain?
+  !
+  if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then
     p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min())
   else
     p_lev(:,nlay+1,:) &
@@ -296,7 +294,6 @@ program rrtmgp_rfmip_sw
     !    via ty_fluxes_broadband
     !
     call stop_on_err(rte_sw(optical_props,   &
-                            top_at_1,        &
                             mu0,             &
                             toa_flux,        &
                             sfc_alb_spec,    &

extensions/mo_compute_bc.F90

@@ -166,7 +166,6 @@ function compute_bc(k_dist,                      &
       ! Compute fluxes
       !
       error_msg = rte_lw(optical_props_1lay, &
-                         top_at_1,           &
                          lw_sources_1lay,    &
                          lower_bc, fluxes_1lev)
     else
@@ -191,8 +190,7 @@ function compute_bc(k_dist,                      &
                                     optical_props_1lay,         &
                                     solar_src)
       error_msg = rte_sw(optical_props_1lay, &
-                         top_at_1, mu0,      &
-                         solar_src,          &
+                         mu0, solar_src,     &
                          lower_bc, lower_bc, fluxes_1lev)
     endif
   end function

extensions/mo_rrtmgp_clr_all_sky.F90

@@ -85,11 +85,6 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev,    &
     ngpt  = k_dist%get_ngpt()
     nband = k_dist%get_nband()
 
-    !$acc kernels copyout(top_at_1)
-    !$omp target map(from:top_at_1)
-    top_at_1 = p_lay(1, 1) < p_lay(1, nlay)
-    !$acc end kernels
-    !$omp end target
 
     ! ------------------------------------------------------------------------------------
     !  Error checking
@@ -161,8 +156,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev,    &
     if(present(aer_props)) error_msg = aer_props%increment(optical_props)
     if(error_msg /= '') return
 
-    error_msg = base_rte_lw(optical_props, top_at_1, sources, &
-                            sfc_emis, clrsky_fluxes,          &
+    error_msg = base_rte_lw(optical_props, sources,  &
+                            sfc_emis, clrsky_fluxes, &
                             inc_flux, n_gauss_angles)
     if(error_msg /= '') return
     ! ------------------------------------------------------------------------------------
@@ -171,8 +166,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev,    &
     error_msg = cloud_props%increment(optical_props)
     if(error_msg /= '') return
 
-    error_msg = base_rte_lw(optical_props, top_at_1, sources, &
-                            sfc_emis, allsky_fluxes,          &
+    error_msg = base_rte_lw(optical_props, sources,  &
+                            sfc_emis, allsky_fluxes, &
                             inc_flux, n_gauss_angles)
 
     call sources%finalize()
@@ -207,7 +202,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     class(ty_optical_props_arry), allocatable :: optical_props
     real(wp), dimension(:,:),     allocatable :: toa_flux
     integer :: ncol, nlay, ngpt, nband, nstr
-    logical :: top_at_1
     ! --------------------------------
     ! Problem sizes
     !
@@ -218,12 +212,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     ngpt  = k_dist%get_ngpt()
     nband = k_dist%get_nband()
 
-    !$acc kernels copyout(top_at_1)
-    !$omp target map(from:top_at_1)
-    top_at_1 = p_lay(1, 1) < p_lay(1, nlay)
-    !$acc end kernels
-    !$omp end target
-
     ! ------------------------------------------------------------------------------------
     !  Error checking
     !
@@ -276,7 +264,7 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     ! Gas optical depth -- pressure need to be expressed as Pa
     !
     error_msg = k_dist%gas_optics(p_lay, p_lev, t_lay, gas_concs,  &
-                                  optical_props, toa_flux,                          &
+                                  optical_props, toa_flux,         &
                                   col_dry)
     if (error_msg /= '') return
     !
@@ -289,9 +277,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     if(present(aer_props)) error_msg = aer_props%increment(optical_props)
     if(error_msg /= '') return
 
-    error_msg = base_rte_sw(optical_props, top_at_1, &
-                               mu0, toa_flux,                   &
-                               sfc_alb_dir, sfc_alb_dif,        &
+    error_msg = base_rte_sw(optical_props, mu0, toa_flux, &
+                               sfc_alb_dir, sfc_alb_dif,  &
                                clrsky_fluxes)
 
     if(error_msg /= '') return
@@ -301,9 +288,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     error_msg = cloud_props%increment(optical_props)
     if(error_msg /= '') return
 
-    error_msg = base_rte_sw(optical_props, top_at_1,  &
-                               mu0, toa_flux,                   &
-                               sfc_alb_dir, sfc_alb_dif,        &
+    error_msg = base_rte_sw(optical_props, mu0, toa_flux, &
+                               sfc_alb_dir, sfc_alb_dif,  &
                                allsky_fluxes)
 
     call optical_props%finalize()

rrtmgp-frontend/mo_gas_optics_rrtmgp.F90

@@ -252,12 +252,17 @@ function gas_optics_int(this,                             &
     integer,     dimension(2,    size(play,dim=1),size(play,dim=2), get_nflav(this)) :: jeta
 
     integer :: ncol, nlay, ngpt, nband
+    logical :: temp
     ! ----------------------------------------------------------
     ncol  = size(play,dim=1)
     nlay  = size(play,dim=2)
     ngpt  = this%get_ngpt()
     nband = this%get_nband()
     !
+    ! Vertical orientation
+    !
+    call optical_props%set_top_at_1(play(1,1) < play(1, nlay))
+    !
     ! Gas optics
     !
     !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
@@ -311,7 +316,7 @@ function gas_optics_int(this,                             &
       !   but isn't with PGI 19.10
       !
       error_msg = source(this,                               &
-                         ncol, nlay, nband, ngpt,            &
+                         ncol, nlay, nband, ngpt, optical_props%top_is_at_1(),  &
                          play, plev, tlay, tsfc,             &
                          jtemp, jpress, jeta, tropo, fmajor, &
                          sources,                            &
@@ -320,7 +325,7 @@ function gas_optics_int(this,                             &
       !$omp target exit data map(release:tlev)
     else
       error_msg = source(this,                               &
-                         ncol, nlay, nband, ngpt,            &
+                         ncol, nlay, nband, ngpt, optical_props%top_is_at_1(),  &
                          play, plev, tlay, tsfc,             &
                          jtemp, jpress, jeta, tropo, fmajor, &
                          sources)
@@ -329,6 +334,7 @@ function gas_optics_int(this,                             &
     !$omp target exit data map(release:tsfc)
     !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta)
     !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
+
   end function gas_optics_int
   !------------------------------------------------------------------------------------------
   !
@@ -372,6 +378,10 @@ function gas_optics_ext(this,                         &
     ngas  = this%get_ngas()
     nflav = get_nflav(this)
     !
+    ! Vertical orientation
+    !
+    call optical_props%set_top_at_1(play(1,1) < play(1, nlay))
+    !
     ! Gas optics
     !
     !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
@@ -384,6 +394,7 @@ function gas_optics_ext(this,                         &
                                  col_dry)
     !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta)
     !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
+
     if(error_msg  /= '') return
 
     ! ----------------------------------------------------------
@@ -824,7 +835,7 @@ end function set_tsi
   ! Compute Planck source functions at layer centers and levels
   !
   function source(this,                               &
-                  ncol, nlay, nbnd, ngpt,             &
+                  ncol, nlay, nbnd, ngpt, top_at_1,   &
                   play, plev, tlay, tsfc,             &
                   jtemp, jpress, jeta, tropo, fmajor, &
                   sources,                            & ! Planck sources
@@ -833,6 +844,7 @@ function source(this,                               &
     ! inputs
     class(ty_gas_optics_rrtmgp),    intent(in ) :: this
     integer,                               intent(in   ) :: ncol, nlay, nbnd, ngpt
+    logical,                               intent(in   ) :: top_at_1
     real(wp), dimension(ncol,nlay),        intent(in   ) :: play   ! layer pressures [Pa, mb]
     real(wp), dimension(ncol,nlay+1),      intent(in   ) :: plev   ! level pressures [Pa, mb]
     real(wp), dimension(ncol,nlay),        intent(in   ) :: tlay   ! layer temperatures [K]
@@ -849,7 +861,6 @@ function source(this,                               &
                                       optional, target :: tlev          ! level temperatures [K]
     character(len=128)                                 :: error_msg
     ! ----------------------------------------------------------
-    logical(wl)                                  :: top_at_1
     integer                                      :: icol, ilay
     ! Variables for temperature at layer edges [K] (ncol, nlay+1)
     real(wp), dimension(   ncol,nlay+1), target  :: tlev_arr
@@ -895,18 +906,12 @@ function source(this,                               &
     !$omp target data                 map(from:sources%lay_source, sources%lev_source) &
     !$omp                             map(from:sources%sfc_source, sources%sfc_source_Jac)
 
-    !$acc kernels copyout(top_at_1)
-    !$omp target map(from:top_at_1)
-    top_at_1 = play(1,1) < play(1, nlay)
-    !$acc end kernels
-    !$omp end target
-
     call compute_Planck_source(ncol, nlay, nbnd, ngpt, &
                 get_nflav(this), this%get_neta(), this%get_npres(), this%get_ntemp(), this%get_nPlanckTemp(), &
-                tlay, tlev_wk, tsfc, merge(nlay, 1, top_at_1), &
-                fmajor, jeta, tropo, jtemp, jpress,                    &
+                tlay, tlev_wk, tsfc, merge(nlay, 1, logical(top_at_1, wl)), &
+                fmajor, jeta, tropo, jtemp, jpress,                         &
                 this%get_gpoint_bands(), this%get_band_lims_gpoint(), this%planck_frac, this%temp_ref_min,&
-                this%totplnk_delta, this%totplnk, this%gpoint_flavor,  &
+                this%totplnk_delta, this%totplnk, this%gpoint_flavor,       &
                 sources%sfc_source, sources%lay_source, sources%lev_source, &
                 sources%sfc_source_Jac)
     !$acc end        data