Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4694)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4708)
@@ -28,4 +28,9 @@
 ! -----------------
 ! $Id$
+! - Bugfix, correct mapping of RRTMG heating rates, as well as incoming/outgoing
+!   radiation onto the topography-following grid
+! - add fill values to the output
+!
+! 4694 2020-09-23 15:09:19Z pavelkrc
 ! Bugfix for tracing of maximum radiative fluxes
 !
@@ -401,5 +406,5 @@
     USE indices,                                                               &
         ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
-               nzb, nzt, topo_top_ind
+               nzb, nzt, topo_top_ind, wall_flags_total_0
 
     USE, INTRINSIC :: iso_c_binding
@@ -4390,9 +4395,14 @@
 !
 !--          Save fluxes
-             DO k = nzb, nzt+1
-                rad_lw_in(k,:,:)  = rrtm_lwdflx(0,k)
-                rad_lw_out(k,:,:) = rrtm_lwuflx(0,k)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   k_topo_l = topo_top_ind(j,i,0)
+                   DO k = k_topo_l, nzt+1
+                      rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k-k_topo_l+k_topo)
+                      rad_lw_out(k,j,i) = rrtm_lwuflx(0,k-k_topo_l+k_topo)
+                   ENDDO
+                ENDDO
              ENDDO
-             rad_lw_in_diff(:,:) = rad_lw_in(k_topo,:,:)
+             rad_lw_in_diff(:,:) = rrtm_lwdflx(0,k_topo)
 !
 !--          Save heating rates (convert from K/d to K/h).
@@ -4404,6 +4414,6 @@
                    k_topo_l = topo_top_ind(j,i,0)
                    DO k = k_topo_l+1, nzt+1
-                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo_l)  * d_hours_day
-                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo_l) * d_hours_day
+                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo_l+k_topo)  * d_hours_day
+                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo_l+k_topo) * d_hours_day
                    ENDDO
                 ENDDO
@@ -4499,7 +4509,12 @@
 !
 !--          Save radiation fluxes for the entire depth of the model domain
-             DO k = nzb, nzt+1
-                rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
-                rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   k_topo_l = topo_top_ind(j,i,0)
+                   DO k = k_topo_l, nzt+1
+                      rad_sw_in(k,j,i)  = rrtm_swdflx(0,k-k_topo_l+k_topo)
+                      rad_sw_out(k,j,i) = rrtm_swuflx(0,k-k_topo_l+k_topo)
+                   ENDDO
+                ENDDO
              ENDDO
 !--          Save direct and diffuse SW radiation at the surface (required by RTM)
@@ -4509,7 +4524,12 @@
 !
 !--          Save heating rates (convert from K/d to K/s)
-             DO k = nzb+1, nzt+1
-                rad_sw_hr(k,:,:)     = rrtm_swhr(0,k)  * d_hours_day
-                rad_sw_cs_hr(k,:,:)  = rrtm_swhrc(0,k) * d_hours_day
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   k_topo_l = topo_top_ind(j,i,0)
+                   DO k = k_topo_l+1, nzt+1
+                      rad_sw_hr(k,j,i)     = rrtm_swhr(0,k-k_topo_l+k_topo)  * d_hours_day
+                      rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k-k_topo_l+k_topo) * d_hours_day
+                   ENDDO
+                ENDDO
              ENDDO
 !
@@ -4523,5 +4543,7 @@
 !
 !--    RRTMG is called for each (j,i) grid point separately, starting at the
-!--    highest topography level. Here no RTM is used since average_radiation is false
+!--    highest topography level. Here no RTM is used since average_radiation is false.
+!--    In fact, this branch is only called for homogeneous flat terrain so that the topography-top
+!--    index is actually alwasy zero.
        ELSE
 !
@@ -4694,5 +4716,5 @@
                 ENDDO
 !
-!--             Obtain topography top index (lower bound of RRTMG)
+!--             Obtain topography top index (lower bound of RRTMG). Is actually always zero.
                 k_topo = topo_top_ind(j,i,0)
 
@@ -5829,6 +5851,8 @@
 !--    Calculate tendency based on heating rate
        DO k = nzb+1, nzt+1
-          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
-                                         * d_exner(k) * d_seconds_hour
+          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))                        &
+                                         * d_exner(k) * d_seconds_hour                             &
+                                         * MERGE( 1.0_wp, 0.0_wp,                                  &
+                                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
        ENDDO
 
@@ -5848,7 +5872,4 @@
  SUBROUTINE radiation_tendency ( tend )
 
-    USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys
-
     IMPLICIT NONE
 
@@ -5863,7 +5884,9 @@
           DO  j = nys, nyn
              DO k = nzb+1, nzt+1
-                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
-                                          +  rad_sw_hr(k,j,i) ) * d_exner(k)   &
-                                          * d_seconds_hour
+                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                                      &
+                                          + rad_sw_hr(k,j,i) ) * d_exner(k)                         &
+                                          * d_seconds_hour                                          &
+                                          * MERGE( 1.0_wp, 0.0_wp,                                  &
+                                                   BTEST( wall_flags_total_0(k,j,i), 0 ) )
              ENDDO
           ENDDO
@@ -11137,5 +11160,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
-                                      local_pf, two_d, nzb_do, nzt_do )
+                                      local_pf, two_d, nzb_do, nzt_do, fill_value )
 
     USE indices
@@ -11161,5 +11184,5 @@
     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
 
-    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
+    REAL(wp) ::  fill_value !< value for the _FillValue attribute
 
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
@@ -11340,5 +11363,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_lw_in(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_lw_in(k,j,i), fill_value,                       &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11364,5 +11388,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_lw_out(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_lw_out(k,j,i), fill_value,                      &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11388,5 +11413,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_lw_cs_hr(k,j,i), fill_value,                    &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11412,5 +11438,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_lw_hr(k,j,i), fill_value,                       &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11436,5 +11463,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_sw_in(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_sw_in(k,j,i), fill_value,                       &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11460,5 +11488,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_sw_out(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_sw_out(k,j,i), fill_value,                      &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11484,5 +11513,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_sw_cs_hr(k,j,i), fill_value,                    &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11508,5 +11538,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                      local_pf(i,j,k) = MERGE( rad_sw_hr(k,j,i), fill_value,                       &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -11542,5 +11573,5 @@
 !> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
+ SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
 
 
@@ -11561,5 +11592,5 @@
     LOGICAL      ::  found       !<
 
-    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
+    REAL(wp)     ::  fill_value  !< value for the _FillValue attribute
 
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
@@ -11626,5 +11657,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_sw_in(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_sw_in(k,j,i), fill_value,                        &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11649,5 +11681,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_sw_out(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_sw_out(k,j,i), fill_value,                       &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11672,5 +11705,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_sw_cs_hr(k,j,i), fill_value,                     &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11695,5 +11729,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_sw_hr(k,j,i), fill_value,                        &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11718,5 +11753,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_lw_in(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_lw_in(k,j,i), fill_value,                        &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11741,5 +11777,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_lw_out(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_lw_out(k,j,i), fill_value,                       &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11764,5 +11801,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_lw_cs_hr(k,j,i), fill_value,                     &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -11787,5 +11825,6 @@
                DO  j = nys, nyn
                   DO  k = nzb_do, nzt_do
-                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                     local_pf(i,j,k) = MERGE( rad_lw_hr(k,j,i), fill_value,                        &
+                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
                   ENDDO
                ENDDO
@@ -12227,5 +12266,6 @@
                 DO  k = 1, mask_size_l(mid,3)
                    local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
-                                      mask_j(mid,j),mask_i(mid,i))
+                                                     mask_j(mid,j), &
+                                                     mask_i(mid,i))
                 ENDDO
              ENDDO
