Changeset 4574 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- Jun 24, 2020 4:33:32 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r4571 r4574 28 28 ! ----------------- 29 29 ! $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 30 34 ! Bugfix in vertical lad_s coordinate 31 35 ! … … 1138 1142 IMPLICIT NONE 1139 1143 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 1146 1151 1147 1152 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 1181 1168 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', & 1185 1201 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out' ) 1186 1202 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 ' // & 1189 1205 'radiation_scheme = "rrtmg"' 1190 1206 CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 ) … … 1192 1208 unit = 'K/h' 1193 1209 1194 CASE ( 'rad_net*', 'r rtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', &1195 'r rtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*', &1196 'rad_sw_out*')1197 IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN1198 ! Workaround for masked output (calls with i=ilen=k=0)1199 unit = 'illegal'1200 RETURN1210 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 ) 1201 1217 ENDIF 1218 unit = 'W/m2' 1219 1220 CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' ) 1202 1221 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 ' // & 1205 1224 'cross sections are allowed for this value' 1206 1225 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) 1207 1226 ENDIF 1208 1227 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 ) 1219 1232 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 = '' 1231 1234 1232 1235 CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', & 1233 1236 'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref') 1234 1237 IF ( .NOT. radiation ) THEN 1235 message_string = 'output of "' // TRIM( var ) // '" require'&1238 message_string = 'output of "' // var // '" require' & 1236 1239 // 's radiation = .TRUE.' 1237 1240 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) … … 1240 1243 1241 1244 CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw' ) 1242 IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN1243 ! Workaround for masked output (calls with i=ilen=k=0)1244 unit = 'illegal'1245 RETURN1246 ENDIF1247 1248 1245 IF ( .NOT. radiation ) THEN 1249 message_string = 'output of "' // TRIM( var ) // '" require'&1246 message_string = 'output of "' // var // '" require' & 1250 1247 // 's radiation = .TRUE.' 1251 1248 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) 1252 1249 ENDIF 1253 1250 IF ( mrt_nlevels == 0 ) THEN 1254 message_string = 'output of "' // TRIM( var ) // '" require'&1251 message_string = 'output of "' // var // '" require' & 1255 1252 // 's mrt_nlevels > 0' 1256 1253 CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 ) 1257 1254 ENDIF 1258 IF ( TRIM( var )== 'rtm_mrt_sw' .AND. .NOT. mrt_include_sw ) THEN1259 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' & 1260 1257 // 's rtm_mrt_sw = .TRUE.' 1261 1258 CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 ) 1262 1259 ENDIF 1263 IF ( TRIM( var )== 'rtm_mrt' ) THEN1260 IF ( var == 'rtm_mrt' ) THEN 1264 1261 unit = 'K' 1265 1262 ELSE … … 1271 1268 1272 1269 END SELECT 1273 END IF1270 END IF 1274 1271 1275 1272 END SUBROUTINE radiation_check_data_output … … 8516 8513 INTEGER(iwp), DIMENSION(3) :: dimnext !< next dimension increments along path 8517 8514 INTEGER(iwp), DIMENSION(3) :: dimdelta !< dimension direction = +- 1 8518 INTEGER(iwp) :: px, py !< number of processors in x and y dir before8519 !< the processor in the question8520 INTEGER(iwp) :: ip !< number of processor where gridbox reside8521 8515 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array 8522 8516 … … 8583 8577 !-- calculate index of the grid with global indices (box(2),box(3)) 8584 8578 !-- 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 ) 8589 8580 IF ( box(1) <= nzterr(ig) ) THEN 8590 8581 visible = .FALSE. … … 8598 8589 crlens(ncsb) = crlen 8599 8590 #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) ) 8602 8594 #endif 8603 8595 ENDIF … … 8716 8708 INTEGER(iwp), DIMENSION(2) :: dimnext !< next dimension increments along path 8717 8709 INTEGER(iwp), DIMENSION(2) :: dimdelta !< dimension direction = +- 1 8718 INTEGER(iwp) :: px, py !< number of processors in x and y dir before8719 !< the processor in the question8720 8710 INTEGER(iwp) :: ip !< number of processor where gridbox reside 8721 8711 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array … … 8821 8811 !-- calculate index of the grid with global indices (column(1),column(2)) 8822 8812 !-- 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 ) 8827 8814 8828 8815 IF ( lastdist == 0._wp ) THEN … … 8898 8885 !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' ) 8899 8886 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 ) 8904 8889 8905 8890 IF ( rad_angular_discretization .AND. calc_svf ) THEN … … 8919 8904 IF ( plantt(ig) < lowest_lad ) CYCLE 8920 8905 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 ) 8922 8909 wcount = plantt(ig)-lowest_lad+1 8923 8910 ! TODO send request ASAP - even during raytracing … … 8942 8929 ELSE ! raytrace_mpi_rma = .F. 8943 8930 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 ) 8948 8933 rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad_g(ig:ig+nly-1) 8949 8934 ENDDO … … 9011 8996 DO i = 1, ntrack ! for each column 9012 8997 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 ) 9016 8999 9017 9000 DO k = 1, nrays ! for each ray … … 9081 9064 DO i = ntrack, 1, -1 ! for each column backwards 9082 9065 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 ) 9086 9067 9087 9068 DO k = 1, nrays ! for each ray … … 9150 9131 INTEGER(iwp), TARGET, INTENT(out) :: isurfl 9151 9132 INTEGER(iwp), INTENT(out) :: iproc 9152 #if defined( __parallel )9153 INTEGER(iwp) :: px, py !< number of processors in x and y direction9154 !< before the processor in the question9155 #endif9156 9133 9157 9134 #if defined( __parallel ) … … 9160 9137 ! 9161 9138 !-- 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 ) 9167 9142 ! 9168 9143 !-- Send MPI_Get request to obtain index target_surfl(i) … … 9854 9829 ! Description: 9855 9830 ! ------------ 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 ! ------------ 9856 9877 !> Subroutine for averaging 3D data 9857 9878 !------------------------------------------------------------------------------!
Note: See TracChangeset
for help on using the changeset viewer.
