Ignore:
Timestamp:
Jun 24, 2020 4:33:32 PM (4 years ago)
Author:
pavelkrc
Message:

Move calculation of MPI global array offsets to a subroutine, restructure code in radiation_check_data_output

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4571 r4574  
    2828! -----------------
    2929! $Id$
     30! - Restructure code in radiation_check_data_output
     31! - Move calculation of MPI global array offsets to a subroutine
     32!
     33! 4571 2020-06-24 08:59:06Z sebschub
    3034! Bugfix in vertical lad_s coordinate
    3135!
     
    11381142       IMPLICIT NONE
    11391143
    1140        CHARACTER (LEN=*) ::  unit          !<
    1141        CHARACTER (LEN=*) ::  variable      !<
    1142 
    1143        INTEGER(iwp) :: i, k
    1144        INTEGER(iwp) :: ilen
    1145        CHARACTER(LEN=varnamelength) :: var  !< TRIM(variable)
     1144       LOGICAL                      ::  directional
     1145       CHARACTER(LEN=*)             ::  unit        !<
     1146       CHARACTER(LEN=*)             ::  variable    !<
     1147       CHARACTER(LEN=varnamelength) ::  var         !< TRIM(variable)
     1148       INTEGER(iwp)                 ::  i, k
     1149       INTEGER(iwp)                 ::  ilast_word
     1150       INTEGER(iwp)                 ::  ilen
    11461151
    11471152       var = TRIM(variable)
    1148 
    1149        IF ( len(var) < 3_iwp  )  THEN
    1150           unit = 'illegal'
    1151           RETURN
    1152        ENDIF
    1153 
    1154        IF ( var(1:3) /= 'rad'  .AND.  var(1:3) /= 'rtm' )  THEN
    1155           unit = 'illegal'
    1156           RETURN
    1157        ENDIF
    1158 
    1159 !--    first process diractional variables
    1160        IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
    1161             var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
    1162             var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
    1163             var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
    1164             var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
    1165             var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
    1166           IF ( .NOT.  radiation ) THEN
    1167                 message_string = 'output of "' // TRIM( var ) // '" require'&
    1168                                  // 's radiation = .TRUE.'
    1169                 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1170           ENDIF
    1171           unit = 'W/m2'
    1172        ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
    1173                  var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft'  .OR.             &
    1174                  var(1:12) == 'rtm_surfalb_'  .OR.  var(1:13) == 'rtm_surfemis_'  ) THEN
    1175           IF ( .NOT.  radiation ) THEN
    1176                 message_string = 'output of "' // TRIM( var ) // '" require'&
    1177                                  // 's radiation = .TRUE.'
    1178                 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1179           ENDIF
    1180           unit = '1'
     1153!
     1154!--    Identify directional variables
     1155       ilast_word = SCAN(var, '_', back=.TRUE.)
     1156       IF ( ilast_word > 0 )  THEN
     1157          SELECT CASE ( var(ilast_word:) )
     1158             CASE ( '_roof', '_south', '_north', '_west', '_east' )
     1159                directional = .TRUE.
     1160                write(9, *) 'vardir', var
     1161                flush(9)
     1162                var = var(1:ilast_word-1)
     1163             CASE DEFAULT
     1164                directional = .FALSE.
     1165                write(9, *) 'varnd', var
     1166                flush(9)
     1167          END SELECT
    11811168       ELSE
    1182 !--       non-directional variables
    1183           SELECT CASE ( TRIM( var ) )
    1184              CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
     1169          directional = .FALSE.
     1170       END IF
     1171
     1172       IF ( directional )  THEN
     1173          SELECT CASE ( var )
     1174             CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw',             &
     1175                    'rtm_rad_inswdir', 'rtm_rad_inswdif', 'rtm_rad_inswref',   &
     1176                    'rtm_rad_inlwdif', 'rtm_rad_inlwref', 'rtm_rad_outsw',     &
     1177                    'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
     1178                IF ( .NOT.  radiation ) THEN
     1179                   message_string = 'output of "' // var // '" require'        &
     1180                                    // 's radiation = .TRUE.'
     1181                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1182                ENDIF
     1183                unit = 'W/m2'
     1184
     1185             CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft',           &
     1186                    'rtm_surfalb', 'rtm_surfemis' )
     1187                IF ( .NOT.  radiation ) THEN
     1188                   message_string = 'output of "' // var // '" require'&
     1189                                    // 's radiation = .TRUE.'
     1190                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1191                ENDIF
     1192                unit = '1'
     1193
     1194             CASE DEFAULT
     1195                unit = 'illegal'
     1196          END SELECT
     1197
     1198       ELSE
     1199          SELECT CASE ( var )
     1200             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out',     &
    11851201                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
    11861202                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    1187                    message_string = '"output of "' // TRIM( var ) // '" requi' // &
    1188                                     'res radiation = .TRUE. and ' //              &
     1203                   message_string = '"output of "' // var // '" requi' //      &
     1204                                    'res radiation = .TRUE. and ' //            &
    11891205                                    'radiation_scheme = "rrtmg"'
    11901206                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
     
    11921208                unit = 'K/h'
    11931209
    1194              CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
    1195                     'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
    1196                     'rad_sw_out*')
    1197                 IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
    1198                    ! Workaround for masked output (calls with i=ilen=k=0)
    1199                    unit = 'illegal'
    1200                    RETURN
     1210             CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',      &
     1211                    'rad_sw_out*' )
     1212                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     1213                   message_string = 'illegal value for data_output: "' //       &
     1214                                    var // '" & only 2d-horizontal ' //         &
     1215                                    'cross sections are allowed for this value'
     1216                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    12011217                ENDIF
     1218                unit = 'W/m2'
     1219
     1220             CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
    12021221                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    1203                    message_string = 'illegal value for data_output: "' //         &
    1204                                     TRIM( var ) // '" & only 2d-horizontal ' //   &
     1222                   message_string = 'illegal value for data_output: "' //       &
     1223                                    var // '" & only 2d-horizontal ' //         &
    12051224                                    'cross sections are allowed for this value'
    12061225                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    12071226                ENDIF
    12081227                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
    1209                    IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
    1210                         TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
    1211                         TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
    1212                         TRIM( var ) == 'rrtm_asdir*'      )                       &
    1213                    THEN
    1214                       message_string = 'output of "' // TRIM( var ) // '" require'&
    1215                                        // 's radiation = .TRUE. and radiation_sch'&
    1216                                        // 'eme = "rrtmg"'
    1217                       CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
    1218                    ENDIF
     1228                   message_string = 'output of "' // var // '" require'         &
     1229                                    // 's radiation = .TRUE. and radiation_sch' &
     1230                                    // 'eme = "rrtmg"'
     1231                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
    12191232                ENDIF
    1220 
    1221                 IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
    1222                 IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
    1223                 IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
    1224                 IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
    1225                 IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
    1226                 IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
    1227                 IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''
    1228                 IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
    1229                 IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
    1230                 IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
     1233                unit = ''
    12311234
    12321235             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
    12331236                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
    12341237                IF ( .NOT.  radiation ) THEN
    1235                    message_string = 'output of "' // TRIM( var ) // '" require'&
     1238                   message_string = 'output of "' // var // '" require'         &
    12361239                                    // 's radiation = .TRUE.'
    12371240                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     
    12401243
    12411244             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
    1242                 IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
    1243                    ! Workaround for masked output (calls with i=ilen=k=0)
    1244                    unit = 'illegal'
    1245                    RETURN
    1246                 ENDIF
    1247 
    12481245                IF ( .NOT.  radiation ) THEN
    1249                    message_string = 'output of "' // TRIM( var ) // '" require'&
     1246                   message_string = 'output of "' // var // '" require'         &
    12501247                                    // 's radiation = .TRUE.'
    12511248                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    12521249                ENDIF
    12531250                IF ( mrt_nlevels == 0 ) THEN
    1254                    message_string = 'output of "' // TRIM( var ) // '" require'&
     1251                   message_string = 'output of "' // var // '" require'         &
    12551252                                    // 's mrt_nlevels > 0'
    12561253                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
    12571254                ENDIF
    1258                 IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
    1259                    message_string = 'output of "' // TRIM( var ) // '" require'&
     1255                IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
     1256                   message_string = 'output of "' // var // '" require'         &
    12601257                                    // 's rtm_mrt_sw = .TRUE.'
    12611258                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
    12621259                ENDIF
    1263                 IF ( TRIM( var ) == 'rtm_mrt' ) THEN
     1260                IF ( var == 'rtm_mrt' ) THEN
    12641261                   unit = 'K'
    12651262                ELSE
     
    12711268
    12721269          END SELECT
    1273        ENDIF
     1270       END IF
    12741271
    12751272    END SUBROUTINE radiation_check_data_output
     
    85168513        INTEGER(iwp), DIMENSION(3)             :: dimnext      !< next dimension increments along path
    85178514        INTEGER(iwp), DIMENSION(3)             :: dimdelta     !< dimension direction = +- 1
    8518         INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before
    8519                                                                !< the processor in the question
    8520         INTEGER(iwp)                           :: ip           !< number of processor where gridbox reside
    85218515        INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
    85228516
     
    85838577!--             calculate index of the grid with global indices (box(2),box(3))
    85848578!--             in the array nzterr and plantt and id of the coresponding processor
    8585                 px = box(3)/nnx
    8586                 py = box(2)/nny
    8587                 ip = px*pdims(2)+py
    8588                 ig = ip*nnx*nny + (box(3)-px*nnx)*nny + box(2)-py*nny
     8579                CALL radiation_calc_global_offset( box(3), box(2), 0, 1, offs_glob=ig )
    85898580                IF ( box(1) <= nzterr(ig) )  THEN
    85908581                    visible = .FALSE.
     
    85988589                        crlens(ncsb) = crlen
    85998590#if defined( __parallel )
    8600                         lad_ip(ncsb) = ip
    8601                         lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz_plant) + (box(2)-py*nny)*nz_plant + box(1)-nz_urban_b
     8591                        CALL radiation_calc_global_offset( box(3), box(2), box(1)-nz_urban_b, &
     8592                                                           nz_plant, iproc=lad_ip(ncsb),      &
     8593                                                           offs_proc=lad_disp(ncsb) )
    86028594#endif
    86038595                    ENDIF
     
    87168708      INTEGER(iwp), DIMENSION(2)             ::  dimnext      !< next dimension increments along path
    87178709      INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
    8718       INTEGER(iwp)                           ::  px, py       !< number of processors in x and y dir before
    8719                                                               !< the processor in the question
    87208710      INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
    87218711      INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
     
    88218811!--         calculate index of the grid with global indices (column(1),column(2))
    88228812!--         in the array nzterr and plantt and id of the coresponding processor
    8823             px = column(2)/nnx
    8824             py = column(1)/nny
    8825             ip = px*pdims(2)+py
    8826             ig = ip*nnx*nny + (column(2)-px*nnx)*nny + column(1)-py*nny
     8813            CALL radiation_calc_global_offset( column(2), column(1), 0, 1, offs_glob=ig )
    88278814
    88288815            IF ( lastdist == 0._wp )  THEN
     
    88988885            !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
    88998886            DO  i = 1, ntrack
    8900                px = rt2_track(2,i)/nnx
    8901                py = rt2_track(1,i)/nny
    8902                ip = px*pdims(2)+py
    8903                ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
     8887               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, &
     8888                                                  offs_glob=ig )
    89048889
    89058890               IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     
    89198904               IF ( plantt(ig) < lowest_lad )  CYCLE
    89208905
    8921                wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant + lowest_lad-nz_urban_b
     8906               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i),                  &
     8907                                                  lowest_lad-nz_urban_b, nz_plant, iproc=ip,       &
     8908                                                  offs_proc=wdisp )
    89228909               wcount = plantt(ig)-lowest_lad+1
    89238910               ! TODO send request ASAP - even during raytracing
     
    89428929         ELSE ! raytrace_mpi_rma = .F.
    89438930            DO  i = 1, ntrack
    8944                px = rt2_track(2,i)/nnx
    8945                py = rt2_track(1,i)/nny
    8946                ip = px*pdims(2)+py
    8947                ig = ip*nnx*nny*nz_plant + (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant
     8931               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, nz_plant, &
     8932                                                  offs_glob=ig )
    89488933               rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
    89498934            ENDDO
     
    90118996            DO  i = 1, ntrack ! for each column
    90128997               dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
    9013                px = rt2_track(2,i)/nnx
    9014                py = rt2_track(1,i)/nny
    9015                ip = px*pdims(2)+py
     8998               CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, iproc=ip )
    90168999
    90179000               DO  k = 1, nrays ! for each ray
     
    90819064         DO  i = ntrack, 1, -1 ! for each column backwards
    90829065            dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
    9083             px = rt2_track(2,i)/nnx
    9084             py = rt2_track(1,i)/nny
    9085             ip = px*pdims(2)+py
     9066            CALL radiation_calc_global_offset( rt2_track(2,i), rt2_track(1,i), 0, 1, iproc=ip )
    90869067
    90879068            DO  k = 1, nrays ! for each ray
     
    91509131         INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
    91519132         INTEGER(iwp), INTENT(out)           ::  iproc
    9152 #if defined( __parallel )
    9153          INTEGER(iwp)                        ::  px, py        !< number of processors in x and y direction
    9154                                                                !< before the processor in the question
    9155 #endif
    91569133
    91579134#if defined( __parallel )
     
    91609137!
    91619138!--      Calculate target processor and index in the remote local target gridsurf array
    9162          px = x / nnx
    9163          py = y / nny
    9164          iproc = px * pdims(2) + py
    9165          target_displ = ((x-px*nnx) * nny + y - py*nny ) * nz_urban * nsurf_type_u +&
    9166                         ( z-nz_urban_b ) * nsurf_type_u + d
     9139         CALL radiation_calc_global_offset( x, y, (z - nz_urban_b) * nsurf_type_u + d, &
     9140                                            nz_urban * nsurf_type_u, iproc=iproc,      &
     9141                                            offs_proc=target_displ )
    91679142!
    91689143!--      Send MPI_Get request to obtain index target_surfl(i)
     
    98549829! Description:
    98559830! ------------
     9831!> For given coordinates, calculates indices within a global 3D (or 2D if
     9832!> nlayers=1) field, e.g. an MPI one-sided window or an array which has been
     9833!> created using e.g. MPI_AllGather.
     9834!------------------------------------------------------------------------------!
     9835 PURE SUBROUTINE radiation_calc_global_offset( i, j, k, nlayers,           &
     9836                                               iproc, offs_proc, offs_glob )
     9837
     9838    IMPLICIT NONE
     9839
     9840    INTEGER(iwp), INTENT(IN)                              ::  i           !< x-coordinate
     9841    INTEGER(iwp), INTENT(IN)                              ::  j           !< y-coordinate
     9842    INTEGER(iwp), INTENT(IN)                              ::  k           !< z-coordinate
     9843    INTEGER(iwp), INTENT(IN)                              ::  nlayers     !< number of z-layers
     9844    INTEGER(iwp), INTENT(OUT), OPTIONAL                   ::  iproc       !< MPI process rank
     9845    INTEGER(kind=MPI_ADDRESS_KIND), INTENT(OUT), OPTIONAL ::  offs_proc   !< offset within MPI proc
     9846    INTEGER(iwp), INTENT(OUT), OPTIONAL                   ::  offs_glob   !< global offset
     9847
     9848    INTEGER(iwp)            ::  ipx         !< process index in x-direction
     9849    INTEGER(iwp)            ::  ipy         !< process index in y-direction
     9850    INTEGER(iwp)            ::  iproc_l     !< local variable for iproc
     9851    INTEGER(iwp)            ::  oproc_l     !< local variable for offs_proc
     9852    INTEGER(iwp)            ::  offs_pstart !< global start of the MPI process
     9853
     9854    ipx = i / nnx
     9855    ipy = j / nny
     9856    iproc_l = ipx * pdims(2) + ipy
     9857    IF ( PRESENT(iproc) )  iproc = iproc_l
     9858
     9859    IF ( PRESENT(offs_proc)  .OR.  PRESENT(offs_glob) )  THEN
     9860       oproc_l = (i - ipx*nnx) * nny * nlayers + & ! columns before
     9861                 (j - ipy*nny) * nlayers       + & ! rows in column
     9862                 k
     9863       IF ( PRESENT(offs_proc) )  offs_proc = oproc_l
     9864       IF ( PRESENT(offs_glob) )  THEN
     9865          offs_pstart = iproc_l * nnx * nny * nlayers
     9866          offs_glob = offs_pstart + oproc_l
     9867       ENDIF
     9868    ENDIF
     9869
     9870 END SUBROUTINE radiation_calc_global_offset
     9871
     9872
     9873!------------------------------------------------------------------------------!
     9874!
     9875! Description:
     9876! ------------
    98569877!> Subroutine for averaging 3D data
    98579878!------------------------------------------------------------------------------!
Note: See TracChangeset for help on using the changeset viewer.