Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 4573)
+++ /palm/trunk/SOURCE/modules.f90	(revision 4574)
@@ -20,9 +20,14 @@
 ! Current revisions:
 ! ------------------
-! Revise description of rotation_angle
-!
+! 
+! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+!
+! 4572 2020-06-24 10:09:41Z gronemeier
+! Revise description of rotation_angle
+!
+! 4565 2020-06-15 08:30:38Z oliver.maas
 ! added pt_surface_heating_rate
 !
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4573)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4574)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! - Restructure code in radiation_check_data_output
+! - Move calculation of MPI global array offsets to a subroutine
+! 
+! 4571 2020-06-24 08:59:06Z sebschub
 ! Bugfix in vertical lad_s coordinate
 ! 
@@ -1138,53 +1142,65 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=*) ::  unit          !<
-       CHARACTER (LEN=*) ::  variable      !<
-
-       INTEGER(iwp) :: i, k
-       INTEGER(iwp) :: ilen
-       CHARACTER(LEN=varnamelength) :: var  !< TRIM(variable)
+       LOGICAL                      ::  directional
+       CHARACTER(LEN=*)             ::  unit        !<
+       CHARACTER(LEN=*)             ::  variable    !<
+       CHARACTER(LEN=varnamelength) ::  var         !< TRIM(variable)
+       INTEGER(iwp)                 ::  i, k
+       INTEGER(iwp)                 ::  ilast_word
+       INTEGER(iwp)                 ::  ilen
 
        var = TRIM(variable)
-
-       IF ( len(var) < 3_iwp  )  THEN
-          unit = 'illegal'
-          RETURN
-       ENDIF
-
-       IF ( var(1:3) /= 'rad'  .AND.  var(1:3) /= 'rtm' )  THEN
-          unit = 'illegal'
-          RETURN
-       ENDIF
-
-!--    first process diractional variables
-       IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
-            var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
-            var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
-            var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
-            var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
-            var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
-          IF ( .NOT.  radiation ) THEN
-                message_string = 'output of "' // TRIM( var ) // '" require'&
-                                 // 's radiation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
-          ENDIF
-          unit = 'W/m2'
-       ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
-                 var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft'  .OR.             &
-                 var(1:12) == 'rtm_surfalb_'  .OR.  var(1:13) == 'rtm_surfemis_'  ) THEN
-          IF ( .NOT.  radiation ) THEN
-                message_string = 'output of "' // TRIM( var ) // '" require'&
-                                 // 's radiation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
-          ENDIF
-          unit = '1'
+!
+!--    Identify directional variables
+       ilast_word = SCAN(var, '_', back=.TRUE.)
+       IF ( ilast_word > 0 )  THEN
+          SELECT CASE ( var(ilast_word:) )
+             CASE ( '_roof', '_south', '_north', '_west', '_east' )
+                directional = .TRUE.
+                write(9, *) 'vardir', var
+                flush(9)
+                var = var(1:ilast_word-1)
+             CASE DEFAULT
+                directional = .FALSE.
+                write(9, *) 'varnd', var
+                flush(9)
+          END SELECT
        ELSE
-!--       non-directional variables
-          SELECT CASE ( TRIM( var ) )
-             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
+          directional = .FALSE.
+       END IF
+
+       IF ( directional )  THEN
+          SELECT CASE ( var )
+             CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw',             &
+                    'rtm_rad_inswdir', 'rtm_rad_inswdif', 'rtm_rad_inswref',   &
+                    'rtm_rad_inlwdif', 'rtm_rad_inlwref', 'rtm_rad_outsw',     &
+                    'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'        &
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W/m2'
+
+             CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft',           &
+                    'rtm_surfalb', 'rtm_surfemis' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'&
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = '1'
+
+             CASE DEFAULT
+                unit = 'illegal'
+          END SELECT
+
+       ELSE
+          SELECT CASE ( var )
+             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out',     &
                     'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
                 IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
-                   message_string = '"output of "' // TRIM( var ) // '" requi' // &
-                                    'res radiation = .TRUE. and ' //              &
+                   message_string = '"output of "' // var // '" requi' //       &
+                                    'res radiation = .TRUE. and ' //            &
                                     'radiation_scheme = "rrtmg"'
                    CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
@@ -1192,46 +1208,33 @@
                 unit = 'K/h'
 
-             CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
-                    'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
-                    'rad_sw_out*')
-                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
-                   ! Workaround for masked output (calls with i=ilen=k=0)
-                   unit = 'illegal'
-                   RETURN
+             CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',      &
+                    'rad_sw_out*' )
+                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                   message_string = 'illegal value for data_output: "' //       &
+                                    var // '" & only 2d-horizontal ' //         &
+                                    'cross sections are allowed for this value'
+                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
                 ENDIF
+                unit = 'W/m2'
+
+             CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
                 IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-                   message_string = 'illegal value for data_output: "' //         &
-                                    TRIM( var ) // '" & only 2d-horizontal ' //   &
+                   message_string = 'illegal value for data_output: "' //       &
+                                    var // '" & only 2d-horizontal ' //         &
                                     'cross sections are allowed for this value'
                    CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
                 ENDIF
                 IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
-                   IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
-                        TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
-                        TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
-                        TRIM( var ) == 'rrtm_asdir*'      )                       &
-                   THEN
-                      message_string = 'output of "' // TRIM( var ) // '" require'&
-                                       // 's radiation = .TRUE. and radiation_sch'&
-                                       // 'eme = "rrtmg"'
-                      CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
-                   ENDIF
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's radiation = .TRUE. and radiation_sch' &
+                                    // 'eme = "rrtmg"'
+                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
                 ENDIF
-
-                IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
-                IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''
-                IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
-                IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
-                IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
+                unit = ''
 
              CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
                     'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
                 IF ( .NOT.  radiation ) THEN
-                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                   message_string = 'output of "' // var // '" require'         &
                                     // 's radiation = .TRUE.'
                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
@@ -1240,26 +1243,20 @@
 
              CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
-                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
-                   ! Workaround for masked output (calls with i=ilen=k=0)
-                   unit = 'illegal'
-                   RETURN
-                ENDIF
-
                 IF ( .NOT.  radiation ) THEN
-                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                   message_string = 'output of "' // var // '" require'         &
                                     // 's radiation = .TRUE.'
                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
                 ENDIF
                 IF ( mrt_nlevels == 0 ) THEN
-                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                   message_string = 'output of "' // var // '" require'         &
                                     // 's mrt_nlevels > 0'
                    CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
                 ENDIF
-                IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
-                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
+                   message_string = 'output of "' // var // '" require'         &
                                     // 's rtm_mrt_sw = .TRUE.'
                    CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
                 ENDIF
-                IF ( TRIM( var ) == 'rtm_mrt' ) THEN
+                IF ( var == 'rtm_mrt' ) THEN
                    unit = 'K'
                 ELSE
@@ -1271,5 +1268,5 @@
 
           END SELECT
-       ENDIF
+       END IF
 
     END SUBROUTINE radiation_check_data_output
@@ -8516,7 +8513,4 @@
         INTEGER(iwp), DIMENSION(3)             :: dimnext      !< next dimension increments along path
         INTEGER(iwp), DIMENSION(3)             :: dimdelta     !< dimension direction = +- 1
-        INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before
-                                                               !< the processor in the question
-        INTEGER(iwp)                           :: ip           !< number of processor where gridbox reside
         INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
 
@@ -8583,8 +8577,5 @@
 !--             calculate index of the grid with global indices (box(2),box(3))
 !--             in the array nzterr and plantt and id of the coresponding processor
-                px = box(3)/nnx
-                py = box(2)/nny
-                ip = px*pdims(2)+py
-                ig = ip*nnx*nny + (box(3)-px*nnx)*nny + box(2)-py*nny
+                CALL radiation_calc_global_offset( box(3), box(2), 0, 1, offs_glob=ig )
                 IF ( box(1) <= nzterr(ig) )  THEN
                     visible = .FALSE.
@@ -8598,6 +8589,7 @@
                         crlens(ncsb) = crlen
 #if defined( __parallel )
-                        lad_ip(ncsb) = ip
-                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz_plant) + (box(2)-py*nny)*nz_plant + box(1)-nz_urban_b
+                        CALL radiation_calc_global_offset( box(3), box(2), box(1)-nz_urban_b, &
+                                                           nz_plant, iproc=lad_ip(ncsb),      &
+                                                           offs_proc=lad_disp(ncsb) )
 #endif
                     ENDIF
@@ -8716,6 +8708,4 @@
       INTEGER(iwp), DIMENSION(2)             ::  dimnext      !< next dimension increments along path
       INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
-      INTEGER(iwp)                           ::  px, py       !< number of processors in x and y dir before
-                                                              !< the processor in the question
       INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
       INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
@@ -8821,8 +8811,5 @@
 !--         calculate index of the grid with global indices (column(1),column(2))
 !--         in the array nzterr and plantt and id of the coresponding processor
-            px = column(2)/nnx
-            py = column(1)/nny
-            ip = px*pdims(2)+py
-            ig = ip*nnx*nny + (column(2)-px*nnx)*nny + column(1)-py*nny
+            CALL radiation_calc_global_offset( column(2), column(1), 0, 1, offs_glob=ig )
 
             IF ( lastdist == 0._wp )  THEN
@@ -8898,8 +8885,6 @@
             !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
             DO  i = 1, ntrack
-               px = rt2_track(2,i)/nnx
-               py = rt2_track(1,i)/nny
-               ip = px*pdims(2)+py
-               ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
+               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, &
+                                                  offs_glob=ig )
 
                IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
@@ -8919,5 +8904,7 @@
                IF ( plantt(ig) < lowest_lad )  CYCLE
 
-               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant + lowest_lad-nz_urban_b
+               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i),                  &
+                                                  lowest_lad-nz_urban_b, nz_plant, iproc=ip,       &
+                                                  offs_proc=wdisp )
                wcount = plantt(ig)-lowest_lad+1
                ! TODO send request ASAP - even during raytracing
@@ -8942,8 +8929,6 @@
          ELSE ! raytrace_mpi_rma = .F.
             DO  i = 1, ntrack
-               px = rt2_track(2,i)/nnx
-               py = rt2_track(1,i)/nny
-               ip = px*pdims(2)+py
-               ig = ip*nnx*nny*nz_plant + (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant
+               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, nz_plant, &
+                                                  offs_glob=ig )
                rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
             ENDDO
@@ -9011,7 +8996,5 @@
             DO  i = 1, ntrack ! for each column
                dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
-               px = rt2_track(2,i)/nnx
-               py = rt2_track(1,i)/nny
-               ip = px*pdims(2)+py
+               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, iproc=ip )
 
                DO  k = 1, nrays ! for each ray
@@ -9081,7 +9064,5 @@
          DO  i = ntrack, 1, -1 ! for each column backwards
             dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
-            px = rt2_track(2,i)/nnx
-            py = rt2_track(1,i)/nny
-            ip = px*pdims(2)+py
+            CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, iproc=ip )
 
             DO  k = 1, nrays ! for each ray
@@ -9150,8 +9131,4 @@
          INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
          INTEGER(iwp), INTENT(out)           ::  iproc
-#if defined( __parallel )
-         INTEGER(iwp)                        ::  px, py        !< number of processors in x and y direction
-                                                               !< before the processor in the question
-#endif
 
 #if defined( __parallel )
@@ -9160,9 +9137,7 @@
 !
 !--      Calculate target processor and index in the remote local target gridsurf array
-         px = x / nnx
-         py = y / nny
-         iproc = px * pdims(2) + py
-         target_displ = ((x-px*nnx) * nny + y - py*nny ) * nz_urban * nsurf_type_u +&
-                        ( z-nz_urban_b ) * nsurf_type_u + d
+         CALL radiation_calc_global_offset( x, y, (z - nz_urban_b) * nsurf_type_u + d, &
+                                            nz_urban * nsurf_type_u, iproc=iproc,      &
+                                            offs_proc=target_displ )
 !
 !--      Send MPI_Get request to obtain index target_surfl(i)
@@ -9854,4 +9829,50 @@
 ! Description:
 ! ------------
+!> For given coordinates, calculates indices within a global 3D (or 2D if
+!> nlayers=1) field, e.g. an MPI one-sided window or an array which has been
+!> created using e.g. MPI_AllGather.
+!------------------------------------------------------------------------------!
+ PURE SUBROUTINE radiation_calc_global_offset( i, j, k, nlayers,           &
+                                               iproc, offs_proc, offs_glob )
+
+    IMPLICIT NONE
+
+    INTEGER(iwp), INTENT(IN)                              ::  i           !< x-coordinate
+    INTEGER(iwp), INTENT(IN)                              ::  j           !< y-coordinate
+    INTEGER(iwp), INTENT(IN)                              ::  k           !< z-coordinate
+    INTEGER(iwp), INTENT(IN)                              ::  nlayers     !< number of z-layers
+    INTEGER(iwp), INTENT(OUT), OPTIONAL                   ::  iproc       !< MPI process rank
+    INTEGER(kind=MPI_ADDRESS_KIND), INTENT(OUT), OPTIONAL ::  offs_proc   !< offset within MPI proc
+    INTEGER(iwp), INTENT(OUT), OPTIONAL                   ::  offs_glob   !< global offset
+
+    INTEGER(iwp)            ::  ipx         !< process index in x-direction
+    INTEGER(iwp)            ::  ipy         !< process index in y-direction
+    INTEGER(iwp)            ::  iproc_l     !< local variable for iproc
+    INTEGER(iwp)            ::  oproc_l     !< local variable for offs_proc
+    INTEGER(iwp)            ::  offs_pstart !< global start of the MPI process
+
+    ipx = i / nnx
+    ipy = j / nny
+    iproc_l = ipx * pdims(2) + ipy
+    IF ( PRESENT(iproc) )  iproc = iproc_l
+
+    IF ( PRESENT(offs_proc)  .OR.  PRESENT(offs_glob) )  THEN
+       oproc_l = (i - ipx*nnx) * nny * nlayers + & ! columns before
+                 (j - ipy*nny) * nlayers       + & ! rows in column
+                 k
+       IF ( PRESENT(offs_proc) )  offs_proc = oproc_l
+       IF ( PRESENT(offs_glob) )  THEN
+          offs_pstart = iproc_l * nnx * nny * nlayers
+          offs_glob = offs_pstart + oproc_l
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE radiation_calc_global_offset
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
 !> Subroutine for averaging 3D data
 !------------------------------------------------------------------------------!
