Changeset 3337 for palm/trunk
- Timestamp:
- Oct 12, 2018 3:17:09 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 16 edited
- 1 copied
-
. (modified) (1 prop)
-
Makefile (modified) (8 diffs, 1 prop)
-
average_3d_data.f90 (modified) (5 diffs)
-
biometeorology_mod.f90 (copied) (copied from palm/branches/resler/SOURCE/biometeorology_mod.f90)
-
check_parameters.f90 (modified) (5 diffs, 1 prop)
-
chem_emissions_mod.f90 (modified) (17 diffs)
-
data_output_3d.f90 (modified) (3 diffs)
-
modules.f90 (modified) (2 diffs)
-
netcdf_data_input_mod.f90 (modified) (7 diffs)
-
netcdf_interface_mod.f90 (modified) (7 diffs)
-
palm.f90 (modified) (2 diffs)
-
plant_canopy_model_mod.f90 (modified) (4 diffs)
-
prognostic_equations.f90 (modified) (2 diffs)
-
radiation_model_mod.f90 (modified) (125 diffs)
-
sum_up_3d_data.f90 (modified) (8 diffs)
-
time_integration_spinup.f90 (modified) (3 diffs, 1 prop)
-
urban_surface_mod.f90 (modified) (70 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/resler/SOURCE (added) merged: 2136-2138,2324-2325,2655,2679,2684,2694-2695,2783-2786,2810,2985,3000,3015,3017,3060-3061,3063,3154,3317,3319-3321,3323,3335-3336
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/Makefile
- Property svn:mergeinfo changed
/palm/branches/resler/SOURCE/Makefile (added) merged: 2136-2138,2684,2783,2785,2810,2985,3000,3317,3320,3336
r3322 r3337 25 25 # ----------------- 26 26 # $Id$ 27 # (from branch resler) 28 # Add biometeorology 29 # 30 # 31 # 3322 2018-10-09 10:02:39Z kanani 27 32 # Formatting and cleanup 28 33 # … … 514 519 gust_mod.f90 \ 515 520 header.f90 \ 521 biometeorology_mod.f90 \ 516 522 inflow_turbulence.f90 \ 517 523 init_3d_model.f90 \ … … 712 718 modules.o 713 719 average_3d_data.o: \ 720 biometeorology_mod.o \ 714 721 bulk_cloud_model_mod.o \ 715 722 chemistry_model_mod.o \ … … 726 733 basic_constants_and_equations_mod.o: \ 727 734 mod_kinds.o 735 biometeorology_mod.o: \ 736 mod_kinds.o \ 737 radiation_model_mod.o 728 738 boundary_conds.o: \ 729 739 basic_constants_and_equations_mod.o \ … … 764 774 check_parameters.o: \ 765 775 basic_constants_and_equations_mod.o \ 776 biometeorology_mod.o \ 766 777 bulk_cloud_model_mod.o \ 767 778 chemistry_model_mod.o \ … … 904 915 data_output_3d.o: \ 905 916 basic_constants_and_equations_mod.o \ 917 biometeorology_mod.o \ 906 918 bulk_cloud_model_mod.o \ 907 919 chemistry_model_mod.o \ … … 1264 1276 netcdf_interface_mod.o: \ 1265 1277 basic_constants_and_equations_mod.o \ 1278 biometeorology_mod.o \ 1266 1279 chemistry_model_mod.o \ 1267 1280 gust_mod.o \ … … 1513 1526 sum_up_3d_data.o: \ 1514 1527 basic_constants_and_equations_mod.o \ 1528 biometeorology_mod.o \ 1515 1529 bulk_cloud_model_mod.o \ 1516 1530 chemistry_model_mod.o \ - Property svn:mergeinfo changed
-
palm/trunk/SOURCE/average_3d_data.f90
r3294 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from resler branch) 28 ! Add biometeorology output, 29 ! fix for chemistry output, 30 ! move of urban surface output 31 ! 32 ! 3294 2018-10-01 02:37:10Z raasch 27 33 ! changes concerning modularization of ocean option 28 34 ! … … 174 180 ONLY: radiation, radiation_3d_data_averaging 175 181 182 USE biometeorology_mod, & 183 ONLY: biometeorology_3d_data_averaging 184 176 185 USE turbulence_closure_mod, & 177 186 ONLY: tcm_3d_data_averaging … … 203 212 DO ii = 1, doav_n 204 213 205 !206 !-- Temporary solution to account for data output within the new urban207 !-- surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )208 214 trimvar = TRIM( doav(ii) ) 209 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN210 trimvar = 'usm_output'211 ENDIF212 215 213 216 ! … … 533 536 ENDIF 534 537 535 CASE ( 'usm_output' )536 !537 !-- Block of urban surface model outputs538 CALL usm_average_3d_data( 'average', doav(ii) )539 540 538 CASE DEFAULT 541 539 ! 542 540 !-- Averaging of data from other modules 543 IF ( air_chemistry ) THEN 541 IF ( air_chemistry .AND. & 542 (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') ) THEN 544 543 CALL chem_3d_data_averaging( 'average', doav(ii) ) 545 544 ENDIF … … 565 564 ENDIF 566 565 566 IF ( radiation .AND. trimvar(1:4) == 'bio_' ) THEN 567 CALL biometeorology_3d_data_averaging( 'average', doav(ii) ) 568 ENDIF 569 567 570 CALL tcm_3d_data_averaging( 'average', doav(ii) ) 571 572 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN 573 CALL usm_average_3d_data( 'average', doav(ii) ) 574 ENDIF 575 568 576 ! 569 577 !-- User-defined quantities -
palm/trunk/SOURCE/check_parameters.f90
- Property svn:mergeinfo changed
/palm/branches/resler/SOURCE/check_parameters.f90 (added) merged: 2136-2138,2324,2684,2783,2785,2810,3000,3060,3317,3320,3323
r3312 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Add biometeorology, 29 ! fix for chemistry output, 30 ! increase of uv_heights dimension 31 ! 32 ! 3312 2018-10-06 14:15:46Z knoop 27 33 ! small code rearrangements 28 34 ! … … 756 762 ONLY: radiation, radiation_check_data_output, & 757 763 radiation_check_data_output_pr, radiation_check_parameters 764 765 USE biometeorology_mod, & 766 ONLY: biometeorology_check_data_output 758 767 759 768 USE spectra_mod, & … … 1662 1671 DO k = 1, nz+1 1663 1672 1664 IF ( kk < 100 ) THEN1673 IF ( kk < 200 ) THEN 1665 1674 DO WHILE ( uv_heights(kk+1) <= zu(k) ) 1666 1675 kk = kk + 1 1667 IF ( kk == 100 ) EXIT1676 IF ( kk == 200 ) EXIT 1668 1677 ENDDO 1669 1678 ENDIF 1670 1679 1671 IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp ) THEN1680 IF ( kk < 200 .AND. uv_heights(kk+1) /= 9999999.9_wp ) THEN 1672 1681 u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) / & 1673 1682 ( uv_heights(kk+1) - uv_heights(kk) ) * & … … 3156 3165 3157 3166 IF ( unit == 'illegal' .AND. air_chemistry & 3158 .AND. var(1:3) == 'kc_') THEN3167 .AND. (var(1:3) == 'kc_' .OR. var(1:3) == 'em_') ) THEN 3159 3168 CALL chem_check_data_output( var, unit, i, ilen, k ) 3160 3169 ENDIF … … 3162 3171 IF ( unit == 'illegal' ) THEN 3163 3172 CALL lsm_check_data_output ( var, unit, i, ilen, k ) 3173 ENDIF 3174 3175 IF ( unit == 'illegal' ) THEN 3176 CALL biometeorology_check_data_output( var, unit, i, ilen, k ) 3164 3177 ENDIF 3165 3178 - Property svn:mergeinfo changed
-
palm/trunk/SOURCE/chem_emissions_mod.f90
r3312 r3337 27 27 ! ----------------- 28 28 ! $Id$ 29 ! (from branch resler) 30 ! Formatting 31 ! 32 ! 3312 2018-10-06 14:15:46Z knoop 29 33 ! Initial revision 30 34 ! … … 719 723 DO ispec = 1 , len_index 720 724 721 IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN 722 723 message_string = 'PARAMETERIZED emissions mode selected:' // & 725 IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. & 726 emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN 727 728 message_string = 'PARAMETERIZED emissions mode selected:' // & 724 729 ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // & 725 730 ' but values of scaling factors for street types' // & … … 1246 1251 IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN 1247 1252 !> Kg/m2 kg/m2*s 1248 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor 1253 delta_emis(nys:nyn,nxl:nxr) = & 1254 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor 1249 1255 1250 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1256 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1257 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1251 1258 1252 1259 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN 1253 1260 1254 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor 1255 1256 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1261 delta_emis(nys:nyn,nxl:nxr) = & 1262 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor 1263 1264 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1265 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1257 1266 1258 1267 !> SOX Compositions … … 1260 1269 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN 1261 1270 1262 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor 1263 1264 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1271 delta_emis(nys:nyn,nxl:nxr) = & 1272 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor 1273 1274 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1275 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1265 1276 1266 1277 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN 1267 1278 1268 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor 1269 1270 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1279 delta_emis(nys:nyn,nxl:nxr) = & 1280 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor 1281 1282 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1283 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1271 1284 1272 1285 … … 1278 1291 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1)) 1279 1292 1280 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor 1293 delta_emis(nys:nyn,nxl:nxr) = & 1294 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor 1281 1295 1282 1296 1283 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1297 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1298 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1284 1299 1285 1300 ENDDO … … 1291 1306 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2)) 1292 1307 1293 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor 1308 delta_emis(nys:nyn,nxl:nxr) = & 1309 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor 1294 1310 1295 1311 1296 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1312 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1313 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1297 1314 1298 1315 ENDDO … … 1304 1321 DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3)) 1305 1322 1306 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor 1323 delta_emis(nys:nyn,nxl:nxr) = & 1324 emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor 1307 1325 1308 1326 1309 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1327 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1328 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1310 1329 1311 1330 ENDDO … … 1319 1338 IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN 1320 1339 1321 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)* 1340 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)* & 1322 1341 emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor 1323 1342 1324 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1343 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1344 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1325 1345 1326 1346 ENDIF … … 1333 1353 delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor 1334 1354 1335 emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1355 emis_distribution(1,nys:nyn,nxl:nxr,ispec) = & 1356 emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 1336 1357 1337 1358 ENDIF ! IF (spc_names==) … … 1375 1396 1376 1397 !> PMs are already in mass units:micrograms:they have to be converted to kilograms 1377 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1398 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1399 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1378 1400 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1379 1401 … … 1389 1411 !> Other Species: inputs are micromoles: they have to be converted in moles 1390 1412 ! ppm/s *m *kg/m^3 1391 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* 1413 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))* & 1392 1414 ! micromoles/(m^2*s) 1393 emis_distribution(1,j,i,ispec) * 1415 emis_distribution(1,j,i,ispec) * & 1394 1416 ! m**3/Nmole 1395 conv_to_ratio(k,j,i)* &1417 conv_to_ratio(k,j,i)* & 1396 1418 ! kg/m^3 1397 1419 rho_air(k) … … 1403 1425 1404 1426 1405 ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. street_type_f%var(j,i) < main_street_id ) THEN 1427 ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. & 1428 street_type_f%var(j,i) < main_street_id ) THEN 1406 1429 1407 1430 !> Cycle over already matched species … … 1409 1432 1410 1433 !> PMs are already in mass units: micrograms 1411 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1434 IF ( TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1435 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1412 1436 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1413 1437 1414 1438 ! kg/(m^2*s) *kg/m^3 1415 surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * 1439 surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) * & 1416 1440 ! kg/(m^2*s) 1417 1441 emis_distribution(1,j,i,ispec)* & … … 1423 1447 !>Other Species: inputs are micromoles 1424 1448 ! ppm/s *m *kg/m^3 1425 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * 1449 surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * & 1426 1450 ! micromoles/(m^2*s) 1427 emis_distribution(1,j,i,ispec) * 1451 emis_distribution(1,j,i,ispec) * & 1428 1452 ! m**3/Nmole 1429 conv_to_ratio(k,j,i)* &1453 conv_to_ratio(k,j,i)* & 1430 1454 ! kg/m^3 1431 1455 rho_air(k) … … 1466 1490 1467 1491 !> PMs 1468 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1469 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1492 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1493 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1494 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1470 1495 1471 1496 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) … … 1522 1547 1523 1548 !> PMs 1524 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1525 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1549 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1550 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1551 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1526 1552 1527 1553 ! kg/(m^2*s) * kg/m^3 kg/(m^2*s) … … 1573 1599 1574 1600 !> PMs 1575 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1576 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1601 IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" & 1602 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25" & 1603 .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN 1577 1604 1578 1605 ! kg/(m^2*s) *kg/m^3 kg/(m^2*s) -
palm/trunk/SOURCE/data_output_3d.f90
r3294 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Add Biometeorology 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 269 273 USE radiation_model_mod, & 270 274 ONLY: nzub, nzut, radiation, radiation_data_output_3d 275 276 USE biometeorology_mod, & 277 ONLY: biometeorology_data_output_3d 271 278 272 279 USE turbulence_closure_mod, & … … 744 751 IF ( .NOT. found .AND. radiation ) THEN 745 752 CALL radiation_data_output_3d( av, do3d(av,if), found, & 753 local_pf, nzb_do, nzt_do ) 754 resorted = .TRUE. 755 ENDIF 756 757 ! 758 !-- Biometeorology quantity 759 IF ( .NOT. found .AND. radiation ) THEN 760 CALL biometeorology_data_output_3d( av, do3d(av,if), found, & 746 761 local_pf, nzb_do, nzt_do ) 747 762 resorted = .TRUE. -
palm/trunk/SOURCE/modules.f90
r3302 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Increase dimension of uv_heights etc. 29 ! 30 ! 3302 2018-10-03 02:39:40Z raasch 27 31 ! +u/v_stokes_zu, u/v_stokes_zw 28 32 ! … … 1566 1570 REAL(wp) :: tsc(10) = (/ 1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, & !< array used for controlling time-integration at different substeps 1567 1571 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp /) 1568 REAL(wp) :: u_profile( 100) = 9999999.9_wp !< namelist parameter1569 REAL(wp) :: uv_heights( 100) = 9999999.9_wp !< namelist parameter1570 REAL(wp) :: v_profile( 100) = 9999999.9_wp !< namelist parameter1572 REAL(wp) :: u_profile(200) = 9999999.9_wp !< namelist parameter 1573 REAL(wp) :: uv_heights(200) = 9999999.9_wp !< namelist parameter 1574 REAL(wp) :: v_profile(200) = 9999999.9_wp !< namelist parameter 1571 1575 REAL(wp) :: ug_vertical_gradient(10) = 0.0_wp !< namelist parameter 1572 1576 REAL(wp) :: ug_vertical_gradient_level(10) = -9999999.9_wp !< namelist parameter -
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r3298 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Formatting 29 ! 30 ! 3298 2018-10-02 12:21:11Z kanani 27 31 ! Modified EXPERT mode to PRE-PROCESSED mode (Russo) 28 32 ! Introduced Chemistry static netcdf file (Russo) … … 1030 1034 ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear)) 1031 1035 !Read-in Variable 1032 CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,emt_att%ncat,1,emt_att%nhoursyear) 1036 CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1, & 1037 emt_att%ncat,1,emt_att%nhoursyear) 1033 1038 1034 1039 !-- MDH … … 1038 1043 ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour)) 1039 1044 !-- Read-in Variable 1040 CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,emt_att%ncat,1,emt_att%nmonthdayhour) 1045 CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1, & 1046 emt_att%ncat,1,emt_att%nmonthdayhour) 1041 1047 1042 1048 ELSE … … 1055 1061 DO ispec=1,emt_att%nspec 1056 1062 1057 IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1)) 1063 IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) & 1064 ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1)) 1058 1065 1059 1066 ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)) … … 1061 1068 CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr) 1062 1069 1063 emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1) 1070 emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1) = & 1071 dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1) 1064 1072 1065 1073 DEALLOCATE (dum_var_3d) … … 1105 1113 1106 1114 !Allocation for the entire domain has to be done only for the first processor between all the subdomains 1107 IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) ) ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1)) 1115 IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) ) & 1116 ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1)) 1108 1117 1109 1118 !> allocate variable where to pass emission values read from netcdf … … 1114 1123 1115 1124 1116 emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1) 1125 emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1) = & 1126 dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1) 1117 1127 1118 1128 DEALLOCATE ( dum_var_4d ) -
palm/trunk/SOURCE/netcdf_interface_mod.f90
r3294 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Add biometeorology 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 610 614 ONLY: radiation, radiation_define_netcdf_grid 611 615 616 USE biometeorology_mod, & 617 ONLY: biometeorology_define_netcdf_grid 618 612 619 USE spectra_mod, & 613 620 ONLY: averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction … … 1016 1023 grid_z ) 1017 1024 ENDIF 1025 ! 1026 !-- Check for biometeorology quantities 1027 IF ( .NOT. found .AND. radiation ) THEN 1028 CALL biometeorology_define_netcdf_grid( domask(mid,av,i),& 1029 found, grid_x, grid_y,& 1030 grid_z ) 1031 ENDIF 1018 1032 1019 1033 CALL tcm_define_netcdf_grid( domask( mid,av,i), found, & … … 1545 1559 grid_y, grid_z ) 1546 1560 ENDIF 1547 1561 1548 1562 ! 1549 1563 !-- Check for radiation quantities 1550 1564 IF ( .NOT. found .AND. radiation ) THEN 1551 1565 CALL radiation_define_netcdf_grid( do3d(av,i), found, & 1566 grid_x, grid_y, & 1567 grid_z ) 1568 ENDIF 1569 1570 ! 1571 !-- Check for biometeorology quantities 1572 IF ( .NOT. found .AND. radiation ) THEN 1573 CALL biometeorology_define_netcdf_grid( do3d(av,i), found,& 1552 1574 grid_x, grid_y, & 1553 1575 grid_z ) … … 2322 2344 2323 2345 ! 2346 !-- Check for biometeorology quantities 2347 IF ( .NOT. found .AND. radiation ) THEN 2348 CALL biometeorology_define_netcdf_grid( do2d(av,i),& 2349 found, grid_x, grid_y,& 2350 grid_z ) 2351 ENDIF 2352 2353 ! 2324 2354 !-- Check for gust module quantities 2325 2355 IF ( .NOT. found .AND. gust_module_enabled ) THEN … … 3034 3064 3035 3065 ! 3066 !-- Check for biometeorology quantities 3067 IF ( .NOT. found .AND. radiation ) THEN 3068 CALL biometeorology_define_netcdf_grid( do2d(av,i), & 3069 found, & 3070 grid_x, grid_y, & 3071 grid_z ) 3072 ENDIF 3073 3074 ! 3036 3075 !-- Check for gust module quantities 3037 3076 IF ( .NOT. found .AND. gust_module_enabled ) THEN … … 3694 3733 3695 3734 ! 3735 !-- Check for biometeorology quantities 3736 IF ( .NOT. found .AND. radiation ) THEN 3737 CALL biometeorology_define_netcdf_grid( do2d(av,i), & 3738 found, & 3739 grid_x, grid_y, & 3740 grid_z ) 3741 ENDIF 3742 3743 ! 3696 3744 !-- Check for gust module quantities 3697 3745 IF ( .NOT. found .AND. gust_module_enabled ) THEN -
palm/trunk/SOURCE/palm.f90
r3298 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Fix chemistry call 29 ! 30 ! 3298 2018-10-02 12:21:11Z kanani 27 31 ! - Minor formatting (kanani) 28 32 ! - Added Call of date_and_time_init (Russo) … … 425 429 !-- --> Needs to be moved!! What is the dependency about? 426 430 ! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 427 IF ( air_chemistry .AND. do_emis) THEN431 IF ( air_chemistry ) THEN 428 432 429 433 IF ( do_emis ) THEN -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r3294 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Fix reading plant canopy over buildings 28 ! 29 ! 3294 2018-10-01 02:37:10Z raasch 27 30 ! ocean renamed ocean_mode 28 31 ! … … 1131 1134 INTEGER(iwp) :: i, j !< running index 1132 1135 INTEGER(iwp) :: nzp !< number of vertical layers of plant canopy 1136 INTEGER(iwp) :: nzpltop !< 1137 INTEGER(iwp) :: nzpl !< 1138 INTEGER(iwp) :: kk !< 1133 1139 1134 1140 REAL(wp), DIMENSION(:), ALLOCATABLE :: col !< vertical column of input data … … 1144 1150 FORM='FORMATTED', ERR=515) 1145 1151 READ(152, *, ERR=516, END=517) nzp !< read first line = number of vertical layers 1146 1152 nzpltop = MIN(nzt+1, nzb+nzp-1) 1153 nzpl = nzpltop - nzb + 1 !< no. of layers to assign 1147 1154 ALLOCATE( col(0:nzp-1) ) 1148 1155 … … 1164 1171 CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 ) 1165 1172 ENDIF 1166 lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype)1167 1173 kk = get_topography_top_index_ji( j, i, 's' ) 1174 lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype) 1168 1175 CASE DEFAULT 1169 1176 WRITE(message_string, '(a,i2,a)') & -
palm/trunk/SOURCE/prognostic_equations.f90
r3302 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Fix for chemistry call 29 ! 30 ! 3302 2018-10-03 02:39:40Z raasch 27 31 ! Stokes drift + wave breaking term added 28 32 ! … … 482 486 lsp_usr = 1 483 487 ! 484 !-- If required, calculate photolysis frequencies -485 !-- UNFINISHED: Why not before the intermediate timestep loop?486 IF ( intermediate_timestep_count == 1 ) THEN487 CALL photolysis_control488 ENDIF489 !490 488 !-- Chemical reactions 491 489 CALL cpu_log( log_point(82), '(chem react + exch_h)', 'start' ) 492 490 493 491 IF ( chem_gasphase_on ) THEN 492 ! 493 !-- If required, calculate photolysis frequencies - 494 !-- UNFINISHED: Why not before the intermediate timestep loop? 495 IF ( intermediate_timestep_count == 1 ) THEN 496 CALL photolysis_control 497 ENDIF 494 498 DO i = nxl, nxr 495 499 DO j = nys, nyn -
palm/trunk/SOURCE/radiation_model_mod.f90
r3274 r3337 15 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 16 16 ! 17 ! Copyright 2015-2018 Czech Technical University in Prague18 17 ! Copyright 2015-2018 Institute of Computer Science of the 19 18 ! Czech Academy of Sciences, Prague 19 ! Copyright 2015-2018 Czech Technical University in Prague 20 20 ! Copyright 1997-2018 Leibniz Universitaet Hannover 21 21 !------------------------------------------------------------------------------! … … 28 28 ! ----------------- 29 29 ! $Id$ 30 ! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague) 31 ! added calculation of the MRT inside the RTM module 32 ! MRT fluxes are consequently used in the new biometeorology module 33 ! for calculation of biological indices (MRT, PET) 34 ! Fixes of v. 2.5 and SVN trunk: 35 ! - proper initialization of rad_net_l 36 ! - make arrays nsurfs and surfstart TARGET to prevent some MPI problems 37 ! - initialization of arrays used in MPI one-sided operation as 1-D arrays 38 ! to prevent problems with some MPI/compiler combinations 39 ! - fix indexing of target displacement in subroutine request_itarget to 40 ! consider nzub 41 ! - fix LAD dimmension range in PCB calculation 42 ! - check ierr in all MPI calls 43 ! - use proper per-gridbox sky and diffuse irradiance 44 ! - fix shading for reflected irradiance 45 ! - clear away the residuals of "atmospheric surfaces" implementation 46 ! - fix rounding bug in raytrace_2d introduced in SVN trunk 47 ! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague) 48 ! can use angular discretization for all SVF 49 ! (i.e. reflected and emitted radiation in addition to direct and diffuse), 50 ! allowing for much better scaling wih high resoltion and/or complex terrain 51 ! - Unite array grow factors 52 ! - Fix slightly shifted terrain height in raytrace_2d 53 ! - Use more efficient MPI_Win_allocate for reverse gridsurf index 54 ! - Fix random MPI RMA bugs on Intel compilers 55 ! - Fix approx. double plant canopy sink values for reflected radiation 56 ! - Fix mostly missing plant canopy sinks for direct radiation 57 ! - Fix discretization errors for plant canopy sink in diffuse radiation 58 ! - Fix rounding errors in raytrace_2d 59 ! 60 ! 3274 2018-09-24 15:42:55Z knoop 30 61 ! Modularization of all bulk cloud physics code components 31 62 ! … … 432 463 !> @todo Output of other rrtm arrays (such as volume mixing ratios) 433 464 !> @todo Check for mis-used NINT() calls in raytrace_2d 465 !> RESULT: Original was correct (carefully verified formula), the change 466 !> to INT broke raytracing -- P. Krc 434 467 !> @todo Optimize radiation_tendency routines 435 468 !> … … 440 473 441 474 USE arrays_3d, & 442 ONLY: dzw, hyp, nc, pt, q, ql, zu, zw, exner, d_exner475 ONLY: dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner 443 476 444 477 USE basic_constants_and_equations_mod, & … … 735 768 rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d) 736 769 rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 737 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m ²)738 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m ²)770 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m) 771 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m) 739 772 740 773 REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation … … 771 804 INTEGER(iwp), PARAMETER :: ndsvf = 2 !< number of dimensions of real values in SVF 772 805 INTEGER(iwp), PARAMETER :: idsvf = 2 !< number of dimensions of integer values in SVF 773 INTEGER(iwp), PARAMETER :: ndcsf = 2!< number of dimensions of real values in CSF806 INTEGER(iwp), PARAMETER :: ndcsf = 1 !< number of dimensions of real values in CSF 774 807 INTEGER(iwp), PARAMETER :: idcsf = 2 !< number of dimensions of integer values in CSF 775 808 INTEGER(iwp), PARAMETER :: kdcsf = 4 !< number of dimensions of integer values in CSF calculation array … … 779 812 INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf 780 813 781 INTEGER(iwp), PARAMETER :: nsurf_type = 1 6!< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1814 INTEGER(iwp), PARAMETER :: nsurf_type = 10 !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1 782 815 783 816 INTEGER(iwp), PARAMETER :: iup_u = 0 !< 0 - index of urban upward surface (ground or roof) … … 794 827 INTEGER(iwp), PARAMETER :: iwest_l = 10 !< 10- index of land westward facing wall 795 828 796 INTEGER(iwp), PARAMETER :: iup_a = 11 !< 11- index of atm. cell ubward virtual surface 797 INTEGER(iwp), PARAMETER :: idown_a = 12 !< 12- index of atm. cell downward virtual surface 798 INTEGER(iwp), PARAMETER :: inorth_a = 13 !< 13- index of atm. cell northward facing virtual surface 799 INTEGER(iwp), PARAMETER :: isouth_a = 14 !< 14- index of atm. cell southward facing virtual surface 800 INTEGER(iwp), PARAMETER :: ieast_a = 15 !< 15- index of atm. cell eastward facing virtual surface 801 INTEGER(iwp), PARAMETER :: iwest_a = 16 !< 16- index of atm. cell westward facing virtual surface 802 803 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1,0, 0,0, 0,1,-1/) !< surface normal direction x indices 804 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0,0, 0,1,-1,0, 0/) !< surface normal direction y indices 805 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0,1,-1,0, 0,0, 0/) !< surface normal direction z indices 829 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/) !< surface normal direction x indices 830 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices 831 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices 806 832 !< parameter but set in the code 807 833 … … 816 842 817 843 !-- indices and sizes of urban and land surface models 818 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x] 819 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x] 844 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfl_l !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x] 845 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surf_l !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x] 846 INTEGER(iwp), DIMENSION(:,:), POINTER :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x] 847 INTEGER(iwp), DIMENSION(:,:), POINTER :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x] 820 848 INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor 821 INTEGER(iwp), DIMENSION(:), ALLOCATABLE 849 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: nsurfs !< array of number of all surfaces in individual processors 822 850 INTEGER(iwp) :: nsurf !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs) 823 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surfstart !< starts of blocks of surfaces for individual processors in array surf824 !< respective block for particular processor is surfstart[iproc ]+1 : surfstart[iproc+1]851 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfstart !< starts of blocks of surfaces for individual processors in array surf (indexed from 1) 852 !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1] 825 853 826 854 !-- block variables needed for calculation of the plant canopy model inside the urban surface model … … 835 863 836 864 !-- configuration parameters (they can be setup in PALM config) 837 LOGICAL :: rma_lad_raytrace = .FALSE. !< use MPI RMA to access LAD for raytracing (instead of global array) 838 LOGICAL :: mrt_factors = .FALSE. !< whether to generate MRT factor files during init 865 LOGICAL :: raytrace_mpi_rma = .TRUE. !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing 866 LOGICAL :: read_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be read from a file at initialization 867 LOGICAL :: write_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be written out to a file 868 LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 869 !< reflected radiation (as opposed to all mutually visible pairs) 870 INTEGER(iwp) :: mrt_nlevels = 0 !< number of vertical boxes above surface for which to calculate MRT 871 LOGICAL :: mrt_skip_roof = .TRUE. !< do not calculate MRT above roof surfaces 872 LOGICAL :: mrt_include_sw = .TRUE. !< should MRT calculation include SW radiation as well? 839 873 INTEGER(iwp) :: nrefsteps = 0 !< number of reflection steps to perform 840 874 REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha) … … 874 908 INTEGER(iwp) :: isurfs !< 875 909 REAL(wp) :: rsvf !< 876 REAL(wp) :: rtransp !<877 910 END TYPE 878 911 879 912 !-- arrays storing the values of USM 880 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of source and targetsurface for svf[isvf]913 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of target and source surface for svf[isvf] 881 914 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors for local surfaces 882 915 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection … … 892 925 INTEGER(iwp) :: ndsidir !< number of apparent solar directions used 893 926 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: dsidir_rev !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present 927 928 INTEGER(iwp) :: nmrtbl !< No. of local grid boxes for which MRT is calculated 929 INTEGER(iwp) :: nmrtf !< number of MRT factors for local processor 930 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtbl !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x] 931 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtfsurf !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf] 932 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtf !< array of MRT factors for each local MRT box 933 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtft !< array of MRT factors including transparency for each local MRT box 934 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtsky !< array of sky view factor for each local MRT box 935 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtskyt !< array of sky view factor including transparency for each local MRT box 936 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mrtdsit !< array of direct solar transparencies for each local MRT box 937 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw !< mean SW radiant flux for each MRT box 938 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw !< mean LW radiant flux for each MRT box 939 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt !< mean radiant temperature for each MRT box 940 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw_av !< time average mean SW radiant flux for each MRT box 941 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw_av !< time average mean LW radiant flux for each MRT box 942 REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt_av !< time average mean radiant temperature for each MRT box 894 943 895 944 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< array of sw radiation falling to local surface including radiation from reflections … … 920 969 TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array 921 970 TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array 971 TYPE(t_svf), DIMENSION(:), POINTER :: amrtf !< pointer to growing mrtf array 922 972 TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array 923 973 TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array 974 TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: amrtf1, amrtf2 !< realizations of mftf array 924 975 INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor 925 976 INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor 926 INTEGER(iwp) :: msvf, mcsf !< mod for swapping the growing array 927 INTEGER(iwp), PARAMETER :: gasize = 10000 !< initial size of growing arrays 977 INTEGER(iwp) :: nmrtfa !< dimmension of array allocated for storage of mrt 978 INTEGER(iwp) :: msvf, mcsf, mmrtf!< mod for swapping the growing array 979 INTEGER(iwp), PARAMETER :: gasize = 100000 !< initial size of growing arrays 980 REAL(wp), PARAMETER :: grow_factor = 1.4_wp !< growth factor of growing arrays 928 981 REAL(wp) :: dist_max_svf = -9999.0 !< maximum distance to calculate the minimum svf to be considered. It is 929 982 !< used to avoid very small SVFs resulting from too far surfaces with mutual visibility … … 932 985 !< needed only during calc_svf but must be here because it is 933 986 !< shared between subroutines calc_svf and raytrace 934 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< index of local pcb[k,j,i] 987 INTEGER(iwp), DIMENSION(:,:,:,:), POINTER :: gridsurf !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization) 988 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< reverse index of local pcbl[k,j,i] 989 INTEGER(iwp), PARAMETER :: nsurf_type_u = 6 !< number of urban surf types (used in gridsurf) 935 990 936 991 !-- temporary arrays for calculation of csf in raytracing … … 942 997 INTEGER(kind=MPI_ADDRESS_KIND), & 943 998 DIMENSION(:), ALLOCATABLE :: lad_disp !< array of displaycements of lad in local array of proc lad_ip 999 INTEGER(iwp) :: win_lad !< MPI RMA window for leaf area density 1000 INTEGER(iwp) :: win_gridsurf !< MPI RMA window for reverse grid surface index 944 1001 #endif 945 1002 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray 1003 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: target_surfl 946 1004 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: rt2_track 947 1005 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rt2_track_lad … … 1083 1141 !-- Public variables and constants / NEEDS SORTING 1084 1142 PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,& 1085 emissivity, force_radiation_call, & 1086 lat, lon, rad_net_av, radiation, radiation_scheme, rad_lw_in, & 1143 emissivity, force_radiation_call, lat, lon, & 1144 mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl, & 1145 rad_net_av, radiation, radiation_scheme, rad_lw_in, & 1087 1146 rad_lw_in_av, rad_lw_out, rad_lw_out_av, & 1088 1147 rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, & … … 1091 1150 skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,& 1092 1151 zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, & 1093 nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf, & 1152 write_svf_on_init, read_svf_on_init, & 1153 nrefsteps, dist_max_svf, nsvfl, svf, & 1094 1154 svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir, & 1095 1155 surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir, & 1096 1156 rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl, & 1097 surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,&1157 surfoutll, idir, jdir, kdir, id, iz, iy, ix, & 1098 1158 surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw, & 1099 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, &1159 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 1100 1160 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 1101 nsurf_type, nzub, nzut, nzu, pch, nsurf, & 1102 iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a, & 1161 nsurf_type, nzub, nzut, nzpt, nzu, pch, nsurf, & 1103 1162 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 1104 1163 radiation_interactions, startwall, startland, endland, endwall, & 1105 skyvf, skyvft, radiation_interactions_on, average_radiation 1164 skyvf, skyvft, radiation_interactions_on, average_radiation, npcbl, & 1165 pcbl 1106 1166 1107 1167 #if defined ( __rrtmg ) … … 1163 1223 SELECT CASE ( TRIM( var ) ) 1164 1224 1165 CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )1225 CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' ) 1166 1226 IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN 1167 1227 message_string = '"output of "' // TRIM( var ) // '" requi' // & … … 1205 1265 IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = '' 1206 1266 1267 CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' ) 1268 IF ( .NOT. radiation ) THEN 1269 message_string = 'output of "' // TRIM( var ) // '" require'& 1270 // 's radiation = .TRUE.' 1271 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) 1272 ENDIF 1273 IF ( mrt_nlevels == 0 ) THEN 1274 message_string = 'output of "' // TRIM( var ) // '" require'& 1275 // 's mrt_nlevels > 0' 1276 CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 ) 1277 ENDIF 1278 IF ( TRIM( var ) == 'rad_mrt_sw' .AND. .NOT. mrt_include_sw ) THEN 1279 message_string = 'output of "' // TRIM( var ) // '" require'& 1280 // 's rad_mrt_sw = .TRUE.' 1281 CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 ) 1282 ENDIF 1283 IF ( TRIM( var ) == 'rad_mrt' ) THEN 1284 unit = 'K' 1285 ELSE 1286 unit = 'W m-2' 1287 ENDIF 1288 1207 1289 CASE DEFAULT 1208 1290 unit = 'illegal' … … 1457 1539 ENDIF 1458 1540 ENDIF 1541 ! 1542 !-- Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented 1543 #if defined( __parallel ) 1544 IF ( rad_angular_discretization .AND. .NOT. raytrace_mpi_rma ) THEN 1545 message_string = 'rad_angular_discretization can only be used ' // & 1546 'together with raytrace_mpi_rma or when ' // & 1547 'no parallelization is applied.' 1548 CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 ) 1549 ENDIF 1550 #endif 1459 1551 1460 1552 IF ( cloud_droplets .AND. radiation_scheme == 'rrtmg' .AND. & … … 2354 2446 !-- In case averaged radiation is used, calculate mean temperature and 2355 2447 !-- liquid water mixing ratio at the urban-layer top. 2356 IF ( average_radiation ) THEN 2448 IF ( average_radiation ) THEN 2357 2449 pt1 = 0.0_wp 2358 2450 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp … … 2364 2456 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 2365 2457 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2366 IF ( bulk_cloud_model .OR. cloud_droplets ) & 2367 CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2458 IF ( ierr /= 0 ) THEN 2459 WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1 2460 FLUSH(9) 2461 ENDIF 2462 2463 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 2464 CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2465 IF ( ierr /= 0 ) THEN 2466 WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1 2467 FLUSH(9) 2468 ENDIF 2469 ENDIF 2368 2470 #else 2369 2471 pt1 = pt1_l … … 2534 2636 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 2535 2637 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2536 IF ( bulk_cloud_model .OR. cloud_droplets ) & 2638 IF ( ierr /= 0 ) THEN 2639 WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1 2640 FLUSH(9) 2641 ENDIF 2642 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 2537 2643 CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2644 IF ( ierr /= 0 ) THEN 2645 WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1 2646 FLUSH(9) 2647 ENDIF 2648 ENDIF 2538 2649 #else 2539 2650 pt1 = pt1_l … … 2756 2867 albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, & 2757 2868 constant_albedo, dt_radiation, emissivity, & 2758 lw_radiation, net_radiation, & 2869 lw_radiation, mrt_nlevels, mrt_skip_roof, & 2870 mrt_include_sw, net_radiation, & 2759 2871 radiation_scheme, skip_time_do_radiation, & 2760 2872 sw_radiation, unscheduled_radiation_calls, & 2873 read_svf_on_init, write_svf_on_init, & 2761 2874 max_raytracing_dist, min_irrf_value, & 2762 nrefsteps, mrt_factors, rma_lad_raytrace,&2875 nrefsteps, raytrace_mpi_rma, & 2763 2876 dist_max_svf, & 2764 2877 surface_reflections, svfnorm_report_thresh, & 2765 radiation_interactions_on 2878 radiation_interactions_on, & 2879 rad_angular_discretization, & 2880 raytrace_discrete_azims, & 2881 raytrace_discrete_elevs 2766 2882 2767 2883 NAMELIST /radiation_parameters/ albedo, albedo_type, albedo_lw_dir, & 2768 2884 albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, & 2769 2885 constant_albedo, dt_radiation, emissivity, & 2770 lw_radiation, net_radiation, & 2886 lw_radiation, mrt_nlevels, mrt_skip_roof, & 2887 mrt_include_sw, net_radiation, & 2771 2888 radiation_scheme, skip_time_do_radiation, & 2772 2889 sw_radiation, unscheduled_radiation_calls, & 2773 2890 max_raytracing_dist, min_irrf_value, & 2774 nrefsteps, mrt_factors, rma_lad_raytrace,&2891 nrefsteps, raytrace_mpi_rma, & 2775 2892 dist_max_svf, & 2776 2893 surface_reflections, svfnorm_report_thresh, & 2777 radiation_interactions_on 2894 radiation_interactions_on, & 2895 rad_angular_discretization, & 2896 raytrace_discrete_azims, & 2897 raytrace_discrete_elevs 2778 2898 2779 2899 line = ' ' … … 4456 4576 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll 4457 4577 INTEGER(iwp) :: isurf, isurfsrc, isvf, icsf, ipcgb 4578 INTEGER(iwp) :: imrt, imrtf 4458 4579 INTEGER(iwp) :: isd !< solar direction number 4459 4580 INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy … … 4547 4668 surfoutsl(:) = 0.0_wp !start-end 4548 4669 surfoutll(:) = 0.0_wp !start-end 4670 IF ( nmrtbl > 0 ) THEN 4671 mrtinsw(:) = 0._wp 4672 mrtinlw(:) = 0._wp 4673 ENDIF 4674 4549 4675 4550 4676 !-- Set up thermal radiation from surfaces … … 4623 4749 CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, & 4624 4750 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global 4751 IF ( ierr /= 0 ) THEN 4752 WRITE(9,*) 'Error MPI_AllGatherv1:', ierr, SIZE(surfoutll), nsurfl, & 4753 SIZE(surfoutl), nsurfs, surfstart 4754 FLUSH(9) 4755 ENDIF 4625 4756 #else 4626 4757 surfoutl(:) = surfoutll(:) !nsurf global … … 4639 4770 ENDDO 4640 4771 ENDIF 4641 4642 !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors 4643 surfinswdif(:) = rad_sw_in_diff(nyn,nxl) * skyvft(:) 4644 surfinlwdif(:) = rad_lw_in_diff(nyn,nxl) * skyvf(:) 4772 ! 4773 !-- diffuse radiation using sky view factor 4774 DO isurf = 1, nsurfl 4775 j = surfl(iy, isurf) 4776 i = surfl(ix, isurf) 4777 surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf) 4778 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 4779 ENDDO 4780 ! 4781 !-- MRT diffuse irradiance 4782 DO imrt = 1, nmrtbl 4783 j = mrtbl(iy, imrt) 4784 i = mrtbl(ix, imrt) 4785 mrtinsw(imrt) = mrtskyt(imrt) * rad_sw_in_diff(j,i) 4786 mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i) 4787 ENDDO 4645 4788 4646 4789 !-- direct radiation … … 4654 4797 isd = dsidir_rev(j, i) 4655 4798 DO isurf = 1, nsurfl 4656 surfinswdir(isurf) = rad_sw_in_dir(nyn,nxl) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0) 4799 j = surfl(iy, isurf) 4800 i = surfl(ix, isurf) 4801 surfinswdir(isurf) = rad_sw_in_dir(j,i) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0) 4802 ENDDO 4803 ! 4804 !-- MRT direct irradiance 4805 DO imrt = 1, nmrtbl 4806 j = mrtbl(iy, imrt) 4807 i = mrtbl(ix, imrt) 4808 mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) & 4809 / zenith(0) / 4._wp ! normal to sphere 4657 4810 ENDDO 4658 4811 ENDIF 4812 ! 4813 !-- MRT first pass thermal 4814 DO imrtf = 1, nmrtf 4815 imrt = mrtfsurf(1, imrtf) 4816 isurfsrc = mrtfsurf(2, imrtf) 4817 mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc) 4818 ENDDO 4659 4819 4660 4820 IF ( npcbl > 0 ) THEN … … 4674 4834 IF ( isurfsrc == -1 ) THEN 4675 4835 !-- Diffuse rad from sky. 4676 pcbinswdif(ipcgb) = csf(1,icsf) * csf(2,icsf) *rad_sw_in_diff(j,i)4836 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 4677 4837 4678 4838 !--Direct rad … … 4685 4845 * pc_abs_frac * dsitransc(ipcgb, isd) 4686 4846 ENDIF 4687 4688 EXIT ! only isurfsrc=-1 is processed here4689 4847 ENDIF 4690 4848 ENDDO … … 4722 4880 CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, & 4723 4881 surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 4882 IF ( ierr /= 0 ) THEN 4883 WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, & 4884 SIZE(surfouts), nsurfs, surfstart 4885 FLUSH(9) 4886 ENDIF 4887 4724 4888 CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, & 4725 4889 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 4890 IF ( ierr /= 0 ) THEN 4891 WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, & 4892 SIZE(surfoutl), nsurfs, surfstart 4893 FLUSH(9) 4894 ENDIF 4895 4726 4896 #else 4727 4897 surfouts = surfoutsl … … 4747 4917 IF ( isurfsrc == -1 ) CYCLE ! sky->face only in 1st pass, not here 4748 4918 4749 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc) 4919 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) 4920 ENDDO 4921 ! 4922 !-- MRT reflected 4923 DO imrtf = 1, nmrtf 4924 imrt = mrtfsurf(1, imrtf) 4925 isurfsrc = mrtfsurf(2, imrtf) 4926 mrtinsw(imrt) = mrtinsw(imrt) + mrtft(imrtf) * surfouts(isurfsrc) 4927 mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc) 4750 4928 ENDDO 4751 4929 … … 4755 4933 surfoutlw = surfoutlw + surfoutll 4756 4934 4757 ENDDO 4935 ENDDO ! refstep 4758 4936 4759 4937 !-- push heat flux absorbed by plant canopy to respective 3D arrays … … 4778 4956 ENDIF 4779 4957 ! 4958 !-- Calculate black body MRT (after all reflections) 4959 IF ( nmrtbl > 0 ) THEN 4960 IF ( mrt_include_sw ) THEN 4961 mrt(:) = ((mrtinsw(:) + mrtinlw(:)) / sigma_sb) ** .25_wp 4962 ELSE 4963 mrt(:) = (mrtinlw(:) / sigma_sb) ** .25_wp 4964 ENDIF 4965 ENDIF 4966 ! 4780 4967 !-- Transfer radiation arrays required for energy balance to the respective data types 4781 4968 DO i = 1, nsurfl … … 4791 4978 surf_usm_h%rad_net(m) = surfinsw(i) - surfoutsw(i) + & 4792 4979 surfinlw(i) - surfoutlw(i) 4980 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net(m) 4793 4981 ! 4794 4982 !-- northward-facding … … 4800 4988 surf_usm_v(0)%rad_net(m) = surfinsw(i) - surfoutsw(i) + & 4801 4989 surfinlw(i) - surfoutlw(i) 4990 surf_usm_v(0)%rad_net_l(m) = surf_usm_v(0)%rad_net(m) 4802 4991 ! 4803 4992 !-- southward-facding … … 4809 4998 surf_usm_v(1)%rad_net(m) = surfinsw(i) - surfoutsw(i) + & 4810 4999 surfinlw(i) - surfoutlw(i) 5000 surf_usm_v(1)%rad_net_l(m) = surf_usm_v(1)%rad_net(m) 4811 5001 ! 4812 5002 !-- eastward-facing … … 4818 5008 surf_usm_v(2)%rad_net(m) = surfinsw(i) - surfoutsw(i) + & 4819 5009 surfinlw(i) - surfoutlw(i) 5010 surf_usm_v(2)%rad_net_l(m) = surf_usm_v(2)%rad_net(m) 4820 5011 ! 4821 5012 !-- westward-facding … … 4827 5018 surf_usm_v(3)%rad_net(m) = surfinsw(i) - surfoutsw(i) + & 4828 5019 surfinlw(i) - surfoutlw(i) 5020 surf_usm_v(3)%rad_net_l(m) = surf_usm_v(3)%rad_net(m) 4829 5021 ! 4830 5022 !-- (2) land surfaces … … 4957 5149 #if defined( __parallel ) 4958 5150 CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 5151 IF ( ierr /= 0 ) THEN 5152 WRITE(9,*) 'Error MPI_AllReduce5:', ierr, pinswl, pinsw 5153 FLUSH(9) 5154 ENDIF 4959 5155 CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 5156 IF ( ierr /= 0 ) THEN 5157 WRITE(9,*) 'Error MPI_AllReduce6:', ierr, pinlwl, pinlw 5158 FLUSH(9) 5159 ENDIF 4960 5160 CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 5161 IF ( ierr /= 0 ) THEN 5162 WRITE(9,*) 'Error MPI_AllReduce7:', ierr, pabsswl, pabssw 5163 FLUSH(9) 5164 ENDIF 4961 5165 CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 5166 IF ( ierr /= 0 ) THEN 5167 WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pabslwl, pabslw 5168 FLUSH(9) 5169 ENDIF 4962 5170 CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 5171 IF ( ierr /= 0 ) THEN 5172 WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pemitlwl, pemitlw 5173 FLUSH(9) 5174 ENDIF 4963 5175 CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 5176 IF ( ierr /= 0 ) THEN 5177 WRITE(9,*) 'Error MPI_AllReduce9:', ierr, emiss_sum_surfl, emiss_sum_surf 5178 FLUSH(9) 5179 ENDIF 4964 5180 CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 5181 IF ( ierr /= 0 ) THEN 5182 WRITE(9,*) 'Error MPI_AllReduce10:', ierr, area_surfl, area_surf 5183 FLUSH(9) 5184 ENDIF 4965 5185 #else 4966 5186 pinsw = pinswl … … 5038 5258 ! t_rad_urb = pt_surf_urb / count_surfaces 5039 5259 5260 5040 5261 CONTAINS 5041 5262 … … 5201 5422 INTEGER(iwp) :: i, j, k, l, m 5202 5423 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 5203 INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb 5424 INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt 5204 5425 REAL(wp) :: mrl 5426 #if defined( __parallel ) 5427 INTEGER(iwp), DIMENSION(:), POINTER :: gridsurf_rma !< fortran pointer, but lower bounds are 1 5428 TYPE(c_ptr) :: gridsurf_rma_p !< allocated c pointer 5429 INTEGER(iwp) :: minfo !< MPI RMA window info handle 5430 #endif 5205 5431 5206 5432 … … 5263 5489 #if defined( __parallel ) 5264 5490 CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr ) 5491 IF ( ierr /= 0 ) THEN 5492 WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nzub 5493 FLUSH(9) 5494 ENDIF 5265 5495 CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5496 IF ( ierr /= 0 ) THEN 5497 WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nzut 5498 FLUSH(9) 5499 ENDIF 5266 5500 CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5501 IF ( ierr /= 0 ) THEN 5502 WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nzpt 5503 FLUSH(9) 5504 ENDIF 5267 5505 #else 5268 5506 nzub = nzubl … … 5353 5591 !-- cycles must not be altered, certain file input routines may depend 5354 5592 !-- on it) 5355 ALLOCATE(surfl(5,nsurfl)) ! is it mecessary to allocate it with (5,nsurfl)? 5593 ALLOCATE(surfl_l(5*nsurfl)) ! is it necessary to allocate it with (5,nsurfl)? 5594 surfl(1:5,1:nsurfl) => surfl_l(1:5*nsurfl) 5356 5595 isurf = 0 5596 IF ( rad_angular_discretization ) THEN 5597 ! 5598 !-- Allocate and fill the reverse indexing array gridsurf 5599 #if defined( __parallel ) 5600 ! 5601 !-- raytrace_mpi_rma is asserted 5602 5603 CALL MPI_Info_create(minfo, ierr) 5604 IF ( ierr /= 0 ) THEN 5605 WRITE(9,*) 'Error MPI_Info_create1:', ierr 5606 FLUSH(9) 5607 ENDIF 5608 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr) 5609 IF ( ierr /= 0 ) THEN 5610 WRITE(9,*) 'Error MPI_Info_set1:', ierr 5611 FLUSH(9) 5612 ENDIF 5613 CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr) 5614 IF ( ierr /= 0 ) THEN 5615 WRITE(9,*) 'Error MPI_Info_set2:', ierr 5616 FLUSH(9) 5617 ENDIF 5618 CALL MPI_Info_set(minfo, 'same_size', 'true', ierr) 5619 IF ( ierr /= 0 ) THEN 5620 WRITE(9,*) 'Error MPI_Info_set3:', ierr 5621 FLUSH(9) 5622 ENDIF 5623 CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr) 5624 IF ( ierr /= 0 ) THEN 5625 WRITE(9,*) 'Error MPI_Info_set4:', ierr 5626 FLUSH(9) 5627 ENDIF 5628 5629 CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx, & 5630 kind=MPI_ADDRESS_KIND), & 5631 INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND), & 5632 minfo, comm2d, gridsurf_rma_p, win_gridsurf, ierr) 5633 IF ( ierr /= 0 ) THEN 5634 WRITE(9,*) 'Error MPI_Win_allocate1:', ierr, & 5635 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx,kind=MPI_ADDRESS_KIND), & 5636 INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND), win_gridsurf 5637 FLUSH(9) 5638 ENDIF 5639 5640 CALL MPI_Info_free(minfo, ierr) 5641 IF ( ierr /= 0 ) THEN 5642 WRITE(9,*) 'Error MPI_Info_free1:', ierr 5643 FLUSH(9) 5644 ENDIF 5645 5646 ! 5647 !-- On Intel compilers, calling c_f_pointer to transform a C pointer 5648 !-- directly to a multi-dimensional Fotran pointer leads to strange 5649 !-- errors on dimension boundaries. However, transforming to a 1D 5650 !-- pointer and then redirecting a multidimensional pointer to it works 5651 !-- fine. 5652 CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nzu*nny*nnx /)) 5653 gridsurf(0:nsurf_type_u-1, nzub:nzut, nys:nyn, nxl:nxr) => & 5654 gridsurf_rma(1:nsurf_type_u*nzu*nny*nnx) 5655 #else 5656 ALLOCATE(gridsurf(0:nsurf_type_u-1,nzub:nzut,nys:nyn,nxl:nxr) ) 5657 #endif 5658 gridsurf(:,:,:,:) = -999 5659 ENDIF 5357 5660 5358 5661 !-- add horizontal surface elements (land and urban surfaces) … … 5362 5665 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 5363 5666 k = surf_usm_h%k(m) 5364 5365 5667 isurf = isurf + 1 5366 5668 surfl(:,isurf) = (/iup_u,k,j,i,m/) 5669 IF ( rad_angular_discretization ) THEN 5670 gridsurf(iup_u,k,j,i) = isurf 5671 ENDIF 5367 5672 ENDDO 5368 5673 5369 5674 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 5370 5675 k = surf_lsm_h%k(m) 5371 5372 5676 isurf = isurf + 1 5373 5677 surfl(:,isurf) = (/iup_l,k,j,i,m/) 5678 IF ( rad_angular_discretization ) THEN 5679 gridsurf(iup_u,k,j,i) = isurf 5680 ENDIF 5374 5681 ENDDO 5375 5682 … … 5384 5691 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 5385 5692 k = surf_usm_v(l)%k(m) 5386 5387 isurf = isurf + 1 5693 isurf = isurf + 1 5388 5694 surfl(:,isurf) = (/inorth_u,k,j,i,m/) 5695 IF ( rad_angular_discretization ) THEN 5696 gridsurf(inorth_u,k,j,i) = isurf 5697 ENDIF 5389 5698 ENDDO 5390 5699 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 5391 5700 k = surf_lsm_v(l)%k(m) 5392 5393 isurf = isurf + 1 5701 isurf = isurf + 1 5394 5702 surfl(:,isurf) = (/inorth_l,k,j,i,m/) 5703 IF ( rad_angular_discretization ) THEN 5704 gridsurf(inorth_u,k,j,i) = isurf 5705 ENDIF 5395 5706 ENDDO 5396 5707 … … 5398 5709 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 5399 5710 k = surf_usm_v(l)%k(m) 5400 5401 isurf = isurf + 1 5711 isurf = isurf + 1 5402 5712 surfl(:,isurf) = (/isouth_u,k,j,i,m/) 5713 IF ( rad_angular_discretization ) THEN 5714 gridsurf(isouth_u,k,j,i) = isurf 5715 ENDIF 5403 5716 ENDDO 5404 5717 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 5405 5718 k = surf_lsm_v(l)%k(m) 5406 5407 isurf = isurf + 1 5719 isurf = isurf + 1 5408 5720 surfl(:,isurf) = (/isouth_l,k,j,i,m/) 5721 IF ( rad_angular_discretization ) THEN 5722 gridsurf(isouth_u,k,j,i) = isurf 5723 ENDIF 5409 5724 ENDDO 5410 5725 … … 5412 5727 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 5413 5728 k = surf_usm_v(l)%k(m) 5414 5415 isurf = isurf + 1 5729 isurf = isurf + 1 5416 5730 surfl(:,isurf) = (/ieast_u,k,j,i,m/) 5731 IF ( rad_angular_discretization ) THEN 5732 gridsurf(ieast_u,k,j,i) = isurf 5733 ENDIF 5417 5734 ENDDO 5418 5735 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 5419 5736 k = surf_lsm_v(l)%k(m) 5420 5421 isurf = isurf + 1 5737 isurf = isurf + 1 5422 5738 surfl(:,isurf) = (/ieast_l,k,j,i,m/) 5739 IF ( rad_angular_discretization ) THEN 5740 gridsurf(ieast_u,k,j,i) = isurf 5741 ENDIF 5423 5742 ENDDO 5424 5743 … … 5426 5745 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 5427 5746 k = surf_usm_v(l)%k(m) 5428 5429 isurf = isurf + 1 5747 isurf = isurf + 1 5430 5748 surfl(:,isurf) = (/iwest_u,k,j,i,m/) 5749 IF ( rad_angular_discretization ) THEN 5750 gridsurf(iwest_u,k,j,i) = isurf 5751 ENDIF 5431 5752 ENDDO 5432 5753 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 5433 5754 k = surf_lsm_v(l)%k(m) 5434 5435 isurf = isurf + 1 5755 isurf = isurf + 1 5436 5756 surfl(:,isurf) = (/iwest_l,k,j,i,m/) 5757 IF ( rad_angular_discretization ) THEN 5758 gridsurf(iwest_u,k,j,i) = isurf 5759 ENDIF 5437 5760 ENDDO 5438 5761 ENDDO 5439 5762 ENDDO 5763 ! 5764 !-- Add local MRT boxes for specified number of levels 5765 nmrtbl = 0 5766 IF ( mrt_nlevels > 0 ) THEN 5767 DO i = nxl, nxr 5768 DO j = nys, nyn 5769 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 5770 ! 5771 !-- Skip roof if requested 5772 IF ( mrt_skip_roof .AND. surf_usm_h%isroof_surf(m) ) CYCLE 5773 ! 5774 !-- Cycle over specified no of levels 5775 nmrtbl = nmrtbl + mrt_nlevels 5776 ENDDO 5777 ! 5778 !-- Dtto for LSM 5779 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 5780 nmrtbl = nmrtbl + mrt_nlevels 5781 ENDDO 5782 ENDDO 5783 ENDDO 5784 5785 ALLOCATE( mrtbl(iz:ix,nmrtbl), mrtsky(nmrtbl), mrtskyt(nmrtbl), & 5786 mrtinsw(nmrtbl), mrtinlw(nmrtbl), mrt(nmrtbl) ) 5787 5788 imrt = 0 5789 DO i = nxl, nxr 5790 DO j = nys, nyn 5791 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 5792 ! 5793 !-- Skip roof if requested 5794 IF ( mrt_skip_roof .AND. surf_usm_h%isroof_surf(m) ) CYCLE 5795 ! 5796 !-- Cycle over specified no of levels 5797 l = surf_usm_h%k(m) 5798 DO k = l, l + mrt_nlevels - 1 5799 imrt = imrt + 1 5800 mrtbl(:,imrt) = (/k,j,i/) 5801 ENDDO 5802 ENDDO 5803 ! 5804 !-- Dtto for LSM 5805 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 5806 l = surf_lsm_h%k(m) 5807 DO k = l, l + mrt_nlevels - 1 5808 imrt = imrt + 1 5809 mrtbl(:,imrt) = (/k,j,i/) 5810 ENDDO 5811 ENDDO 5812 ENDDO 5813 ENDDO 5814 ENDIF 5440 5815 5441 5816 ! … … 5458 5833 #if defined( __parallel ) 5459 5834 CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr) 5835 IF ( ierr /= 0 ) THEN 5836 WRITE(9,*) 'Error MPI_AllGather1:', ierr, nsurfl, nsurfs 5837 FLUSH(9) 5838 ENDIF 5839 5460 5840 #else 5461 5841 nsurfs(0) = nsurfl … … 5469 5849 surfstart(numprocs) = k 5470 5850 nsurf = k 5471 ALLOCATE(surf(5,nsurf)) 5851 ALLOCATE(surf_l(5*nsurf)) 5852 surf(1:5,1:nsurf) => surf_l(1:5*nsurf) 5472 5853 5473 5854 #if defined( __parallel ) 5474 CALL MPI_AllGatherv(surfl , nsurfl*5, MPI_INTEGER, surf, nsurfs*5, &5855 CALL MPI_AllGatherv(surfl_l, nsurfl*5, MPI_INTEGER, surf_l, nsurfs*5, & 5475 5856 surfstart(0:numprocs-1)*5, MPI_INTEGER, comm2d, ierr) 5857 IF ( ierr /= 0 ) THEN 5858 WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_l), nsurfl*5, & 5859 SIZE(surf_l), nsurfs*5, surfstart(0:numprocs-1)*5 5860 FLUSH(9) 5861 ENDIF 5476 5862 #else 5477 5863 surf = surfl … … 5534 5920 5535 5921 INTEGER(iwp) :: i, j, k, d, ip, jp 5536 INTEGER(iwp) :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt t, imrtf, ipcgb5922 INTEGER(iwp) :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt, imrtf, ipcgb 5537 5923 INTEGER(iwp) :: sd, td 5538 5924 INTEGER(iwp) :: iaz, izn !< azimuth, zenith counters … … 5542 5928 REAL(wp) :: az1, az2 !< relative azimuth of section borders 5543 5929 REAL(wp) :: azmid !< ray (center) azimuth 5544 REAL(wp) :: horizon !< computed horizon height (tangent of elevation)5545 REAL(wp) :: azen !< zenith angle5546 5930 REAL(wp), DIMENSION(2) :: yxdir !< y,x *unit* vector of ray direction (in grid units) 5547 5931 REAL(wp), DIMENSION(:), ALLOCATABLE :: zdirs !< directions in z (tangent of elevation) 5548 5932 REAL(wp), DIMENSION(:), ALLOCATABLE :: zbdry !< zenith angle boundaries 5549 5933 REAL(wp), DIMENSION(:), ALLOCATABLE :: vffrac !< view factor fractions for individual rays 5934 REAL(wp), DIMENSION(:), ALLOCATABLE :: vffrac0 !< dtto (original values) 5550 5935 REAL(wp), DIMENSION(:), ALLOCATABLE :: ztransp !< array of transparency in z steps 5936 INTEGER(iwp) :: lowest_free_ray !< index into zdirs 5937 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: itarget !< face indices of detected obstacles 5938 INTEGER(iwp) :: itarg0, itarg1 5939 #if defined( __parallel ) 5940 #endif 5941 5942 5943 5551 5944 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 5552 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzterrl 5553 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csflt, pcsflt 5554 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: kcsflt,kpcsflt 5945 INTEGER(iwp) :: udim 5946 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l 5947 INTEGER(iwp), DIMENSION(:,:), POINTER :: nzterrl 5948 REAL(wp), DIMENSION(:), ALLOCATABLE,TARGET:: csflt_l, pcsflt_l 5949 REAL(wp), DIMENSION(:,:), POINTER :: csflt, pcsflt 5950 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: kcsflt_l,kpcsflt_l 5951 INTEGER(iwp), DIMENSION(:,:), POINTER :: kcsflt,kpcsflt 5555 5952 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: icsflt,dcsflt,ipcsflt,dpcsflt 5556 5953 REAL(wp), DIMENSION(3) :: uv … … 5561 5958 INTEGER(idp) :: ray_skip_maxdist, ray_skip_minval !< skipped raytracing counts 5562 5959 INTEGER(iwp) :: max_track_len !< maximum 2d track length 5563 INTEGER(iwp) :: win_lad,minfo5564 REAL(wp), DIMENSION(: ,:,:), POINTER :: lad_s_rma !< fortran pointer, but lower bounds are 15960 INTEGER(iwp) :: minfo 5961 REAL(wp), DIMENSION(:), POINTER :: lad_s_rma !< fortran 1D pointer 5565 5962 TYPE(c_ptr) :: lad_s_rma_p !< allocated c pointer 5566 5963 #if defined( __parallel ) … … 5569 5966 ! 5570 5967 INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts 5571 !CHARACTER(200) :: msg5968 CHARACTER(200) :: msg 5572 5969 5573 5970 !-- calculation of the SVF 5574 5971 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 5575 !CALL radiation_write_debug_log('Start calculation of SVF and CSF')5972 CALL radiation_write_debug_log('Start calculation of SVF and CSF') 5576 5973 5577 5974 !-- precalculate face areas for different face directions using normal vector … … 5596 5993 acsf => acsf1 5597 5994 ENDIF 5995 nmrtf = 0 5996 IF ( mrt_nlevels > 0 ) THEN 5997 nmrtfa = gasize 5998 mmrtf = 1 5999 ALLOCATE ( amrtf1(nmrtfa) ) 6000 amrtf => amrtf1 6001 ENDIF 5598 6002 ray_skip_maxdist = 0 5599 6003 ray_skip_minval = 0 … … 5602 6006 ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) ) 5603 6007 #if defined( __parallel ) 5604 ALLOCATE( nzterrl(nys:nyn,nxl:nxr) ) 6008 !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) ) 6009 ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) ) 6010 nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1)) 5605 6011 nzterrl = get_topography_top_index( 's' ) 5606 CALL MPI_AllGather( nzterrl , nnx*nny, MPI_INTEGER, &6012 CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, & 5607 6013 nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr ) 5608 DEALLOCATE(nzterrl) 6014 IF ( ierr /= 0 ) THEN 6015 WRITE(9,*) 'Error MPI_AllGather1:', ierr, SIZE(nzterrl_l), nnx*nny, & 6016 SIZE(nzterr), nnx*nny 6017 FLUSH(9) 6018 ENDIF 6019 DEALLOCATE(nzterrl_l) 5609 6020 #else 5610 6021 nzterr = RESHAPE( get_topography_top_index( 's' ), (/(nx+1)*(ny+1)/) ) … … 5621 6032 CALL MPI_AllGather( pct, nnx*nny, MPI_INTEGER, & 5622 6033 plantt, nnx*nny, MPI_INTEGER, comm2d, ierr ) 5623 6034 IF ( ierr /= 0 ) THEN 6035 WRITE(9,*) 'Error MPI_AllGather2:', ierr, SIZE(pct), nnx*nny, & 6036 SIZE(plantt), nnx*nny 6037 FLUSH(9) 6038 ENDIF 6039 5624 6040 !-- temporary arrays storing values for csf calculation during raytracing 5625 6041 ALLOCATE( lad_ip(maxboxesg) ) 5626 6042 ALLOCATE( lad_disp(maxboxesg) ) 5627 6043 5628 IF ( r ma_lad_raytrace) THEN6044 IF ( raytrace_mpi_rma ) THEN 5629 6045 ALLOCATE( lad_s_ray(maxboxesg) ) 5630 6046 5631 6047 ! set conditions for RMA communication 5632 6048 CALL MPI_Info_create(minfo, ierr) 6049 IF ( ierr /= 0 ) THEN 6050 WRITE(9,*) 'Error MPI_Info_create2:', ierr 6051 FLUSH(9) 6052 ENDIF 5633 6053 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr) 6054 IF ( ierr /= 0 ) THEN 6055 WRITE(9,*) 'Error MPI_Info_set5:', ierr 6056 FLUSH(9) 6057 ENDIF 5634 6058 CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr) 6059 IF ( ierr /= 0 ) THEN 6060 WRITE(9,*) 'Error MPI_Info_set6:', ierr 6061 FLUSH(9) 6062 ENDIF 5635 6063 CALL MPI_Info_set(minfo, 'same_size', 'true', ierr) 6064 IF ( ierr /= 0 ) THEN 6065 WRITE(9,*) 'Error MPI_Info_set7:', ierr 6066 FLUSH(9) 6067 ENDIF 5636 6068 CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr) 6069 IF ( ierr /= 0 ) THEN 6070 WRITE(9,*) 'Error MPI_Info_set8:', ierr 6071 FLUSH(9) 6072 ENDIF 5637 6073 5638 6074 !-- Allocate and initialize the MPI RMA window … … 5643 6079 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, & 5644 6080 lad_s_rma_p, win_lad, ierr) 5645 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp, nny, nnx /)) 5646 sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:) 6081 IF ( ierr /= 0 ) THEN 6082 WRITE(9,*) 'Error MPI_Win_allocate2:', ierr, size_lad_rma, & 6083 STORAGE_SIZE(1.0_wp)/8, win_lad 6084 FLUSH(9) 6085 ENDIF 6086 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp*nny*nnx /)) 6087 sub_lad(nzub:nzpt, nys:nyn, nxl:nxr) => lad_s_rma(1:nzp*nny*nnx) 5647 6088 ELSE 5648 6089 ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr)) … … 5666 6107 5667 6108 #if defined( __parallel ) 5668 IF ( r ma_lad_raytrace) THEN6109 IF ( raytrace_mpi_rma ) THEN 5669 6110 CALL MPI_Info_free(minfo, ierr) 6111 IF ( ierr /= 0 ) THEN 6112 WRITE(9,*) 'Error MPI_Info_free2:', ierr 6113 FLUSH(9) 6114 ENDIF 5670 6115 CALL MPI_Win_lock_all(0, win_lad, ierr) 6116 IF ( ierr /= 0 ) THEN 6117 WRITE(9,*) 'Error MPI_Win_lock_all1:', ierr, win_lad 6118 FLUSH(9) 6119 ENDIF 6120 5671 6121 ELSE 5672 6122 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) ) 5673 6123 CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, & 5674 6124 sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr ) 6125 IF ( ierr /= 0 ) THEN 6126 WRITE(9,*) 'Error MPI_AllGather3:', ierr, SIZE(sub_lad), & 6127 nnx*nny*nzp, SIZE(sub_lad_g), nnx*nny*nzp 6128 FLUSH(9) 6129 ENDIF 5675 6130 ENDIF 5676 6131 #endif 5677 6132 ENDIF 5678 6133 5679 IF ( mrt_factors ) THEN 5680 OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', & 5681 action='READ', status='OLD', form='FORMATTED', err=524) 5682 OPEN(154, file='MRT_FACTORS'//myid_char, access='DIRECT', recl=(5*4+2*8), & 5683 action='WRITE', status='REPLACE', form='UNFORMATTED', err=525) 5684 imrtf = 1 5685 DO 5686 READ(153, *, end=526, err=524) imrtt, i, j, k 5687 IF ( i < nxl .OR. i > nxr & 5688 .OR. j < nys .OR. j > nyn ) CYCLE 5689 ta = (/ REAL(k), REAL(j), REAL(i) /) 5690 5691 DO isurfs = 1, nsurf 5692 IF ( .NOT. surface_facing(i, j, k, -1, & 5693 surf(ix, isurfs), surf(iy, isurfs), & 5694 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5695 CYCLE 5696 ENDIF 5697 5698 sd = surf(id, isurfs) 5699 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 5700 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 5701 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 5702 5703 !-- unit vector source -> target 5704 uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5705 sqdist = SUM(uv(:)**2) 5706 uv = uv / SQRT(sqdist) 5707 5708 !-- irradiance factor - see svf. Here we consider that target face is always normal, 5709 !-- i.e. the second dot product equals 1 5710 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & 5711 / (pi * sqdist) * facearea(sd) 5712 5713 !-- raytrace while not creating any canopy sink factors 5714 CALL raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., & 5715 visible, transparency, win_lad) 5716 IF ( .NOT. visible ) CYCLE 5717 5718 !rsvf = rirrf * transparency 5719 WRITE(154, rec=imrtf, err=525) INT(imrtt, kind=4), & 5720 INT(surf(id, isurfs), kind=4), & 5721 INT(surf(iz, isurfs), kind=4), & 5722 INT(surf(iy, isurfs), kind=4), & 5723 INT(surf(ix, isurfs), kind=4), & 5724 REAL(rirrf, kind=8), REAL(transparency, kind=8) 5725 imrtf = imrtf + 1 5726 5727 ENDDO !< isurfs 5728 ENDDO !< MRT_TARGETS record 5729 5730 524 message_string = 'error reading file MRT_TARGETS' 5731 CALL message( 'radiation_calc_svf', 'PA0524', 1, 2, 0, 6, 0 ) 5732 5733 525 message_string = 'error writing file MRT_FACTORS'//myid_char 5734 CALL message( 'radiation_calc_svf', 'PA0525', 1, 2, 0, 6, 0 ) 5735 5736 526 CLOSE(153) 5737 CLOSE(154) 5738 ENDIF !< mrt_factors 5739 6134 !-- prepare the MPI_Win for collecting the surface indices 6135 !-- from the reverse index arrays gridsurf from processors of target surfaces 6136 #if defined( __parallel ) 6137 IF ( rad_angular_discretization ) THEN 6138 ! 6139 !-- raytrace_mpi_rma is asserted 6140 CALL MPI_Win_lock_all(0, win_gridsurf, ierr) 6141 IF ( ierr /= 0 ) THEN 6142 WRITE(9,*) 'Error MPI_Win_lock_all2:', ierr, win_gridsurf 6143 FLUSH(9) 6144 ENDIF 6145 ENDIF 6146 #endif 6147 6148 5740 6149 !--Directions opposite to face normals are not even calculated, 5741 6150 !--they must be preset to 0 … … 5798 6207 END SELECT 5799 6208 5800 ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn), ztransp(1:nzn) ) 6209 ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), & 6210 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6211 !in case of rad_angular_discretization 6212 6213 itarg0 = 1 6214 itarg1 = nzn 5801 6215 zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/) 5802 6216 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 5803 6217 IF ( td == iup_u .OR. td == iup_l ) THEN 5804 !-- For horizontal target, vf fractions are constant per azimuth 5805 vffrac(:) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp) 5806 !--sum of vffrac for all iaz equals 1, verified 6218 vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp) 6219 ! 6220 !-- For horizontal target, vf fractions are constant per azimuth 6221 DO iaz = 1, naz-1 6222 vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn) 6223 ENDDO 6224 !-- sum of whole vffrac equals 1, verified 5807 6225 ENDIF 5808 5809 !--Calculate sky-view factor and direct solar visibility using 2D raytracing6226 ! 6227 !-- Calculate sky-view factor and direct solar visibility using 2D raytracing 5810 6228 DO iaz = 1, naz 5811 6229 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs … … 5814 6232 az1 = az2 - azs 5815 6233 !TODO precalculate after 1st line 5816 vffrac( :) = (SIN(az2) - SIN(az1))&6234 vffrac(itarg0:itarg1) = (SIN(az2) - SIN(az1)) & 5817 6235 * (zbdry(1:nzn) - zbdry(0:nzn-1) & 5818 6236 + SIN(zbdry(0:nzn-1))*COS(zbdry(0:nzn-1)) & 5819 6237 - SIN(zbdry(1:nzn))*COS(zbdry(1:nzn))) & 5820 6238 / (2._wp * pi) 5821 !--sum of vffrac for all iazequals 1, verified6239 !-- sum of whole vffrac equals 1, verified 5822 6240 ENDIF 5823 6241 yxdir = (/ COS(azmid), SIN(azmid) /) 5824 CALL raytrace_2d(ta, yxdir, zdirs, & 5825 surfstart(myid) + isurflt, facearea(td), & 5826 vffrac, .TRUE., .FALSE., win_lad, horizon, & 5827 ztransp) 5828 5829 5830 azen = pi/2 - ATAN(horizon) 5831 IF ( td == iup_u .OR. td == iup_l ) THEN 5832 azen = MIN(azen, pi/2) !only above horizontal direction 5833 skyvf(isurflt) = skyvf(isurflt) + (1._wp - COS(2*azen)) / & 5834 (2._wp * raytrace_discrete_azims) 5835 ELSE 5836 skyvf(isurflt) = skyvf(isurflt) + (SIN(az2) - SIN(az1)) * & 5837 (azen - SIN(azen)*COS(azen)) / (2._wp*pi) 5838 ENDIF 5839 skyvft(isurflt) = skyvft(isurflt) + SUM(ztransp(:) * vffrac(:)) 6242 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6243 surfstart(myid) + isurflt, facearea(td), & 6244 vffrac(itarg0:itarg1), .TRUE., .TRUE., & 6245 .FALSE., lowest_free_ray, & 6246 ztransp(itarg0:itarg1), & 6247 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6248 !FIXME itarget available only in 6249 !case of rad_angular_discretization 6250 skyvf(isurflt) = skyvf(isurflt) + & 6251 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) 6252 skyvft(isurflt) = skyvft(isurflt) + & 6253 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) & 6254 * vffrac(itarg0:itarg0+lowest_free_ray-1)) 5840 6255 5841 !--Save direct solar transparency6256 !-- Save direct solar transparency 5842 6257 j = MODULO(NINT(azmid/ & 5843 6258 (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), & … … 5846 6261 DO k = 1, raytrace_discrete_elevs/2 5847 6262 i = dsidir_rev(k-1, j) 5848 IF ( i /= -1 ) dsitrans(isurflt, i) = ztransp(k) 6263 IF ( i /= -1 .AND. k <= lowest_free_ray ) & 6264 dsitrans(isurflt, i) = ztransp(itarg0+k-1) 5849 6265 ENDDO 6266 6267 ! 6268 !--Advance itarget indices 6269 itarg0 = itarg1 + 1 6270 itarg1 = itarg1 + nzn 5850 6271 ENDDO 5851 6272 5852 DEALLOCATE ( zdirs, zbdry, vffrac, ztransp ) 5853 ! 5854 !-- Following calculations only required for surface_reflections 5855 IF ( surface_reflections ) THEN 5856 5857 DO isurfs = 1, nsurf 5858 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 5859 surfl(iz, isurflt), surfl(id, isurflt), & 5860 surf(ix, isurfs), surf(iy, isurfs), & 5861 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5862 CYCLE 6273 IF ( rad_angular_discretization ) THEN 6274 !-- sort itarget by face id 6275 CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz) 6276 ! 6277 !-- find the first valid position 6278 itarg0 = 1 6279 DO WHILE ( itarg0 <= nzn*naz ) 6280 IF ( itarget(itarg0) /= -1 ) EXIT 6281 itarg0 = itarg0 + 1 6282 ENDDO 6283 6284 DO i = itarg0, nzn*naz 6285 ! 6286 !-- For duplicate values, only sum up vf fraction value 6287 IF ( i < nzn*naz ) THEN 6288 IF ( itarget(i+1) == itarget(i) ) THEN 6289 vffrac(i+1) = vffrac(i+1) + vffrac(i) 6290 CYCLE 6291 ENDIF 5863 6292 ENDIF 5864 5865 sd = surf(id, isurfs) 5866 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 5867 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 5868 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 5869 5870 !-- unit vector source -> target 5871 uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5872 sqdist = SUM(uv(:)**2) 5873 uv = uv / SQRT(sqdist) 5874 5875 !-- reject raytracing above max distance 5876 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN 5877 ray_skip_maxdist = ray_skip_maxdist + 1 5878 CYCLE 5879 ENDIF 5880 5881 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 5882 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 5883 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 5884 / (pi * sqdist) & ! square of distance between centers 5885 * facearea(sd) 5886 5887 !-- reject raytracing for potentially too small view factor values 5888 IF ( rirrf < min_irrf_value ) THEN 5889 ray_skip_minval = ray_skip_minval + 1 5890 CYCLE 5891 ENDIF 5892 5893 !-- raytrace + process plant canopy sinks within 5894 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., & 5895 visible, transparency, win_lad) 5896 5897 IF ( .NOT. visible ) CYCLE 5898 ! rsvf = rirrf * transparency 5899 6293 ! 5900 6294 !-- write to the svf array 5901 6295 nsvfl = nsvfl + 1 5902 6296 !-- check dimmension of asvf array and enlarge it if needed 5903 6297 IF ( nsvfla < nsvfl ) THEN 5904 k = nsvfla * 26298 k = CEILING(REAL(nsvfla, kind=wp) * grow_factor) 5905 6299 IF ( msvf == 0 ) THEN 5906 6300 msvf = 1 … … 5917 6311 ENDIF 5918 6312 5919 ! WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 5920 ! CALL radiation_write_debug_log( msg ) 6313 WRITE (msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 6314 CALL radiation_write_debug_log( msg ) 6315 6316 nsvfla = k 6317 ENDIF 6318 !-- write svf values into the array 6319 asvf(nsvfl)%isurflt = isurflt 6320 asvf(nsvfl)%isurfs = itarget(i) 6321 asvf(nsvfl)%rsvf = vffrac(i) 6322 asvf(nsvfl)%rtransp = ztransp(i) 6323 END DO 6324 6325 ENDIF ! rad_angular_discretization 6326 6327 DEALLOCATE ( zdirs, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only 6328 !in case of rad_angular_discretization 6329 ! 6330 !-- Following calculations only required for surface_reflections 6331 IF ( surface_reflections .AND. .NOT. rad_angular_discretization ) THEN 6332 6333 DO isurfs = 1, nsurf 6334 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 6335 surfl(iz, isurflt), surfl(id, isurflt), & 6336 surf(ix, isurfs), surf(iy, isurfs), & 6337 surf(iz, isurfs), surf(id, isurfs)) ) THEN 6338 CYCLE 6339 ENDIF 6340 6341 sd = surf(id, isurfs) 6342 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 6343 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 6344 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 6345 6346 !-- unit vector source -> target 6347 uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 6348 sqdist = SUM(uv(:)**2) 6349 uv = uv / SQRT(sqdist) 6350 6351 !-- reject raytracing above max distance 6352 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN 6353 ray_skip_maxdist = ray_skip_maxdist + 1 6354 CYCLE 6355 ENDIF 6356 6357 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 6358 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 6359 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 6360 / (pi * sqdist) & ! square of distance between centers 6361 * facearea(sd) 6362 6363 !-- reject raytracing for potentially too small view factor values 6364 IF ( rirrf < min_irrf_value ) THEN 6365 ray_skip_minval = ray_skip_minval + 1 6366 CYCLE 6367 ENDIF 6368 6369 !-- raytrace + process plant canopy sinks within 6370 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., & 6371 visible, transparency) 6372 6373 IF ( .NOT. visible ) CYCLE 6374 ! rsvf = rirrf * transparency 6375 6376 !-- write to the svf array 6377 nsvfl = nsvfl + 1 6378 !-- check dimmension of asvf array and enlarge it if needed 6379 IF ( nsvfla < nsvfl ) THEN 6380 k = CEILING(REAL(nsvfla, kind=wp) * grow_factor) 6381 IF ( msvf == 0 ) THEN 6382 msvf = 1 6383 ALLOCATE( asvf1(k) ) 6384 asvf => asvf1 6385 asvf1(1:nsvfla) = asvf2 6386 DEALLOCATE( asvf2 ) 6387 ELSE 6388 msvf = 0 6389 ALLOCATE( asvf2(k) ) 6390 asvf => asvf2 6391 asvf2(1:nsvfla) = asvf1 6392 DEALLOCATE( asvf1 ) 6393 ENDIF 6394 6395 WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 6396 CALL radiation_write_debug_log( msg ) 5921 6397 5922 6398 nsvfla = k … … 5931 6407 ENDDO 5932 6408 5933 !--Raytrace to canopy boxes to fill dsitransc TODO optimize 5934 !-- 5935 dsitransc(:,:) = -999._wp !FIXME6409 !-- 6410 !-- Raytrace to canopy boxes to fill dsitransc TODO optimize 6411 dsitransc(:,:) = 0._wp 5936 6412 az0 = 0._wp 5937 6413 naz = raytrace_discrete_azims … … 5940 6416 nzn = raytrace_discrete_elevs / 2 5941 6417 zns = pi / 2._wp / REAL(nzn, wp) 5942 ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn) ) 6418 ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn), & 6419 itarget(1:nzn) ) 5943 6420 zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/) 5944 6421 vffrac(:) = 0._wp 5945 6422 5946 DO ipcgb = 1, npcbl6423 DO ipcgb = 1, npcbl 5947 6424 ta = (/ REAL(pcbl(iz, ipcgb), wp), & 5948 6425 REAL(pcbl(iy, ipcgb), wp), & 5949 6426 REAL(pcbl(ix, ipcgb), wp) /) 5950 !--Calculate sky-view factor anddirect solar visibility using 2D raytracing5951 DO iaz = 1, naz6427 !-- Calculate direct solar visibility using 2D raytracing 6428 DO iaz = 1, naz 5952 6429 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 5953 6430 yxdir = (/ COS(azmid), SIN(azmid) /) 5954 CALL raytrace_2d(ta, yxdir, zdirs,&5955 -999, -999._wp, vffrac, .FALSE., . TRUE., &5956 win_lad, horizon, ztransp)5957 5958 !--Save direct solar transparency6431 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6432 -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., & 6433 lowest_free_ray, ztransp, itarget) !FIXME unit vect in grid units + zdirs 6434 6435 !-- Save direct solar transparency 5959 6436 j = MODULO(NINT(azmid/ & 5960 6437 (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), & 5961 6438 raytrace_discrete_azims) 5962 DO k = 1, raytrace_discrete_elevs/26439 DO k = 1, raytrace_discrete_elevs/2 5963 6440 i = dsidir_rev(k-1, j) 5964 IF ( i /= -1 ) dsitransc(ipcgb, i) = ztransp(k) 6441 IF ( i /= -1 .AND. k <= lowest_free_ray ) & 6442 dsitransc(ipcgb, i) = ztransp(k) 5965 6443 ENDDO 5966 6444 ENDDO 5967 6445 ENDDO 5968 DEALLOCATE ( zdirs, vffrac, ztransp ) 5969 5970 ! CALL radiation_write_debug_log( 'End of calculation SVF' ) 5971 ! WRITE(msg, *) 'Raytracing skipped for maximum distance of ', & 5972 ! max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.' 5973 ! CALL radiation_write_debug_log( msg ) 5974 ! WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', & 5975 ! min_irrf_value , ' on ', ray_skip_minval, ' pairs.' 5976 ! CALL radiation_write_debug_log( msg ) 6446 DEALLOCATE ( zdirs, vffrac, ztransp, itarget ) 6447 !-- 6448 !-- Raytrace to MRT boxes 6449 IF ( nmrtbl > 0 ) THEN 6450 mrtdsit(:,:) = 0._wp 6451 mrtsky(:) = 0._wp 6452 mrtskyt(:) = 0._wp 6453 az0 = 0._wp 6454 naz = raytrace_discrete_azims 6455 azs = 2._wp * pi / REAL(naz, wp) 6456 zn0 = 0._wp 6457 nzn = raytrace_discrete_elevs 6458 zns = pi / REAL(nzn, wp) 6459 ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), & 6460 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6461 !in case of rad_angular_discretization 6462 6463 zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/) 6464 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 6465 vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp) 6466 6467 DO imrt = 1, nmrtbl 6468 ta = (/ REAL(mrtbl(iz, imrt), wp), & 6469 REAL(mrtbl(iy, imrt), wp), & 6470 REAL(mrtbl(ix, imrt), wp) /) 6471 ! 6472 !-- vf fractions are constant per azimuth 6473 DO iaz = 0, naz-1 6474 vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac0(:) 6475 ENDDO 6476 !-- sum of whole vffrac equals 1, verified 6477 itarg0 = 1 6478 itarg1 = nzn 6479 ! 6480 !-- Calculate sky-view factor and direct solar visibility using 2D raytracing 6481 DO iaz = 1, naz 6482 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 6483 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), nzn, zdirs, & 6484 -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., & 6485 .FALSE., .TRUE., lowest_free_ray, & 6486 ztransp(itarg0:itarg1), & 6487 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6488 !FIXME itarget available only in 6489 !case of rad_angular_discretization 6490 6491 !-- Sky view factors for MRT 6492 mrtsky(imrt) = mrtsky(imrt) + & 6493 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) 6494 mrtskyt(imrt) = mrtskyt(imrt) + & 6495 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) & 6496 * vffrac(itarg0:itarg0+lowest_free_ray-1)) 6497 !-- Direct solar transparency for MRT 6498 j = MODULO(NINT(azmid/ & 6499 (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), & 6500 raytrace_discrete_azims) 6501 DO k = 1, raytrace_discrete_elevs/2 6502 i = dsidir_rev(k-1, j) 6503 IF ( i /= -1 .AND. k <= lowest_free_ray ) & 6504 mrtdsit(imrt, i) = ztransp(itarg0+k-1) 6505 ENDDO 6506 ! 6507 !-- Advance itarget indices 6508 itarg0 = itarg1 + 1 6509 itarg1 = itarg1 + nzn 6510 ENDDO 6511 6512 !-- sort itarget by face id 6513 CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz) 6514 ! 6515 !-- find the first valid position 6516 itarg0 = 1 6517 DO WHILE ( itarg0 <= nzn*naz ) 6518 IF ( itarget(itarg0) /= -1 ) EXIT 6519 itarg0 = itarg0 + 1 6520 ENDDO 6521 6522 DO i = itarg0, nzn*naz 6523 ! 6524 !-- For duplicate values, only sum up vf fraction value 6525 IF ( i < nzn*naz ) THEN 6526 IF ( itarget(i+1) == itarget(i) ) THEN 6527 vffrac(i+1) = vffrac(i+1) + vffrac(i) 6528 CYCLE 6529 ENDIF 6530 ENDIF 6531 ! 6532 !-- write to the mrtf array 6533 nmrtf = nmrtf + 1 6534 !-- check dimmension of mrtf array and enlarge it if needed 6535 IF ( nmrtfa < nmrtf ) THEN 6536 k = CEILING(REAL(nmrtfa, kind=wp) * grow_factor) 6537 IF ( mmrtf == 0 ) THEN 6538 mmrtf = 1 6539 ALLOCATE( amrtf1(k) ) 6540 amrtf => amrtf1 6541 amrtf1(1:nmrtfa) = amrtf2 6542 DEALLOCATE( amrtf2 ) 6543 ELSE 6544 mmrtf = 0 6545 ALLOCATE( amrtf2(k) ) 6546 amrtf => amrtf2 6547 amrtf2(1:nmrtfa) = amrtf1 6548 DEALLOCATE( amrtf1 ) 6549 ENDIF 6550 6551 WRITE (msg,'(A,3I12)') 'Grow amrtf:', nmrtf, nmrtfa, k 6552 CALL radiation_write_debug_log( msg ) 6553 6554 nmrtfa = k 6555 ENDIF 6556 !-- write mrtf values into the array 6557 amrtf(nmrtf)%isurflt = imrt 6558 amrtf(nmrtf)%isurfs = itarget(i) 6559 amrtf(nmrtf)%rsvf = vffrac(i) 6560 amrtf(nmrtf)%rtransp = ztransp(i) 6561 ENDDO ! itarg 6562 6563 ENDDO ! imrt 6564 DEALLOCATE ( zdirs, zbdry, vffrac, vffrac0, ztransp, itarget ) 6565 ! 6566 !-- Move MRT factors to final arrays 6567 ALLOCATE ( mrtf(nmrtf), mrtft(nmrtf), mrtfsurf(2,nmrtf) ) 6568 DO imrtf = 1, nmrtf 6569 mrtf(imrtf) = amrtf(imrtf)%rsvf 6570 mrtft(imrtf) = amrtf(imrtf)%rsvf * amrtf(imrtf)%rtransp 6571 mrtfsurf(:,imrtf) = (/amrtf(imrtf)%isurflt, amrtf(imrtf)%isurfs /) 6572 ENDDO 6573 IF ( ALLOCATED(amrtf1) ) DEALLOCATE( amrtf1 ) 6574 IF ( ALLOCATED(amrtf2) ) DEALLOCATE( amrtf2 ) 6575 ENDIF ! nmrtbl > 0 6576 6577 IF ( rad_angular_discretization ) THEN 6578 #if defined( __parallel ) 6579 !-- finalize MPI_RMA communication established to get global index of the surface from grid indices 6580 !-- flush all MPI window pending requests 6581 CALL MPI_Win_flush_all(win_gridsurf, ierr) 6582 IF ( ierr /= 0 ) THEN 6583 WRITE(9,*) 'Error MPI_Win_flush_all1:', ierr, win_gridsurf 6584 FLUSH(9) 6585 ENDIF 6586 !-- unlock MPI window 6587 CALL MPI_Win_unlock_all(win_gridsurf, ierr) 6588 IF ( ierr /= 0 ) THEN 6589 WRITE(9,*) 'Error MPI_Win_unlock_all1:', ierr, win_gridsurf 6590 FLUSH(9) 6591 ENDIF 6592 !-- free MPI window 6593 CALL MPI_Win_free(win_gridsurf, ierr) 6594 IF ( ierr /= 0 ) THEN 6595 WRITE(9,*) 'Error MPI_Win_free1:', ierr, win_gridsurf 6596 FLUSH(9) 6597 ENDIF 6598 #else 6599 DEALLOCATE ( gridsurf ) 6600 #endif 6601 ENDIF 6602 6603 CALL radiation_write_debug_log( 'End of calculation SVF' ) 6604 WRITE(msg, *) 'Raytracing skipped for maximum distance of ', & 6605 max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.' 6606 CALL radiation_write_debug_log( msg ) 6607 WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', & 6608 min_irrf_value , ' on ', ray_skip_minval, ' pairs.' 6609 CALL radiation_write_debug_log( msg ) 5977 6610 5978 6611 CALL location_message( ' waiting for completion of SVF and CSF calculation in all processes', .TRUE. ) … … 5983 6616 !-- finalize mpi_rma communication and deallocate temporary arrays 5984 6617 #if defined( __parallel ) 5985 IF ( r ma_lad_raytrace) THEN6618 IF ( raytrace_mpi_rma ) THEN 5986 6619 CALL MPI_Win_flush_all(win_lad, ierr) 6620 IF ( ierr /= 0 ) THEN 6621 WRITE(9,*) 'Error MPI_Win_flush_all2:', ierr, win_lad 6622 FLUSH(9) 6623 ENDIF 5987 6624 !-- unlock MPI window 5988 6625 CALL MPI_Win_unlock_all(win_lad, ierr) 6626 IF ( ierr /= 0 ) THEN 6627 WRITE(9,*) 'Error MPI_Win_unlock_all2:', ierr, win_lad 6628 FLUSH(9) 6629 ENDIF 5989 6630 !-- free MPI window 5990 6631 CALL MPI_Win_free(win_lad, ierr) 5991 6632 IF ( ierr /= 0 ) THEN 6633 WRITE(9,*) 'Error MPI_Win_free2:', ierr, win_lad 6634 FLUSH(9) 6635 ENDIF 5992 6636 !-- deallocate temporary arrays storing values for csf calculation during raytracing 5993 6637 DEALLOCATE( lad_s_ray ) 5994 !-- sub_lad is the pointer to lad_s_rma in case of r ma_lad_raytrace6638 !-- sub_lad is the pointer to lad_s_rma in case of raytrace_mpi_rma 5995 6639 !-- and must not be deallocated here 5996 6640 ELSE … … 6009 6653 CALL location_message( ' calculation of the complete SVF array', .TRUE. ) 6010 6654 6011 ! CALL radiation_write_debug_log( 'Start SVF sort' ) 6012 !-- sort svf ( a version of quicksort ) 6013 CALL quicksort_svf(asvf,1,nsvfl) 6014 6015 !< load svf from the structure array to plain arrays 6016 ! CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' ) 6017 ALLOCATE( svf(ndsvf,nsvfl) ) 6018 ALLOCATE( svfsurf(idsvf,nsvfl) ) 6019 svfnorm_counts(:) = 0._wp 6020 isurflt_prev = -1 6021 ksvf = 1 6022 svfsum = 0._wp 6023 DO isvf = 1, nsvfl 6024 !-- normalize svf per target face 6025 IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN 6026 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 6027 !< update histogram of logged svf normalization values 6028 i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev))) 6029 svfnorm_counts(i) = svfnorm_counts(i) + 1 6030 6031 svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev)) 6032 ENDIF 6033 isurflt_prev = asvf(ksvf)%isurflt 6034 isvf_surflt = isvf 6035 svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 6036 ELSE 6037 svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 6038 ENDIF 6039 6040 svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /) 6041 svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /) 6042 6043 !-- next element 6044 ksvf = ksvf + 1 6045 ENDDO 6046 6047 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 6048 i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev))) 6049 svfnorm_counts(i) = svfnorm_counts(i) + 1 6050 6051 svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev)) 6052 ENDIF 6053 !TODO we should be able to deallocate skyvf, from now on we only need skyvft 6655 IF ( rad_angular_discretization ) THEN 6656 CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' ) 6657 ALLOCATE( svf(ndsvf,nsvfl) ) 6658 ALLOCATE( svfsurf(idsvf,nsvfl) ) 6659 6660 DO isvf = 1, nsvfl 6661 svf(:, isvf) = (/ asvf(isvf)%rsvf, asvf(isvf)%rtransp /) 6662 svfsurf(:, isvf) = (/ asvf(isvf)%isurflt, asvf(isvf)%isurfs /) 6663 ENDDO 6664 ELSE 6665 CALL radiation_write_debug_log( 'Start SVF sort' ) 6666 !-- sort svf ( a version of quicksort ) 6667 CALL quicksort_svf(asvf,1,nsvfl) 6668 6669 !< load svf from the structure array to plain arrays 6670 CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' ) 6671 ALLOCATE( svf(ndsvf,nsvfl) ) 6672 ALLOCATE( svfsurf(idsvf,nsvfl) ) 6673 svfnorm_counts(:) = 0._wp 6674 isurflt_prev = -1 6675 ksvf = 1 6676 svfsum = 0._wp 6677 DO isvf = 1, nsvfl 6678 !-- normalize svf per target face 6679 IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN 6680 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 6681 !< update histogram of logged svf normalization values 6682 i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev))) 6683 svfnorm_counts(i) = svfnorm_counts(i) + 1 6684 6685 svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev)) 6686 ENDIF 6687 isurflt_prev = asvf(ksvf)%isurflt 6688 isvf_surflt = isvf 6689 svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 6690 ELSE 6691 svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp 6692 ENDIF 6693 6694 svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /) 6695 svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /) 6696 6697 !-- next element 6698 ksvf = ksvf + 1 6699 ENDDO 6700 6701 IF ( isurflt_prev /= -1 .AND. svfsum /= 0._wp ) THEN 6702 i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev))) 6703 svfnorm_counts(i) = svfnorm_counts(i) + 1 6704 6705 svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev)) 6706 ENDIF 6707 WRITE(9, *) 'SVF normalization histogram:', svfnorm_counts, & 6708 'on thresholds:', svfnorm_report_thresh(1:svfnorm_report_num), '(val < thresh <= val)' 6709 !TODO we should be able to deallocate skyvf, from now on we only need skyvft 6710 ENDIF ! rad_angular_discretization 6054 6711 6055 6712 !-- deallocate temporary asvf array … … 6067 6724 6068 6725 CALL location_message( ' calculation of the complete CSF array', .TRUE. ) 6069 !CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )6726 CALL radiation_write_debug_log( 'Calculation of the complete CSF array' ) 6070 6727 !-- sort and merge csf for the last time, keeping the array size to minimum 6071 6728 CALL merge_and_grow_csf(-1) … … 6073 6730 !-- aggregate csb among processors 6074 6731 !-- allocate necessary arrays 6075 ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) ) 6076 ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) ) 6732 udim = max(ncsfl,1) 6733 ALLOCATE( csflt_l(ndcsf*udim) ) 6734 csflt(1:ndcsf,1:udim) => csflt_l(1:ndcsf*udim) 6735 ALLOCATE( kcsflt_l(kdcsf*udim) ) 6736 kcsflt(1:kdcsf,1:udim) => kcsflt_l(1:kdcsf*udim) 6077 6737 ALLOCATE( icsflt(0:numprocs-1) ) 6078 6738 ALLOCATE( dcsflt(0:numprocs-1) ) … … 6108 6768 !-- fill out real values of rsvf, rtransp 6109 6769 csflt(1,kcsf) = acsf(kcsf)%rsvf 6110 csflt(2,kcsf) = acsf(kcsf)%rtransp6111 6770 !-- fill out integer values of itz,ity,itx,isurfs 6112 6771 kcsflt(1,kcsf) = acsf(kcsf)%itz … … 6138 6797 !-- scatter and gather the number of elements to and from all processor 6139 6798 !-- and calculate displacements 6140 !CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )6799 CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' ) 6141 6800 CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr) 6142 6801 IF ( ierr /= 0 ) THEN 6802 WRITE(9,*) 'Error MPI_AlltoAll1:', ierr, SIZE(icsflt), SIZE(ipcsflt) 6803 FLUSH(9) 6804 ENDIF 6805 6143 6806 npcsfl = SUM(ipcsflt) 6144 6807 d = 0 … … 6149 6812 6150 6813 !-- exchange csf fields between processors 6151 ! CALL radiation_write_debug_log( 'Exchange csf fields between processors' ) 6152 ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) ) 6153 ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) ) 6154 CALL MPI_AlltoAllv(csflt, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, & 6155 pcsflt, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr) 6156 CALL MPI_AlltoAllv(kcsflt, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, & 6157 kpcsflt, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr) 6814 CALL radiation_write_debug_log( 'Exchange csf fields between processors' ) 6815 udim = max(npcsfl,1) 6816 ALLOCATE( pcsflt_l(ndcsf*udim) ) 6817 pcsflt(1:ndcsf,1:udim) => pcsflt_l(1:ndcsf*udim) 6818 ALLOCATE( kpcsflt_l(kdcsf*udim) ) 6819 kpcsflt(1:kdcsf,1:udim) => kpcsflt_l(1:kdcsf*udim) 6820 CALL MPI_AlltoAllv(csflt_l, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, & 6821 pcsflt_l, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr) 6822 IF ( ierr /= 0 ) THEN 6823 WRITE(9,*) 'Error MPI_AlltoAllv1:', ierr, SIZE(ipcsflt), ndcsf*icsflt, & 6824 ndcsf*dcsflt, SIZE(pcsflt_l),ndcsf*ipcsflt, ndcsf*dpcsflt 6825 FLUSH(9) 6826 ENDIF 6827 6828 CALL MPI_AlltoAllv(kcsflt_l, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, & 6829 kpcsflt_l, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr) 6830 IF ( ierr /= 0 ) THEN 6831 WRITE(9,*) 'Error MPI_AlltoAllv2:', ierr, SIZE(kcsflt_l),kdcsf*icsflt, & 6832 kdcsf*dcsflt, SIZE(kpcsflt_l), kdcsf*ipcsflt, kdcsf*dpcsflt 6833 FLUSH(9) 6834 ENDIF 6158 6835 6159 6836 #else … … 6166 6843 6167 6844 !-- deallocate temporary arrays 6168 DEALLOCATE( csflt )6169 DEALLOCATE( kcsflt )6845 DEALLOCATE( csflt_l ) 6846 DEALLOCATE( kcsflt_l ) 6170 6847 DEALLOCATE( icsflt ) 6171 6848 DEALLOCATE( dcsflt ) … … 6174 6851 6175 6852 !-- sort csf ( a version of quicksort ) 6176 !CALL radiation_write_debug_log( 'Sort csf' )6853 CALL radiation_write_debug_log( 'Sort csf' ) 6177 6854 CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl) 6178 6855 6179 6856 !-- aggregate canopy sink factor records with identical box & source 6180 6857 !-- againg across all values from all processors 6181 !CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )6858 CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' ) 6182 6859 6183 6860 IF ( npcsfl > 0 ) THEN … … 6190 6867 kpcsflt(1,icsf) == kpcsflt(1,icsf+1) .AND. & 6191 6868 kpcsflt(4,icsf) == kpcsflt(4,icsf+1) ) THEN 6192 !-- We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray 6193 !-- probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that 6194 !-- might mean that the traced beam passes longer through the canopy box. 6195 IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) ) THEN 6196 pcsflt(2,kcsf) = pcsflt(2,icsf+1) 6197 ENDIF 6869 6198 6870 pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1) 6199 6871 … … 6224 6896 6225 6897 !-- deallocation of temporary arrays 6226 DEALLOCATE( pcsflt ) 6227 DEALLOCATE( kpcsflt ) 6228 ! CALL radiation_write_debug_log( 'End of aggregate csf' ) 6898 IF ( npcbl > 0 ) DEALLOCATE( gridpcbl ) 6899 DEALLOCATE( pcsflt_l ) 6900 DEALLOCATE( kpcsflt_l ) 6901 CALL radiation_write_debug_log( 'End of aggregate csf' ) 6229 6902 6230 6903 ENDIF … … 6233 6906 CALL MPI_BARRIER( comm2d, ierr ) 6234 6907 #endif 6235 !CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )6908 CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' ) 6236 6909 6237 6910 RETURN … … 6261 6934 !> doesn't need to be the same in all three dimensions). 6262 6935 !------------------------------------------------------------------------------! 6263 SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency , win_lad)6936 SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency) 6264 6937 IMPLICIT NONE 6265 6938 … … 6271 6944 LOGICAL, INTENT(out) :: visible 6272 6945 REAL(wp), INTENT(out) :: transparency !< along whole path 6273 INTEGER(iwp), INTENT(in) :: win_lad6274 6946 INTEGER(iwp) :: i, k, d 6275 6947 INTEGER(iwp) :: seldim !< dimension to be incremented … … 6294 6966 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array 6295 6967 REAL(wp) :: lad_s_target !< recieved lad_s of particular grid box 6296 REAL(wp), PARAMETER :: grow_factor = 1.5_wp !< factor of expansion of grow arrays6297 6968 6298 6969 ! … … 6385 7056 IF ( plant_canopy ) THEN 6386 7057 #if defined( __parallel ) 6387 IF ( r ma_lad_raytrace) THEN7058 IF ( raytrace_mpi_rma ) THEN 6388 7059 !-- send requests for lad_s to appropriate processor 6389 CALL cpu_log( log_point_s(77), 'rad_ init_rma', 'start' )7060 CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'start' ) 6390 7061 DO i = 1, ncsb 6391 7062 CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), & 6392 7063 1, MPI_REAL, win_lad, ierr) 6393 7064 IF ( ierr /= 0 ) THEN 6394 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get' 6395 CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 ) 7065 WRITE(9,*) 'Error MPI_Get1:', ierr, lad_s_ray(i), & 7066 lad_ip(i), lad_disp(i), win_lad 7067 FLUSH(9) 6396 7068 ENDIF 6397 7069 ENDDO … … 6400 7072 CALL MPI_Win_flush_local_all(win_lad, ierr) 6401 7073 IF ( ierr /= 0 ) THEN 6402 WRITE( message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'6403 CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0)7074 WRITE(9,*) 'Error MPI_Win_flush_local_all1:', ierr, win_lad 7075 FLUSH(9) 6404 7076 ENDIF 6405 CALL cpu_log( log_point_s(77), 'rad_ init_rma', 'stop' )7077 CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'stop' ) 6406 7078 6407 7079 ENDIF … … 6411 7083 DO i = 1, ncsb 6412 7084 #if defined( __parallel ) 6413 IF ( r ma_lad_raytrace) THEN7085 IF ( raytrace_mpi_rma ) THEN 6414 7086 lad_s_target = lad_s_ray(i) 6415 7087 ELSE … … 6429 7101 acsf(ncsfl)%itz = boxes(1,i) 6430 7102 acsf(ncsfl)%isurfs = isrc 6431 acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency 6432 acsf(ncsfl)%rtransp = REAL(transparency, wp) 7103 acsf(ncsfl)%rsvf = cursink*transparency*rirrf*atarg 6433 7104 ENDIF !< create_csf 6434 7105 … … 6452 7123 !> vertical_delta / horizontal_distance 6453 7124 !------------------------------------------------------------------------------! 6454 SUBROUTINE raytrace_2d(origin, yxdir, zdirs, iorig, aorig, vffrac,&6455 c reate_csf, skip_1st_pcb, win_lad, horizon,&6456 transparency)7125 SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac, & 7126 calc_svf, create_csf, skip_1st_pcb, & 7127 lowest_free_ray, transparency, itarget) 6457 7128 IMPLICIT NONE 6458 7129 6459 7130 REAL(wp), DIMENSION(3), INTENT(IN) :: origin !< z,y,x coordinates of ray origin 6460 7131 REAL(wp), DIMENSION(2), INTENT(IN) :: yxdir !< y,x *unit* vector of ray direction (in grid units) 6461 REAL(wp), DIMENSION(:), INTENT(IN) :: zdirs !< list of z directions to raytrace (z/hdist, in grid) 7132 INTEGER(iwp) :: nrays !< number of rays (z directions) to raytrace 7133 REAL(wp), DIMENSION(nrays), INTENT(IN) :: zdirs !< list of z directions to raytrace (z/hdist, grid, zenith->nadir) 6462 7134 INTEGER(iwp), INTENT(in) :: iorig !< index of origin face for csf 6463 7135 REAL(wp), INTENT(in) :: aorig !< origin face area for csf 6464 REAL(wp), DIMENSION( LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(in) :: vffrac !<6465 !< view factor fractions of each ray for csf6466 LOGICAL, INTENT(in) :: create_csf !< whether to generate new CSFs during raytracing7136 REAL(wp), DIMENSION(nrays), INTENT(in) :: vffrac !< view factor fractions of each ray for csf 7137 LOGICAL, INTENT(in) :: calc_svf !< whether to calculate SFV (identify obstacle surfaces) 7138 LOGICAL, INTENT(in) :: create_csf !< whether to create canopy sink factors 6467 7139 LOGICAL, INTENT(in) :: skip_1st_pcb !< whether to skip first plant canopy box during raytracing 6468 INTEGER(iwp), INTENT( in) :: win_lad !< leaf area density MPI window6469 REAL(wp), INTENT(OUT) :: horizon !< highest horizon found after raytracing (z/hdist)6470 REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) :: transparency !<6471 !< transparencies of zdirs paths 6472 !--INTEGER(iwp), DIMENSION(3, LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) :: itarget !<6473 !< (z,y,x) coordinates of target faces for zdirs7140 INTEGER(iwp), INTENT(out) :: lowest_free_ray !< index into zdirs 7141 REAL(wp), DIMENSION(nrays), INTENT(OUT) :: transparency !< transparencies of zdirs paths 7142 INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) :: itarget !< global indices of target faces for zdirs 7143 7144 INTEGER(iwp), DIMENSION(nrays) :: target_procs 7145 REAL(wp) :: horizon !< highest horizon found after raytracing (z/hdist) 6474 7146 INTEGER(iwp) :: i, k, l, d 6475 7147 INTEGER(iwp) :: seldim !< dimension to be incremented … … 6506 7178 INTEGER(iwp) :: iz 6507 7179 INTEGER(iwp) :: zsgn 6508 REAL(wp), PARAMETER :: grow_factor = 1.5_wp !< factor of expansion of grow arrays 7180 INTEGER(iwp) :: lowest_lad !< lowest column cell for which we need LAD 7181 INTEGER(iwp) :: lastdir !< wall direction before hitting this column 7182 INTEGER(iwp), DIMENSION(2) :: lastcolumn 6509 7183 6510 7184 #if defined( __parallel ) … … 6515 7189 transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing 6516 7190 horizon = -HUGE(1._wp) 6517 6518 !--Determine distance to boundary (in 2D xy) 7191 lowest_free_ray = nrays 7192 IF ( rad_angular_discretization .AND. calc_svf ) THEN 7193 ALLOCATE(target_surfl(nrays)) 7194 target_surfl(:) = -1 7195 lastdir = -999 7196 lastcolumn(:) = -999 7197 ENDIF 7198 7199 !-- Determine distance to boundary (in 2D xy) 6519 7200 IF ( yxdir(1) > 0._wp ) THEN 6520 7201 bdydim = ny + .5_wp !< north global boundary … … 6548 7229 !-- Since all face coordinates have values *.5 and we'd like to use 6549 7230 !-- integers, all these have .5 added 6550 DO d = 1, 27231 DO d = 1, 2 6551 7232 IF ( yxdir(d) == 0._wp ) THEN 6552 7233 dimnext(d) = HUGE(1_iwp) … … 6569 7250 seldim = minloc(dimnextdist, 1) 6570 7251 nextdist = dimnextdist(seldim) 6571 IF ( nextdist > distance ) nextdist = distance7252 IF ( nextdist > distance ) nextdist = distance 6572 7253 6573 7254 IF ( nextdist > lastdist ) THEN 6574 7255 ntrack = ntrack + 1 6575 7256 crmid = (lastdist + nextdist) * .5_wp 6576 column = INT(yxorigin(:) + yxdir(:) * crmid, iwp) !NINT(yxorigin(:) + yxdir(:) * crmid, iwp)7257 column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp) 6577 7258 6578 7259 !-- calculate index of the grid with global indices (column(1),column(2)) … … 6586 7267 horz_entry = -HUGE(1._wp) 6587 7268 ELSE 6588 horz_entry = ( nzterr(ig)- origin(1)) / lastdist7269 horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist 6589 7270 ENDIF 6590 horz_exit = (nzterr(ig) - origin(1)) / nextdist 7271 horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist 7272 7273 IF ( rad_angular_discretization .AND. calc_svf ) THEN 7274 ! 7275 !-- Identify vertical obstacles hit by rays in current column 7276 DO WHILE ( lowest_free_ray > 0 ) 7277 IF ( zdirs(lowest_free_ray) > horz_entry ) EXIT 7278 ! 7279 !-- This may only happen after 1st column, so lastdir and lastcolumn are valid 7280 CALL request_itarget(lastdir, & 7281 CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), & 7282 lastcolumn(1), lastcolumn(2), & 7283 target_surfl(lowest_free_ray), target_procs(lowest_free_ray)) 7284 lowest_free_ray = lowest_free_ray - 1 7285 ENDDO 7286 ! 7287 !-- Identify horizontal obstacles hit by rays in current column 7288 DO WHILE ( lowest_free_ray > 0 ) 7289 IF ( zdirs(lowest_free_ray) > horz_exit ) EXIT 7290 CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), & 7291 target_surfl(lowest_free_ray), & 7292 target_procs(lowest_free_ray)) 7293 lowest_free_ray = lowest_free_ray - 1 7294 ENDDO 7295 ENDIF 7296 6591 7297 horizon = MAX(horizon, horz_entry, horz_exit) 6592 7298 … … 6598 7304 6599 7305 IF ( nextdist >= distance ) EXIT 7306 7307 IF ( rad_angular_discretization .AND. calc_svf ) THEN 7308 ! 7309 !-- Save wall direction of coming building column (= this air column) 7310 IF ( seldim == 1 ) THEN 7311 IF ( dimdelta(seldim) == 1 ) THEN 7312 lastdir = isouth_u 7313 ELSE 7314 lastdir = inorth_u 7315 ENDIF 7316 ELSE 7317 IF ( dimdelta(seldim) == 1 ) THEN 7318 lastdir = iwest_u 7319 ELSE 7320 lastdir = ieast_u 7321 ENDIF 7322 ENDIF 7323 lastcolumn = column 7324 ENDIF 6600 7325 lastdist = nextdist 6601 7326 dimnext(seldim) = dimnext(seldim) + dimdelta(seldim) … … 6604 7329 6605 7330 IF ( plant_canopy ) THEN 6606 !--Request LAD WHERE applicable6607 !-- 7331 !-- Request LAD WHERE applicable 7332 !-- 6608 7333 #if defined( __parallel ) 6609 IF ( r ma_lad_raytrace) THEN7334 IF ( raytrace_mpi_rma ) THEN 6610 7335 !-- send requests for lad_s to appropriate processor 6611 7336 !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' ) 6612 DO i = 1, ntrack7337 DO i = 1, ntrack 6613 7338 px = rt2_track(2,i)/nnx 6614 7339 py = rt2_track(1,i)/nny 6615 7340 ip = px*pdims(2)+py 6616 7341 ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny 6617 IF ( plantt(ig) <= nzterr(ig) ) CYCLE 6618 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + nzterr(ig)+1-nzub 6619 wcount = plantt(ig)-nzterr(ig) 7342 7343 IF ( rad_angular_discretization .AND. calc_svf ) THEN 7344 ! 7345 !-- For fixed view resolution, we need plant canopy even for rays 7346 !-- to opposing surfaces 7347 lowest_lad = nzterr(ig) + 1 7348 ELSE 7349 ! 7350 !-- We only need LAD for rays directed above horizon (to sky) 7351 lowest_lad = CEILING( -0.5_wp + origin(1) + & 7352 MIN( horizon * rt2_track_dist(i-1), & ! entry 7353 horizon * rt2_track_dist(i) ) ) ! exit 7354 ENDIF 7355 ! 7356 !-- Skip asking for LAD where all plant canopy is under requested level 7357 IF ( plantt(ig) < lowest_lad ) CYCLE 7358 7359 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + lowest_lad-nzub 7360 wcount = plantt(ig)-lowest_lad+1 6620 7361 ! TODO send request ASAP - even during raytracing 6621 CALL MPI_Get(rt2_track_lad( nzterr(ig)+1:plantt(ig), i), wcount, MPI_REAL, ip, &7362 CALL MPI_Get(rt2_track_lad(lowest_lad:plantt(ig), i), wcount, MPI_REAL, ip, & 6622 7363 wdisp, wcount, MPI_REAL, win_lad, ierr) 6623 7364 IF ( ierr /= 0 ) THEN 6624 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get' 6625 CALL message( 'raytrace_2d', 'PA0526', 1, 2, 0, 6, 0 ) 7365 WRITE(9,*) 'Error MPI_Get2:', ierr, rt2_track_lad(lowest_lad:plantt(ig), i), & 7366 wcount, ip, wdisp, win_lad 7367 FLUSH(9) 6626 7368 ENDIF 6627 7369 ENDDO … … 6631 7373 CALL MPI_Win_flush_local_all(win_lad, ierr) 6632 7374 IF ( ierr /= 0 ) THEN 6633 WRITE( message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'6634 CALL message( 'raytrace', 'PA0527', 1, 2, 0, 6, 0)7375 WRITE(9,*) 'Error MPI_Win_flush_local_all2:', ierr, win_lad 7376 FLUSH(9) 6635 7377 ENDIF 6636 7378 !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' ) 6637 ELSE ! rma_lad_raytrace 6638 DO i = 1, ntrack 7379 7380 ELSE ! raytrace_mpi_rma = .F. 7381 DO i = 1, ntrack 6639 7382 px = rt2_track(2,i)/nnx 6640 7383 py = rt2_track(1,i)/nny … … 6645 7388 ENDIF 6646 7389 #else 6647 DO i = 1, ntrack7390 DO i = 1, ntrack 6648 7391 rt2_track_lad(nzub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max) 6649 7392 ENDDO 6650 7393 #endif 6651 6652 !--Skip the PCB around origin if requested 6653 !-- 6654 IF ( skip_1st_pcb ) THEN 7394 ENDIF ! plant_canopy 7395 7396 IF ( rad_angular_discretization .AND. calc_svf ) THEN 7397 #if defined( __parallel ) 7398 !-- wait for all gridsurf requests to complete 7399 CALL MPI_Win_flush_local_all(win_gridsurf, ierr) 7400 IF ( ierr /= 0 ) THEN 7401 WRITE(9,*) 'Error MPI_Win_flush_local_all3:', ierr, win_gridsurf 7402 FLUSH(9) 7403 ENDIF 7404 #endif 7405 ! 7406 !-- recalculate local surf indices into global ones 7407 DO i = 1, nrays 7408 IF ( target_surfl(i) == -1 ) THEN 7409 itarget(i) = -1 7410 ELSE 7411 itarget(i) = target_surfl(i) + surfstart(target_procs(i)) 7412 ENDIF 7413 ENDDO 7414 7415 DEALLOCATE( target_surfl ) 7416 7417 ELSE 7418 itarget(:) = -1 7419 ENDIF ! rad_angular_discretization 7420 7421 IF ( plant_canopy ) THEN 7422 !-- Skip the PCB around origin if requested (for MRT, the PCB might not be there) 7423 !-- 7424 IF ( skip_1st_pcb .AND. NINT(origin(1)) <= plantt_max ) THEN 6655 7425 rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp 6656 7426 ENDIF 6657 7427 6658 !--Assert that we have space allocated for CSFs 6659 !-- 6660 maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1) 7428 !-- Assert that we have space allocated for CSFs 7429 !-- 7430 maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub, & 7431 nzpt - CEILING(origin(1)-.5_wp))) * nrays 6661 7432 IF ( ncsfl + maxboxes > ncsfla ) THEN 6662 7433 !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) … … 6668 7439 ENDIF 6669 7440 6670 !--Calculate transparencies and store new CSFs6671 !-- 7441 !-- Calculate transparencies and store new CSFs 7442 !-- 6672 7443 zbottom = REAL(nzub, wp) - .5_wp 6673 7444 ztop = REAL(plantt_max, wp) + .5_wp 6674 7445 6675 !--Reverse direction of radiation (face->sky), only when create_csf6676 !-- 6677 IF ( c reate_csf ) THEN6678 DO i = 1, ntrack ! for each column7446 !-- Reverse direction of radiation (face->sky), only when calc_svf 7447 !-- 7448 IF ( calc_svf ) THEN 7449 DO i = 1, ntrack ! for each column 6679 7450 dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2 6680 7451 px = rt2_track(2,i)/nnx … … 6682 7453 ip = px*pdims(2)+py 6683 7454 6684 DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray 6685 IF ( zdirs(k) <= horizon ) THEN 6686 CYCLE 7455 DO k = 1, nrays ! for each ray 7456 ! 7457 !-- NOTE 6778: 7458 !-- With traditional svf discretization, CSFs under the horizon 7459 !-- (i.e. for surface to surface radiation) are created in 7460 !-- raytrace(). With rad_angular_discretization, we must create 7461 !-- CSFs under horizon only for one direction, otherwise we would 7462 !-- have duplicate amount of energy. Although we could choose 7463 !-- either of the two directions (they differ only by 7464 !-- discretization error with no bias), we choose the the backward 7465 !-- direction, because it tends to cumulate high canopy sink 7466 !-- factors closer to raytrace origin, i.e. it should potentially 7467 !-- cause less moiree. 7468 IF ( .NOT. rad_angular_discretization ) THEN 7469 IF ( zdirs(k) <= horizon ) CYCLE 6687 7470 ENDIF 6688 7471 6689 zorig = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)6690 IF ( zorig <= zbottom .OR.zorig >= ztop ) CYCLE7472 zorig = origin(1) + zdirs(k) * rt2_track_dist(i-1) 7473 IF ( zorig <= zbottom .OR. zorig >= ztop ) CYCLE 6691 7474 6692 7475 zsgn = INT(SIGN(1._wp, zdirs(k)), iwp) … … 6695 7478 nz = 2 6696 7479 rt2_dist(nz) = SQRT(dxxyy) 6697 iz = NINT(zorig, iwp)7480 iz = CEILING(-.5_wp + zorig, iwp) 6698 7481 ELSE 6699 zexit = MIN(MAX( REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)7482 zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop) 6700 7483 6701 7484 zb0 = FLOOR( zorig * zsgn - .5_wp) + 1 ! because it must be greater than orig … … 6708 7491 ENDIF 6709 7492 6710 DO l = 2, nz7493 DO l = 2, nz 6711 7494 IF ( rt2_track_lad(iz, i) > 0._wp ) THEN 6712 7495 curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1))) 6713 ncsfl = ncsfl + 1 6714 acsf(ncsfl)%ip = ip 6715 acsf(ncsfl)%itx = rt2_track(2,i) 6716 acsf(ncsfl)%ity = rt2_track(1,i) 6717 acsf(ncsfl)%itz = iz 6718 acsf(ncsfl)%isurfs = iorig 6719 acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency 6720 acsf(ncsfl)%rtransp = REAL(transparency(k), wp) 7496 7497 IF ( create_csf ) THEN 7498 ncsfl = ncsfl + 1 7499 acsf(ncsfl)%ip = ip 7500 acsf(ncsfl)%itx = rt2_track(2,i) 7501 acsf(ncsfl)%ity = rt2_track(1,i) 7502 acsf(ncsfl)%itz = iz 7503 acsf(ncsfl)%isurfs = iorig 7504 acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k) 7505 ENDIF 6721 7506 6722 7507 transparency(k) = transparency(k) * curtrans … … 6724 7509 iz = iz + zsgn 6725 7510 ENDDO ! l = 1, nz - 1 6726 ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)7511 ENDDO ! k = 1, nrays 6727 7512 ENDDO ! i = 1, ntrack 6728 7513 6729 transparency( :) = 1._wp !-- Reset all rays to transparent7514 transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778) 6730 7515 ENDIF 6731 7516 6732 !--Forward direction of radiation (sky->face), always6733 !-- 6734 DO i = ntrack, 1, -1 ! for each column backwards7517 !-- Forward direction of radiation (sky->face), always 7518 !-- 7519 DO i = ntrack, 1, -1 ! for each column backwards 6735 7520 dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2 6736 7521 px = rt2_track(2,i)/nnx … … 6738 7523 ip = px*pdims(2)+py 6739 7524 6740 DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray 6741 IF ( zdirs(k) <= horizon ) THEN 6742 transparency(k) = 0._wp 6743 CYCLE 6744 ENDIF 6745 6746 zexit = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1) 6747 IF ( zexit <= zbottom .OR. zexit >= ztop ) CYCLE 7525 DO k = 1, nrays ! for each ray 7526 ! 7527 !-- See NOTE 6778 above 7528 IF ( zdirs(k) <= horizon ) CYCLE 7529 7530 zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1) 7531 IF ( zexit <= zbottom .OR. zexit >= ztop ) CYCLE 6748 7532 6749 7533 zsgn = -INT(SIGN(1._wp, zdirs(k)), iwp) … … 6754 7538 iz = NINT(zexit, iwp) 6755 7539 ELSE 6756 zorig = MIN(MAX( REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)7540 zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop) 6757 7541 6758 7542 zb0 = FLOOR( zorig * zsgn - .5_wp) + 1 ! because it must be greater than orig … … 6765 7549 ENDIF 6766 7550 6767 DO l = 2, nz7551 DO l = 2, nz 6768 7552 IF ( rt2_track_lad(iz, i) > 0._wp ) THEN 6769 7553 curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1))) … … 6775 7559 acsf(ncsfl)%ity = rt2_track(1,i) 6776 7560 acsf(ncsfl)%itz = iz 6777 acsf(ncsfl)%isurfs = -1 ! a special ID indicating sky6778 acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency6779 acsf(ncsfl)%r transp = REAL(transparency(k), wp)6780 ENDIF ! <create_csf7561 acsf(ncsfl)%isurfs = itarget(k) !if above horizon, then itarget(k)==-1, which 7562 !is also a special ID indicating sky 7563 acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k) 7564 ENDIF ! create_csf 6781 7565 6782 7566 transparency(k) = transparency(k) * curtrans … … 6784 7568 iz = iz + zsgn 6785 7569 ENDDO ! l = 1, nz - 1 6786 ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)7570 ENDDO ! k = 1, nrays 6787 7571 ENDDO ! i = 1, ntrack 6788 6789 ELSE ! not plant_canopy 6790 DO k = UBOUND(zdirs, 1), LBOUND(zdirs, 1), -1 ! TODO make more generic 6791 IF ( zdirs(k) > horizon ) EXIT 6792 transparency(k) = 0._wp 7572 ENDIF ! plant_canopy 7573 7574 IF ( .NOT. (rad_angular_discretization .AND. calc_svf) ) THEN 7575 ! 7576 !-- Just update lowest_free_ray according to horizon 7577 DO WHILE ( lowest_free_ray > 0 ) 7578 IF ( zdirs(lowest_free_ray) > horizon ) EXIT 7579 lowest_free_ray = lowest_free_ray - 1 6793 7580 ENDDO 6794 7581 ENDIF 7582 7583 CONTAINS 7584 SUBROUTINE request_itarget(d, z, y, x, isurfl, iproc) 7585 INTEGER(iwp), INTENT(in) :: d, z, y, x 7586 INTEGER(iwp), TARGET, INTENT(out) :: isurfl 7587 INTEGER(iwp), INTENT(out) :: iproc 7588 INTEGER(kind=MPI_ADDRESS_KIND) :: target_displ !< index of the grid in the local gridsurf array 7589 INTEGER(iwp) :: px, py !< number of processors in x and y direction 7590 !< before the processor in the question 7591 7592 !-- calculate target processor and index in the remote local target gridsurf array 7593 px = x/nnx 7594 py = y/nny 7595 iproc = px*pdims(2)+py 7596 target_displ = ((x-px*nnx)*nny + y-py*nny)*nzu*nsurf_type_u + (z-nzub)*nsurf_type_u + d 7597 7598 #if defined( __parallel ) 7599 !-- send MPI_Get request to obtain index target_surfl(i) 7600 CALL MPI_Get(isurfl, 1, MPI_INTEGER, iproc, target_displ, & 7601 1, MPI_INTEGER, win_gridsurf, ierr) 7602 IF ( ierr /= 0 ) THEN 7603 WRITE(9,*) 'Error MPI_Get3:', ierr, isurfl, iproc, target_displ, win_gridsurf 7604 FLUSH(9) 7605 ENDIF 7606 #else 7607 !-- set index target_surfl(i) 7608 isurfl = gridsurf(d,z,y,x) 7609 #endif 7610 END SUBROUTINE request_itarget 6795 7611 6796 7612 END SUBROUTINE raytrace_2d … … 6849 7665 ALLOCATE ( dsitrans(nsurfl, ndsidir) ) 6850 7666 ALLOCATE ( dsitransc(npcbl, ndsidir) ) 7667 IF ( nmrtbl > 0 ) ALLOCATE ( mrtdsit(nmrtbl, ndsidir) ) 6851 7668 6852 7669 WRITE ( message_string, * ) 'Precalculated', ndsidir, ' solar positions', & … … 6905 7722 6906 7723 !-- first check: are the two surfaces directed in the same direction 6907 IF ( (d==iup_u .OR. d==iup_l .OR. d==iup_a) &7724 IF ( (d==iup_u .OR. d==iup_l ) & 6908 7725 .AND. (d2==iup_u .OR. d2==iup_l) ) RETURN 6909 IF ( (d==isouth_u .OR. d==isouth_l .OR. d==isouth_a) &7726 IF ( (d==isouth_u .OR. d==isouth_l ) & 6910 7727 .AND. (d2==isouth_u .OR. d2==isouth_l) ) RETURN 6911 IF ( (d==inorth_u .OR. d==inorth_l .OR. d==inorth_a) &7728 IF ( (d==inorth_u .OR. d==inorth_l ) & 6912 7729 .AND. (d2==inorth_u .OR. d2==inorth_l) ) RETURN 6913 IF ( (d==iwest_u .OR. d==iwest_l .OR. d==iwest_a) &7730 IF ( (d==iwest_u .OR. d==iwest_l ) & 6914 7731 .AND. (d2==iwest_u .OR. d2==iwest_l ) ) RETURN 6915 IF ( (d==ieast_u .OR. d==ieast_l .OR. d==ieast_a) &7732 IF ( (d==ieast_u .OR. d==ieast_l ) & 6916 7733 .AND. (d2==ieast_u .OR. d2==ieast_l ) ) RETURN 6917 7734 6918 7735 !-- second check: are surfaces facing away from each other 6919 7736 SELECT CASE (d) 6920 CASE (iup_u, iup_l , iup_a)!< upward facing surfaces7737 CASE (iup_u, iup_l) !< upward facing surfaces 6921 7738 IF ( z2 < z ) RETURN 6922 CASE (idown_a) !< downward facing surfaces 6923 IF ( z2 > z ) RETURN 6924 CASE (isouth_u, isouth_l, isouth_a) !< southward facing surfaces 7739 CASE (isouth_u, isouth_l) !< southward facing surfaces 6925 7740 IF ( y2 > y ) RETURN 6926 CASE (inorth_u, inorth_l , inorth_a)!< northward facing surfaces7741 CASE (inorth_u, inorth_l) !< northward facing surfaces 6927 7742 IF ( y2 < y ) RETURN 6928 CASE (iwest_u, iwest_l , iwest_a)!< westward facing surfaces7743 CASE (iwest_u, iwest_l) !< westward facing surfaces 6929 7744 IF ( x2 > x ) RETURN 6930 CASE (ieast_u, ieast_l , ieast_a)!< eastward facing surfaces7745 CASE (ieast_u, ieast_l) !< eastward facing surfaces 6931 7746 IF ( x2 < x ) RETURN 6932 7747 END SELECT … … 7143 7958 !> array for csf 7144 7959 !------------------------------------------------------------------------------! 7960 !-- quicksort.f -*-f90-*- 7961 !-- Author: t-nissie, adaptation J.Resler 7962 !-- License: GPLv3 7963 !-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea 7964 RECURSIVE SUBROUTINE quicksort_itarget(itarget, vffrac, ztransp, first, last) 7965 IMPLICIT NONE 7966 INTEGER(iwp), DIMENSION(:), INTENT(INOUT) :: itarget 7967 REAL(wp), DIMENSION(:), INTENT(INOUT) :: vffrac, ztransp 7968 INTEGER(iwp), INTENT(IN) :: first, last 7969 INTEGER(iwp) :: x, t 7970 INTEGER(iwp) :: i, j 7971 REAL(wp) :: tr 7972 7973 IF ( first>=last ) RETURN 7974 x = itarget((first+last)/2) 7975 i = first 7976 j = last 7977 DO 7978 DO WHILE ( itarget(i) < x ) 7979 i=i+1 7980 ENDDO 7981 DO WHILE ( x < itarget(j) ) 7982 j=j-1 7983 ENDDO 7984 IF ( i >= j ) EXIT 7985 t = itarget(i); itarget(i) = itarget(j); itarget(j) = t 7986 tr = vffrac(i); vffrac(i) = vffrac(j); vffrac(j) = tr 7987 tr = ztransp(i); ztransp(i) = ztransp(j); ztransp(j) = tr 7988 i=i+1 7989 j=j-1 7990 ENDDO 7991 IF ( first < i-1 ) CALL quicksort_itarget(itarget, vffrac, ztransp, first, i-1) 7992 IF ( j+1 < last ) CALL quicksort_itarget(itarget, vffrac, ztransp, j+1, last) 7993 END SUBROUTINE quicksort_itarget 7994 7145 7995 PURE FUNCTION svf_lt(svf1,svf2) result (res) 7146 7996 TYPE (t_svf), INTENT(in) :: svf1,svf2 … … 7153 8003 ENDIF 7154 8004 END FUNCTION svf_lt 7155 7156 8005 8006 7157 8007 !-- quicksort.f -*-f90-*- 7158 8008 !-- Author: t-nissie, adaptation J.Resler … … 7186 8036 END SUBROUTINE quicksort_svf 7187 8037 7188 7189 8038 PURE FUNCTION csf_lt(csf1,csf2) result (res) 7190 8039 TYPE (t_csf), INTENT(in) :: csf1,csf2 … … 7241 8090 INTEGER(iwp) :: iread, iwrite 7242 8091 TYPE(t_csf), DIMENSION(:), POINTER :: acsfnew 7243 !CHARACTER(100) :: msg8092 CHARACTER(100) :: msg 7244 8093 7245 8094 IF ( newsize == -1 ) THEN … … 7270 8119 .AND. acsfnew(iwrite)%itz == acsf(iread)%itz & 7271 8120 .AND. acsfnew(iwrite)%isurfs == acsf(iread)%isurfs ) THEN 7272 !-- We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray 7273 !-- probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that 7274 !-- might mean that the traced beam passes longer through the canopy box. 7275 IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf ) THEN 7276 acsfnew(iwrite)%rtransp = acsf(iread)%rtransp 7277 ENDIF 8121 7278 8122 acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf 7279 8123 !-- advance reading index, keep writing index … … 7310 8154 ncsfla = newsize 7311 8155 7312 !WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla7313 !CALL radiation_write_debug_log( msg )8156 WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla 8157 CALL radiation_write_debug_log( msg ) 7314 8158 7315 8159 END SUBROUTINE merge_and_grow_csf … … 7443 8287 SELECT CASE ( d ) 7444 8288 7445 CASE (iup_u,iup_l ,iup_a)!- gridbox up_facing face8289 CASE (iup_u,iup_l) !- gridbox up_facing face 7446 8290 sw_gridbox(1) = surfinsw(l) 7447 8291 lw_gridbox(1) = surfinlw(l) 7448 8292 swd_gridbox(1) = surfinswdif(l) 7449 8293 7450 CASE (idown_a) !- gridbox down_facing face 7451 sw_gridbox(2) = surfinsw(l) 7452 lw_gridbox(2) = surfinlw(l) 7453 swd_gridbox(2) = surfinswdif(l) 7454 7455 CASE (inorth_u,inorth_l,inorth_a) !- gridbox north_facing face 8294 CASE (inorth_u,inorth_l) !- gridbox north_facing face 7456 8295 sw_gridbox(3) = surfinsw(l) 7457 8296 lw_gridbox(3) = surfinlw(l) 7458 8297 swd_gridbox(3) = surfinswdif(l) 7459 8298 7460 CASE (isouth_u,isouth_l ,isouth_a)!- gridbox south_facing face8299 CASE (isouth_u,isouth_l) !- gridbox south_facing face 7461 8300 sw_gridbox(4) = surfinsw(l) 7462 8301 lw_gridbox(4) = surfinlw(l) 7463 8302 swd_gridbox(4) = surfinswdif(l) 7464 8303 7465 CASE (ieast_u,ieast_l ,ieast_a)!- gridbox east_facing face8304 CASE (ieast_u,ieast_l) !- gridbox east_facing face 7466 8305 sw_gridbox(5) = surfinsw(l) 7467 8306 lw_gridbox(5) = surfinlw(l) 7468 8307 swd_gridbox(5) = surfinswdif(l) 7469 8308 7470 CASE (iwest_u,iwest_l ,iwest_a)!- gridbox west_facing face8309 CASE (iwest_u,iwest_l) !- gridbox west_facing face 7471 8310 sw_gridbox(6) = surfinsw(l) 7472 8311 lw_gridbox(6) = surfinlw(l) … … 7520 8359 INTEGER(iwp) :: j !< 7521 8360 INTEGER(iwp) :: k !< 7522 INTEGER(iwp) :: m !< index of current surface element8361 INTEGER(iwp) :: l, m !< index of current surface element 7523 8362 7524 8363 IF ( mode == 'allocate' ) THEN … … 7604 8443 rad_sw_hr_av = 0.0_wp 7605 8444 8445 CASE ( 'rad_mrt_sw' ) 8446 IF ( .NOT. ALLOCATED( mrtinsw_av ) ) THEN 8447 ALLOCATE( mrtinsw_av(nmrtbl) ) 8448 ENDIF 8449 mrtinsw_av = 0.0_wp 8450 8451 CASE ( 'rad_mrt_lw' ) 8452 IF ( .NOT. ALLOCATED( mrtinlw_av ) ) THEN 8453 ALLOCATE( mrtinlw_av(nmrtbl) ) 8454 ENDIF 8455 mrtinlw_av = 0.0_wp 8456 8457 CASE ( 'rad_mrt' ) 8458 IF ( .NOT. ALLOCATED( mrt_av ) ) THEN 8459 ALLOCATE( mrt_av(nmrtbl) ) 8460 ENDIF 8461 mrt_av = 0.0_wp 8462 7606 8463 CASE DEFAULT 7607 8464 CONTINUE … … 7819 8676 ENDIF 7820 8677 8678 CASE ( 'rad_mrt_sw' ) 8679 IF ( ALLOCATED( mrtinsw_av ) ) THEN 8680 mrtinsw_av(:) = mrtinsw_av(:) + mrtinsw(:) 8681 ENDIF 8682 8683 CASE ( 'rad_mrt_lw' ) 8684 IF ( ALLOCATED( mrtinlw_av ) ) THEN 8685 mrtinlw_av(:) = mrtinlw_av(:) + mrtinlw(:) 8686 ENDIF 8687 8688 CASE ( 'rad_mrt' ) 8689 IF ( ALLOCATED( mrt_av ) ) THEN 8690 mrt_av(:) = mrt_av(:) + mrt(:) 8691 ENDIF 8692 7821 8693 CASE DEFAULT 7822 8694 CONTINUE … … 7828 8700 SELECT CASE ( TRIM( variable ) ) 7829 8701 7830 CASE ( 'rad_net*' )8702 CASE ( 'rad_net*' ) 7831 8703 IF ( ALLOCATED( rad_net_av ) ) THEN 7832 8704 DO i = nxlg, nxrg … … 7974 8846 ENDIF 7975 8847 8848 CASE ( 'rad_mrt_sw' ) 8849 IF ( ALLOCATED( mrtinsw_av ) ) THEN 8850 mrtinsw_av(:) = mrtinsw_av(:) / REAL( average_count_3d, KIND=wp ) 8851 ENDIF 8852 8853 CASE ( 'rad_mrt_lw' ) 8854 IF ( ALLOCATED( mrtinlw_av ) ) THEN 8855 mrtinlw_av(:) = mrtinlw_av(:) / REAL( average_count_3d, KIND=wp ) 8856 ENDIF 8857 8858 CASE ( 'rad_mrt' ) 8859 IF ( ALLOCATED( mrt_av ) ) THEN 8860 mrt_av(:) = mrt_av(:) / REAL( average_count_3d, KIND=wp ) 8861 ENDIF 8862 7976 8863 END SELECT 7977 8864 … … 8009 8896 'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz', & 8010 8897 'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz', & 8011 'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' ) 8898 'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz', & 8899 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' ) 8012 8900 grid_x = 'x' 8013 8901 grid_y = 'y' … … 8037 8925 ! Description: 8038 8926 ! ------------ 8039 !> Subroutine defining 3D output variables8927 !> Subroutine defining 2D output variables 8040 8928 !------------------------------------------------------------------------------! 8041 8929 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, & … … 8456 9344 CHARACTER (LEN=*) :: variable !< 8457 9345 8458 INTEGER(iwp) :: av !< 8459 INTEGER(iwp) :: i !< 8460 INTEGER(iwp) :: j !< 8461 INTEGER(iwp) :: k !< 8462 INTEGER(iwp) :: nzb_do !< 8463 INTEGER(iwp) :: nzt_do !< 8464 8465 LOGICAL :: found !< 8466 8467 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 8468 8469 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 8470 9346 INTEGER(iwp) :: av !< 9347 INTEGER(iwp) :: i, j, k, l !< 9348 INTEGER(iwp) :: nzb_do !< 9349 INTEGER(iwp) :: nzt_do !< 9350 9351 LOGICAL :: found !< 9352 9353 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 9354 9355 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 8471 9356 8472 9357 found = .TRUE. … … 8657 9542 ENDDO 8658 9543 ENDDO 9544 ENDIF 9545 9546 CASE ( 'rad_mrt_sw' ) 9547 local_pf = REAL( fill_value, KIND = wp ) 9548 IF ( av == 0 ) THEN 9549 DO l = 1, nmrtbl 9550 i = mrtbl(ix,l) 9551 j = mrtbl(iy,l) 9552 k = mrtbl(iz,l) 9553 local_pf(i,j,k) = mrtinsw(l) 9554 ENDDO 9555 ELSE 9556 IF ( ALLOCATED( mrtinsw_av ) ) THEN 9557 DO l = 1, nmrtbl 9558 i = mrtbl(ix,l) 9559 j = mrtbl(iy,l) 9560 k = mrtbl(iz,l) 9561 local_pf(i,j,k) = mrtinsw_av(l) 9562 ENDDO 9563 ENDIF 9564 ENDIF 9565 9566 CASE ( 'rad_mrt_lw' ) 9567 local_pf = REAL( fill_value, KIND = wp ) 9568 IF ( av == 0 ) THEN 9569 DO l = 1, nmrtbl 9570 i = mrtbl(ix,l) 9571 j = mrtbl(iy,l) 9572 k = mrtbl(iz,l) 9573 local_pf(i,j,k) = mrtinlw(l) 9574 ENDDO 9575 ELSE 9576 IF ( ALLOCATED( mrtinlw_av ) ) THEN 9577 DO l = 1, nmrtbl 9578 i = mrtbl(ix,l) 9579 j = mrtbl(iy,l) 9580 k = mrtbl(iz,l) 9581 local_pf(i,j,k) = mrtinlw_av(l) 9582 ENDDO 9583 ENDIF 9584 ENDIF 9585 9586 CASE ( 'rad_mrt' ) 9587 local_pf = REAL( fill_value, KIND = wp ) 9588 IF ( av == 0 ) THEN 9589 DO l = 1, nmrtbl 9590 i = mrtbl(ix,l) 9591 j = mrtbl(iy,l) 9592 k = mrtbl(iz,l) 9593 local_pf(i,j,k) = mrt(l) 9594 ENDDO 9595 ELSE 9596 IF ( ALLOCATED( mrt_av ) ) THEN 9597 DO l = 1, nmrtbl 9598 i = mrtbl(ix,l) 9599 j = mrtbl(iy,l) 9600 k = mrtbl(iz,l) 9601 local_pf(i,j,k) = mrt_av(l) 9602 ENDDO 9603 ENDIF 8659 9604 ENDIF 8660 9605 … … 9342 10287 !> Subroutine writes debug information 9343 10288 !------------------------------------------------------------------------------! 9344 !SUBROUTINE radiation_write_debug_log ( message )10289 SUBROUTINE radiation_write_debug_log ( message ) 9345 10290 !> it writes debug log with time stamp 9346 !CHARACTER(*) :: message9347 !CHARACTER(15) :: dtc9348 !CHARACTER(8) :: date9349 !CHARACTER(10) :: time9350 !CHARACTER(5) :: zone9351 !CALL date_and_time(date, time, zone)9352 !dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)9353 !WRITE(9,'(2A)') dtc, TRIM(message)9354 !FLUSH(9)9355 !END SUBROUTINE radiation_write_debug_log10291 CHARACTER(*) :: message 10292 CHARACTER(15) :: dtc 10293 CHARACTER(8) :: date 10294 CHARACTER(10) :: time 10295 CHARACTER(5) :: zone 10296 CALL date_and_time(date, time, zone) 10297 dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10) 10298 WRITE(9,'(2A)') dtc, TRIM(message) 10299 FLUSH(9) 10300 END SUBROUTINE radiation_write_debug_log 9356 10301 9357 10302 END MODULE radiation_model_mod -
palm/trunk/SOURCE/sum_up_3d_data.f90
r3294 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Add biometeorology, 29 ! fix chemistry output call, 30 ! move usm calls 31 ! 32 ! 3294 2018-10-01 02:37:10Z raasch 27 33 ! changes concerning modularization of ocean option 28 34 ! … … 260 266 USE radiation_model_mod, & 261 267 ONLY: radiation, radiation_3d_data_averaging 268 269 USE biometeorology_mod, & 270 ONLY: biometeorology_3d_data_averaging 262 271 263 272 USE surface_mod, & … … 305 314 306 315 DO ii = 1, doav_n 307 ! 308 !-- Temporary solution to account for data output within the new urban 309 !-- surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar ) 316 310 317 trimvar = TRIM( doav(ii) ) 311 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN 312 trimvar = 'usm_output' 313 ENDIF 314 318 315 319 SELECT CASE ( trimvar ) 316 320 … … 495 499 z0q_av = 0.0_wp 496 500 497 CASE ( 'usm_output' )498 !499 !-- Block of urban surface model outputs500 CALL usm_average_3d_data( 'allocate', doav(ii) )501 502 501 503 502 CASE DEFAULT … … 505 504 ! 506 505 !-- Allocating and initializing data arrays for other modules 506 507 IF ( air_chemistry .AND. & 508 (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') ) THEN 509 CALL chem_3d_data_averaging( 'allocate', doav(ii) ) 510 ENDIF 511 507 512 IF ( bulk_cloud_model ) THEN 508 513 CALL bcm_3d_data_averaging( 'allocate', doav(ii) ) 509 514 ENDIF 510 515 511 IF ( air_chemistry .AND. trimvar(1:3) == 'kc_') THEN512 CALL chem_3d_data_averaging( 'allocate', doav(ii) )516 IF ( radiation .AND. trimvar(1:4) == 'bio_' ) THEN 517 CALL biometeorology_3d_data_averaging( 'allocate', doav(ii) ) 513 518 ENDIF 514 519 … … 529 534 ENDIF 530 535 536 CALL tcm_3d_data_averaging( 'allocate', doav(ii) ) 537 538 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN 539 CALL usm_average_3d_data( 'allocate', doav(ii) ) 540 ENDIF 541 531 542 IF ( uv_exposure .AND. trimvar(1:5) == 'uvem_') THEN 532 543 CALL uvem_3d_data_averaging( 'allocate', doav(ii) ) 533 544 ENDIF 534 545 535 CALL tcm_3d_data_averaging( 'allocate', doav(ii) )536 537 546 ! 538 547 !-- User-defined quantities … … 548 557 !-- Loop of all variables to be averaged. 549 558 DO ii = 1, doav_n 550 ! 551 !-- Temporary solution to account for data output within the new urban 552 !-- surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar ) 553 trimvar = TRIM( doav(ii) ) 554 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN 555 trimvar = 'usm_output' 556 ENDIF 559 560 trimvar = TRIM( doav(ii) ) 557 561 ! 558 562 !-- Store the array chosen on the temporary array. … … 1146 1150 ENDIF 1147 1151 1148 CASE ( 'usm_output' ) 1149 ! 1150 !-- Block of urban surface model outputs. 1152 CASE DEFAULT 1153 ! 1154 !-- Summing up data from other modules 1155 IF ( bulk_cloud_model ) THEN 1156 CALL bcm_3d_data_averaging( 'sum', doav(ii) ) 1157 ENDIF 1158 1159 IF ( air_chemistry .AND. & 1160 (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') ) THEN 1161 CALL chem_3d_data_averaging( 'sum',doav(ii) ) 1162 ENDIF 1163 1164 IF ( radiation .AND. trimvar(1:4) == 'bio_') THEN 1165 CALL biometeorology_3d_data_averaging( 'sum', doav(ii) ) 1166 ENDIF 1167 1168 IF ( gust_module_enabled ) THEN 1169 CALL gust_3d_data_averaging( 'sum', doav(ii) ) 1170 ENDIF 1171 1172 IF ( land_surface ) THEN 1173 CALL lsm_3d_data_averaging( 'sum', doav(ii) ) 1174 ENDIF 1175 1176 IF ( ocean_mode ) THEN 1177 CALL ocean_3d_data_averaging( 'sum', doav(ii) ) 1178 ENDIF 1179 1180 IF ( radiation ) THEN 1181 CALL radiation_3d_data_averaging( 'sum', doav(ii) ) 1182 ENDIF 1183 1184 CALL tcm_3d_data_averaging( 'sum', doav(ii) ) 1185 1151 1186 !-- In case of urban surface variables it should be always checked 1152 1187 !-- if respective arrays are allocated, at least in case of a restart 1153 1188 !-- run, as averaged usm arrays are not read from file at the moment. 1154 CALL usm_average_3d_data( 'allocate', doav(ii) ) 1155 CALL usm_average_3d_data( 'sum', doav(ii) ) 1156 1157 CASE DEFAULT 1158 ! 1159 !-- Summing up data from other modules 1160 IF ( bulk_cloud_model ) THEN 1161 CALL bcm_3d_data_averaging( 'sum', doav(ii) ) 1162 ENDIF 1163 1164 IF ( air_chemistry .AND. trimvar(1:3) == 'kc_') THEN 1165 CALL chem_3d_data_averaging( 'sum',doav(ii) ) 1166 ENDIF 1167 1168 IF ( gust_module_enabled ) THEN 1169 CALL gust_3d_data_averaging( 'sum', doav(ii) ) 1170 ENDIF 1171 1172 IF ( land_surface ) THEN 1173 CALL lsm_3d_data_averaging( 'sum', doav(ii) ) 1174 ENDIF 1175 1176 IF ( ocean_mode ) THEN 1177 CALL ocean_3d_data_averaging( 'sum', doav(ii) ) 1178 ENDIF 1179 1180 IF ( radiation ) THEN 1181 CALL radiation_3d_data_averaging( 'sum', doav(ii) ) 1182 ENDIF 1183 1184 CALL tcm_3d_data_averaging( 'sum', doav(ii) ) 1189 IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN 1190 CALL usm_average_3d_data( 'allocate', doav(ii) ) 1191 CALL usm_average_3d_data( 'sum', doav(ii) ) 1192 ENDIF 1185 1193 1186 1194 IF ( uv_exposure ) THEN -
palm/trunk/SOURCE/time_integration_spinup.f90
-
Property
svn:mergeinfo
set to
(toggle deleted branches)
/palm/branches/chemistry/SOURCE/time_integration_spinup.f90 2047-3190,3218-3297 /palm/branches/palm4u/SOURCE/time_integration_spinup.f90 2540-2692 /palm/branches/rans/SOURCE/time_integration_spinup.f90 2078-3128 /palm/branches/resler/SOURCE/time_integration_spinup.f90 2023-3336 /palm/branches/forwind/SOURCE/time_integration_spinup.f90 1564-1913 /palm/branches/fricke/SOURCE/time_integration_spinup.f90 942-977 /palm/branches/hoffmann/SOURCE/time_integration_spinup.f90 989-1052 /palm/branches/letzel/masked_output/SOURCE/time_integration_spinup.f90 296-409 /palm/branches/suehring/time_integration_spinup.f90 423-666
r3274 r3337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! (from branch resler) 28 ! Add pt1 initialization 29 ! 30 ! 3274 2018-09-24 15:42:55Z knoop 27 31 ! Modularization of all bulk cloud physics code components 28 32 ! … … 295 299 k = surf_usm_h%k(m) 296 300 pt(k,j,i) = pt_spinup 301 !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!! 302 surf_usm_h%pt1 = pt_spinup 303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 297 304 ENDDO 298 305 … … 303 310 k = surf_usm_v(l)%k(m) 304 311 pt(k,j,i) = pt_spinup 312 !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!! 313 surf_usm_v(l)%pt1 = pt_spinup 314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 305 315 ENDDO 306 316 ENDDO -
Property
svn:mergeinfo
set to
(toggle deleted branches)
-
palm/trunk/SOURCE/urban_surface_mod.f90
r3274 r3337 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw* 31 ! 32 ! 3274 2018-09-24 15:42:55Z knoop 30 33 ! Modularization of all bulk cloud physics code components 31 34 ! … … 349 352 spinup_pt_mean, spinup_time, time_do3d, dt_do3d, & 350 353 average_count_3d, varnamelength, urban_surface, & 351 plant_canopy 354 plant_canopy, dz 352 355 353 356 USE cpulog, & … … 381 384 surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff, & 382 385 rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, & 383 surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,&384 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,&386 surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir, & 387 pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, & 385 388 inorth_l, isouth_l, ieast_l, iwest_l, id, & 386 iz, iy, ix, idir, jdir, kdir, nsurf_type, nsurf, idsvf, ndsvf, & 387 iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a, & 389 iz, iy, ix, nsurf, idsvf, ndsvf, & 388 390 idcsf, ndcsf, kdcsf, pct, & 389 startland, endland, startwall, endwall, skyvf, skyvft 391 startland, endland, startwall, endwall, skyvf, skyvft, nzub, & 392 nzut, nzpt, npcbl, pcbl 390 393 391 394 USE statistics, & … … 414 417 INTEGER(iwp) :: pedestrian_category = 2 !< default category for wall surface in pedestrian zone 415 418 INTEGER(iwp) :: roof_category = 2 !< default category for root surface 419 REAL(wp) :: roughness_concrete = 0.001_wp !< roughness length of average concrete surface 416 420 ! 417 421 !-- Indices of input attributes for (above) ground floor level … … 605 609 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection 606 610 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection 611 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw_av !< Average of pcbinlw 612 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw_av !< Average of pcbinsw 613 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdir_av !< Average of pcbinswdir 614 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdif_av !< Average of pcbinswdif 615 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswref_av !< Average of pcbinswref 607 616 608 617 … … 1260 1269 !-- find the real name of the variable 1261 1270 ids = -1 1271 l = -1 1262 1272 var = TRIM(variable) 1263 1273 DO i = 0, nd-1 1264 1274 k = len(TRIM(var)) 1265 1275 j = len(TRIM(dirname(i))) 1266 IF ( var(k-j+1:k) == dirname(i) ) THEN1276 IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN 1267 1277 ids = i 1268 1278 idsint = dirint(ids) … … 1271 1281 ENDIF 1272 1282 ENDDO 1283 l = idsint - 2 ! horisontal direction index - terible hack ! 1284 IF ( l < 0 .OR. l > 3 ) THEN 1285 l = -1 1286 END IF 1273 1287 IF ( ids == -1 ) THEN 1274 1288 var = TRIM(variable) … … 1311 1325 CASE ( 'usm_rad_net' ) 1312 1326 !-- array of complete radiation balance 1313 IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN1327 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN 1314 1328 ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) ) 1315 1329 surf_usm_h%rad_net_av = 0.0_wp 1330 ELSE 1331 IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN 1332 ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) ) 1333 surf_usm_v(l)%rad_net_av = 0.0_wp 1334 ENDIF 1316 1335 ENDIF 1317 DO l = 0, 31318 IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN1319 ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )1320 surf_usm_v(l)%rad_net_av = 0.0_wp1321 ENDIF1322 ENDDO1323 1336 1324 1337 CASE ( 'usm_rad_insw' ) … … 1398 1411 ENDIF 1399 1412 1413 CASE ( 'usm_rad_pc_inlw' ) 1414 !-- array of of lw radiation absorbed in plant canopy 1415 IF ( .NOT. ALLOCATED(pcbinlw_av) ) THEN 1416 ALLOCATE( pcbinlw_av(1:npcbl) ) 1417 pcbinlw_av = 0.0_wp 1418 ENDIF 1419 1420 CASE ( 'usm_rad_pc_insw' ) 1421 !-- array of of sw radiation absorbed in plant canopy 1422 IF ( .NOT. ALLOCATED(pcbinsw_av) ) THEN 1423 ALLOCATE( pcbinsw_av(1:npcbl) ) 1424 pcbinsw_av = 0.0_wp 1425 ENDIF 1426 1427 CASE ( 'usm_rad_pc_inswdir' ) 1428 !-- array of of direct sw radiation absorbed in plant canopy 1429 IF ( .NOT. ALLOCATED(pcbinswdir_av) ) THEN 1430 ALLOCATE( pcbinswdir_av(1:npcbl) ) 1431 pcbinswdir_av = 0.0_wp 1432 ENDIF 1433 1434 CASE ( 'usm_rad_pc_inswdif' ) 1435 !-- array of of diffuse sw radiation absorbed in plant canopy 1436 IF ( .NOT. ALLOCATED(pcbinswdif_av) ) THEN 1437 ALLOCATE( pcbinswdif_av(1:npcbl) ) 1438 pcbinswdif_av = 0.0_wp 1439 ENDIF 1440 1441 CASE ( 'usm_rad_pc_inswref' ) 1442 !-- array of of reflected sw radiation absorbed in plant canopy 1443 IF ( .NOT. ALLOCATED(pcbinswref_av) ) THEN 1444 ALLOCATE( pcbinswref_av(1:npcbl) ) 1445 pcbinswref_av = 0.0_wp 1446 ENDIF 1447 1400 1448 CASE ( 'usm_rad_hf' ) 1401 1449 !-- array of heat flux from radiation for surfaces after i-th reflection 1402 IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN1450 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN 1403 1451 ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) ) 1404 1452 surf_usm_h%surfhf_av = 0.0_wp 1405 ENDIF 1406 DO l = 0, 3 1453 ELSE 1407 1454 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN 1408 1455 ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) ) 1409 1456 surf_usm_v(l)%surfhf_av = 0.0_wp 1410 1457 ENDIF 1411 END DO1458 ENDIF 1412 1459 1413 1460 CASE ( 'usm_wshf' ) 1414 1461 !-- array of sensible heat flux from surfaces 1415 1462 !-- land surfaces 1416 IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN1463 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN 1417 1464 ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) ) 1418 1465 surf_usm_h%wshf_eb_av = 0.0_wp 1419 ENDIF 1420 DO l = 0, 3 1466 ELSE 1421 1467 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN 1422 1468 ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) ) 1423 1469 surf_usm_v(l)%wshf_eb_av = 0.0_wp 1424 1470 ENDIF 1425 END DO1471 ENDIF 1426 1472 ! 1427 1473 !-- Please note, the following output quantities belongs to the … … 1431 1477 CASE ( 'usm_wghf' ) 1432 1478 !-- array of heat flux from ground (wall, roof, land) 1433 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN1479 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN 1434 1480 ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) ) 1435 1481 surf_usm_h%wghf_eb_av = 0.0_wp 1436 ENDIF 1437 DO l = 0, 3 1482 ELSE 1438 1483 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN 1439 1484 ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) ) 1440 1485 surf_usm_v(l)%wghf_eb_av = 0.0_wp 1441 1486 ENDIF 1442 END DO1487 ENDIF 1443 1488 1444 1489 CASE ( 'usm_wghf_window' ) 1445 1490 !-- array of heat flux from window ground (wall, roof, land) 1446 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN1491 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN 1447 1492 ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) ) 1448 1493 surf_usm_h%wghf_eb_window_av = 0.0_wp 1449 ENDIF 1450 DO l = 0, 3 1494 ELSE 1451 1495 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN 1452 1496 ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) ) 1453 1497 surf_usm_v(l)%wghf_eb_window_av = 0.0_wp 1454 1498 ENDIF 1455 END DO1499 ENDIF 1456 1500 1457 1501 CASE ( 'usm_wghf_green' ) 1458 1502 !-- array of heat flux from green ground (wall, roof, land) 1459 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN1503 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN 1460 1504 ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) ) 1461 1505 surf_usm_h%wghf_eb_green_av = 0.0_wp 1462 ENDIF 1463 DO l = 0, 3 1506 ELSE 1464 1507 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN 1465 1508 ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) ) 1466 1509 surf_usm_v(l)%wghf_eb_green_av = 0.0_wp 1467 1510 ENDIF 1468 END DO1511 ENDIF 1469 1512 1470 1513 CASE ( 'usm_iwghf' ) 1471 1514 !-- array of heat flux from indoor ground (wall, roof, land) 1472 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN1515 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN 1473 1516 ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) ) 1474 1517 surf_usm_h%iwghf_eb_av = 0.0_wp 1475 ENDIF 1476 DO l = 0, 3 1518 ELSE 1477 1519 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN 1478 1520 ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) ) 1479 1521 surf_usm_v(l)%iwghf_eb_av = 0.0_wp 1480 1522 ENDIF 1481 END DO1523 ENDIF 1482 1524 1483 1525 CASE ( 'usm_iwghf_window' ) 1484 1526 !-- array of heat flux from indoor window ground (wall, roof, land) 1485 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN1527 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN 1486 1528 ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) ) 1487 1529 surf_usm_h%iwghf_eb_window_av = 0.0_wp 1488 ENDIF 1489 DO l = 0, 3 1530 ELSE 1490 1531 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN 1491 1532 ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) ) 1492 1533 surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp 1493 1534 ENDIF 1494 END DO1535 ENDIF 1495 1536 1496 1537 CASE ( 'usm_t_surf' ) 1497 1538 !-- surface temperature for surfaces 1498 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN1539 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN 1499 1540 ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) ) 1500 1541 surf_usm_h%t_surf_av = 0.0_wp 1501 ENDIF 1502 DO l = 0, 3 1542 ELSE 1503 1543 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_av) ) THEN 1504 1544 ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) ) 1505 1545 surf_usm_v(l)%t_surf_av = 0.0_wp 1506 1546 ENDIF 1507 END DO1547 ENDIF 1508 1548 1509 1549 CASE ( 'usm_t_surf_window' ) 1510 1550 !-- surface temperature for window surfaces 1511 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN1551 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN 1512 1552 ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) ) 1513 1553 surf_usm_h%t_surf_window_av = 0.0_wp 1514 ENDIF 1515 DO l = 0, 3 1554 ELSE 1516 1555 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN 1517 1556 ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) ) 1518 1557 surf_usm_v(l)%t_surf_window_av = 0.0_wp 1519 1558 ENDIF 1520 END DO1559 ENDIF 1521 1560 1522 1561 CASE ( 'usm_t_surf_green' ) 1523 1562 !-- surface temperature for green surfaces 1524 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN1563 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN 1525 1564 ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) ) 1526 1565 surf_usm_h%t_surf_green_av = 0.0_wp 1527 ENDIF 1528 DO l = 0, 3 1566 ELSE 1529 1567 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN 1530 1568 ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) ) 1531 1569 surf_usm_v(l)%t_surf_green_av = 0.0_wp 1532 1570 ENDIF 1533 END DO1571 ENDIF 1534 1572 1535 1573 CASE ( 'usm_t_surf_10cm' ) 1536 1574 !-- near surface temperature for whole surfaces 1537 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN1575 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN 1538 1576 ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) ) 1539 1577 surf_usm_h%t_surf_10cm_av = 0.0_wp 1540 ENDIF 1541 DO l = 0, 3 1578 ELSE 1542 1579 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) ) THEN 1543 1580 ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) ) 1544 1581 surf_usm_v(l)%t_surf_10cm_av = 0.0_wp 1545 1582 ENDIF 1546 END DO1583 ENDIF 1547 1584 1548 1585 CASE ( 'usm_t_wall' ) 1549 1586 !-- wall temperature for iwl layer of walls and land 1550 IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN1587 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN 1551 1588 ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1552 1589 surf_usm_h%t_wall_av = 0.0_wp 1553 ENDIF 1554 DO l = 0, 3 1590 ELSE 1555 1591 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN 1556 1592 ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1557 1593 surf_usm_v(l)%t_wall_av = 0.0_wp 1558 1594 ENDIF 1559 END DO1595 ENDIF 1560 1596 1561 1597 CASE ( 'usm_t_window' ) 1562 1598 !-- window temperature for iwl layer of walls and land 1563 IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN1599 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN 1564 1600 ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1565 1601 surf_usm_h%t_window_av = 0.0_wp 1566 ENDIF 1567 DO l = 0, 3 1602 ELSE 1568 1603 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN 1569 1604 ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1570 1605 surf_usm_v(l)%t_window_av = 0.0_wp 1571 1606 ENDIF 1572 END DO1607 ENDIF 1573 1608 1574 1609 CASE ( 'usm_t_green' ) 1575 1610 !-- green temperature for iwl layer of walls and land 1576 IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN1611 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN 1577 1612 ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1578 1613 surf_usm_h%t_green_av = 0.0_wp 1579 ENDIF 1580 DO l = 0, 3 1614 ELSE 1581 1615 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN 1582 1616 ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1583 1617 surf_usm_v(l)%t_green_av = 0.0_wp 1584 1618 ENDIF 1585 END DO1619 ENDIF 1586 1620 1587 1621 CASE DEFAULT … … 1596 1630 CASE ( 'usm_rad_net' ) 1597 1631 !-- array of complete radiation balance 1598 DO m = 1, surf_usm_h%ns 1599 surf_usm_h%rad_net_av(m) = & 1600 surf_usm_h%rad_net_av(m) + & 1601 surf_usm_h%rad_net_l(m) 1602 ENDDO 1603 DO l = 0, 3 1632 IF ( l == -1 ) THEN 1633 DO m = 1, surf_usm_h%ns 1634 surf_usm_h%rad_net_av(m) = & 1635 surf_usm_h%rad_net_av(m) + & 1636 surf_usm_h%rad_net_l(m) 1637 ENDDO 1638 ELSE 1604 1639 DO m = 1, surf_usm_v(l)%ns 1605 1640 surf_usm_v(l)%rad_net_av(m) = & … … 1607 1642 surf_usm_v(l)%rad_net_l(m) 1608 1643 ENDDO 1609 END DO1610 1644 ENDIF 1645 1611 1646 CASE ( 'usm_rad_insw' ) 1612 1647 !-- array of sw radiation falling to surface after i-th reflection … … 1700 1735 ENDDO 1701 1736 1737 CASE ( 'usm_rad_pc_inlw' ) 1738 pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:) 1739 1740 CASE ( 'usm_rad_pc_insw' ) 1741 pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:) 1742 1743 CASE ( 'usm_rad_pc_inswdir' ) 1744 pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:) 1745 1746 CASE ( 'usm_rad_pc_inswdif' ) 1747 pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:) 1748 1749 CASE ( 'usm_rad_pc_inswref' ) 1750 pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:) & 1751 - pcbinswdir(:) & 1752 - pcbinswdif(:) 1753 1702 1754 CASE ( 'usm_rad_hf' ) 1703 1755 !-- array of heat flux from radiation for surfaces after i-th reflection 1704 DO m = 1, surf_usm_h%ns 1705 surf_usm_h%surfhf_av(m) = & 1706 surf_usm_h%surfhf_av(m) + & 1707 surf_usm_h%surfhf(m) 1708 ENDDO 1709 DO l = 0, 3 1756 IF ( l == -1 ) THEN 1757 DO m = 1, surf_usm_h%ns 1758 surf_usm_h%surfhf_av(m) = & 1759 surf_usm_h%surfhf_av(m) + & 1760 surf_usm_h%surfhf(m) 1761 ENDDO 1762 ELSE 1710 1763 DO m = 1, surf_usm_v(l)%ns 1711 1764 surf_usm_v(l)%surfhf_av(m) = & … … 1713 1766 surf_usm_v(l)%surfhf(m) 1714 1767 ENDDO 1715 END DO1768 ENDIF 1716 1769 1717 1770 CASE ( 'usm_wshf' ) 1718 1771 !-- array of sensible heat flux from surfaces (land, roof, wall) 1719 DO m = 1, surf_usm_h%ns 1720 surf_usm_h%wshf_eb_av(m) = & 1721 surf_usm_h%wshf_eb_av(m) + & 1722 surf_usm_h%wshf_eb(m) 1723 ENDDO 1724 DO l = 0, 3 1772 IF ( l == -1 ) THEN 1773 DO m = 1, surf_usm_h%ns 1774 surf_usm_h%wshf_eb_av(m) = & 1775 surf_usm_h%wshf_eb_av(m) + & 1776 surf_usm_h%wshf_eb(m) 1777 ENDDO 1778 ELSE 1725 1779 DO m = 1, surf_usm_v(l)%ns 1726 1780 surf_usm_v(l)%wshf_eb_av(m) = & … … 1728 1782 surf_usm_v(l)%wshf_eb(m) 1729 1783 ENDDO 1730 END DO1784 ENDIF 1731 1785 1732 1786 CASE ( 'usm_wghf' ) 1733 1787 !-- array of heat flux from ground (wall, roof, land) 1734 DO m = 1, surf_usm_h%ns 1735 surf_usm_h%wghf_eb_av(m) = & 1736 surf_usm_h%wghf_eb_av(m) + & 1737 surf_usm_h%wghf_eb(m) 1738 ENDDO 1739 DO l = 0, 3 1788 IF ( l == -1 ) THEN 1789 DO m = 1, surf_usm_h%ns 1790 surf_usm_h%wghf_eb_av(m) = & 1791 surf_usm_h%wghf_eb_av(m) + & 1792 surf_usm_h%wghf_eb(m) 1793 ENDDO 1794 ELSE 1740 1795 DO m = 1, surf_usm_v(l)%ns 1741 1796 surf_usm_v(l)%wghf_eb_av(m) = & … … 1743 1798 surf_usm_v(l)%wghf_eb(m) 1744 1799 ENDDO 1745 END DO1800 ENDIF 1746 1801 1747 1802 CASE ( 'usm_wghf_window' ) 1748 1803 !-- array of heat flux from window ground (wall, roof, land) 1749 DO m = 1, surf_usm_h%ns 1750 surf_usm_h%wghf_eb_window_av(m) = & 1751 surf_usm_h%wghf_eb_window_av(m) + & 1752 surf_usm_h%wghf_eb_window(m) 1753 ENDDO 1754 DO l = 0, 3 1804 IF ( l == -1 ) THEN 1805 DO m = 1, surf_usm_h%ns 1806 surf_usm_h%wghf_eb_window_av(m) = & 1807 surf_usm_h%wghf_eb_window_av(m) + & 1808 surf_usm_h%wghf_eb_window(m) 1809 ENDDO 1810 ELSE 1755 1811 DO m = 1, surf_usm_v(l)%ns 1756 1812 surf_usm_v(l)%wghf_eb_window_av(m) = & … … 1758 1814 surf_usm_v(l)%wghf_eb_window(m) 1759 1815 ENDDO 1760 END DO1816 ENDIF 1761 1817 1762 1818 CASE ( 'usm_wghf_green' ) 1763 1819 !-- array of heat flux from green ground (wall, roof, land) 1764 DO m = 1, surf_usm_h%ns 1765 surf_usm_h%wghf_eb_green_av(m) = & 1766 surf_usm_h%wghf_eb_green_av(m) + & 1767 surf_usm_h%wghf_eb_green(m) 1768 ENDDO 1769 DO l = 0, 3 1820 IF ( l == -1 ) THEN 1821 DO m = 1, surf_usm_h%ns 1822 surf_usm_h%wghf_eb_green_av(m) = & 1823 surf_usm_h%wghf_eb_green_av(m) + & 1824 surf_usm_h%wghf_eb_green(m) 1825 ENDDO 1826 ELSE 1770 1827 DO m = 1, surf_usm_v(l)%ns 1771 1828 surf_usm_v(l)%wghf_eb_green_av(m) = & … … 1773 1830 surf_usm_v(l)%wghf_eb_green(m) 1774 1831 ENDDO 1775 END DO1832 ENDIF 1776 1833 1777 1834 CASE ( 'usm_iwghf' ) 1778 1835 !-- array of heat flux from indoor ground (wall, roof, land) 1779 DO m = 1, surf_usm_h%ns 1780 surf_usm_h%iwghf_eb_av(m) = & 1781 surf_usm_h%iwghf_eb_av(m) + & 1782 surf_usm_h%iwghf_eb(m) 1783 ENDDO 1784 DO l = 0, 3 1836 IF ( l == -1 ) THEN 1837 DO m = 1, surf_usm_h%ns 1838 surf_usm_h%iwghf_eb_av(m) = & 1839 surf_usm_h%iwghf_eb_av(m) + & 1840 surf_usm_h%iwghf_eb(m) 1841 ENDDO 1842 ELSE 1785 1843 DO m = 1, surf_usm_v(l)%ns 1786 1844 surf_usm_v(l)%iwghf_eb_av(m) = & … … 1788 1846 surf_usm_v(l)%iwghf_eb(m) 1789 1847 ENDDO 1790 END DO1848 ENDIF 1791 1849 1792 1850 CASE ( 'usm_iwghf_window' ) 1793 1851 !-- array of heat flux from indoor window ground (wall, roof, land) 1794 DO m = 1, surf_usm_h%ns 1795 surf_usm_h%iwghf_eb_window_av(m) = & 1796 surf_usm_h%iwghf_eb_window_av(m) + & 1797 surf_usm_h%iwghf_eb_window(m) 1798 ENDDO 1799 DO l = 0, 3 1852 IF ( l == -1 ) THEN 1853 DO m = 1, surf_usm_h%ns 1854 surf_usm_h%iwghf_eb_window_av(m) = & 1855 surf_usm_h%iwghf_eb_window_av(m) + & 1856 surf_usm_h%iwghf_eb_window(m) 1857 ENDDO 1858 ELSE 1800 1859 DO m = 1, surf_usm_v(l)%ns 1801 1860 surf_usm_v(l)%iwghf_eb_window_av(m) = & … … 1803 1862 surf_usm_v(l)%iwghf_eb_window(m) 1804 1863 ENDDO 1805 END DO1864 ENDIF 1806 1865 1807 1866 CASE ( 'usm_t_surf' ) 1808 1867 !-- surface temperature for surfaces 1809 DO m = 1, surf_usm_h%ns 1810 surf_usm_h%t_surf_av(m) = & 1811 surf_usm_h%t_surf_av(m) + & 1812 t_surf_h(m) 1813 ENDDO 1814 DO l = 0, 3 1868 IF ( l == -1 ) THEN 1869 DO m = 1, surf_usm_h%ns 1870 surf_usm_h%t_surf_av(m) = & 1871 surf_usm_h%t_surf_av(m) + & 1872 t_surf_h(m) 1873 ENDDO 1874 ELSE 1815 1875 DO m = 1, surf_usm_v(l)%ns 1816 1876 surf_usm_v(l)%t_surf_av(m) = & … … 1818 1878 t_surf_v(l)%t(m) 1819 1879 ENDDO 1820 END DO1880 ENDIF 1821 1881 1822 1882 CASE ( 'usm_t_surf_window' ) 1823 1883 !-- surface temperature for window surfaces 1824 DO m = 1, surf_usm_h%ns 1825 surf_usm_h%t_surf_window_av(m) = & 1826 surf_usm_h%t_surf_window_av(m) + & 1827 t_surf_window_h(m) 1828 ENDDO 1829 DO l = 0, 3 1884 IF ( l == -1 ) THEN 1885 DO m = 1, surf_usm_h%ns 1886 surf_usm_h%t_surf_window_av(m) = & 1887 surf_usm_h%t_surf_window_av(m) + & 1888 t_surf_window_h(m) 1889 ENDDO 1890 ELSE 1830 1891 DO m = 1, surf_usm_v(l)%ns 1831 1892 surf_usm_v(l)%t_surf_window_av(m) = & … … 1833 1894 t_surf_window_v(l)%t(m) 1834 1895 ENDDO 1835 END DO1896 ENDIF 1836 1897 1837 1898 CASE ( 'usm_t_surf_green' ) 1838 1899 !-- surface temperature for green surfaces 1839 DO m = 1, surf_usm_h%ns 1840 surf_usm_h%t_surf_green_av(m) = & 1841 surf_usm_h%t_surf_green_av(m) + & 1842 t_surf_green_h(m) 1843 ENDDO 1844 DO l = 0, 3 1900 IF ( l == -1 ) THEN 1901 DO m = 1, surf_usm_h%ns 1902 surf_usm_h%t_surf_green_av(m) = & 1903 surf_usm_h%t_surf_green_av(m) + & 1904 t_surf_green_h(m) 1905 ENDDO 1906 ELSE 1845 1907 DO m = 1, surf_usm_v(l)%ns 1846 1908 surf_usm_v(l)%t_surf_green_av(m) = & … … 1848 1910 t_surf_green_v(l)%t(m) 1849 1911 ENDDO 1850 END DO1912 ENDIF 1851 1913 1852 1914 CASE ( 'usm_t_surf_10cm' ) 1853 1915 !-- near surface temperature for whole surfaces 1854 DO m = 1, surf_usm_h%ns 1855 surf_usm_h%t_surf_10cm_av(m) = & 1856 surf_usm_h%t_surf_10cm_av(m) + & 1857 t_surf_10cm_h(m) 1858 ENDDO 1859 DO l = 0, 3 1916 IF ( l == -1 ) THEN 1917 DO m = 1, surf_usm_h%ns 1918 surf_usm_h%t_surf_10cm_av(m) = & 1919 surf_usm_h%t_surf_10cm_av(m) + & 1920 t_surf_10cm_h(m) 1921 ENDDO 1922 ELSE 1860 1923 DO m = 1, surf_usm_v(l)%ns 1861 1924 surf_usm_v(l)%t_surf_10cm_av(m) = & … … 1863 1926 t_surf_10cm_v(l)%t(m) 1864 1927 ENDDO 1865 END DO1928 ENDIF 1866 1929 1867 1930 1868 1931 CASE ( 'usm_t_wall' ) 1869 1932 !-- wall temperature for iwl layer of walls and land 1870 DO m = 1, surf_usm_h%ns 1871 surf_usm_h%t_wall_av(iwl,m) = & 1872 surf_usm_h%t_wall_av(iwl,m) + & 1873 t_wall_h(iwl,m) 1874 ENDDO 1875 DO l = 0, 3 1933 IF ( l == -1 ) THEN 1934 DO m = 1, surf_usm_h%ns 1935 surf_usm_h%t_wall_av(iwl,m) = & 1936 surf_usm_h%t_wall_av(iwl,m) + & 1937 t_wall_h(iwl,m) 1938 ENDDO 1939 ELSE 1876 1940 DO m = 1, surf_usm_v(l)%ns 1877 1941 surf_usm_v(l)%t_wall_av(iwl,m) = & … … 1879 1943 t_wall_v(l)%t(iwl,m) 1880 1944 ENDDO 1881 END DO1945 ENDIF 1882 1946 1883 1947 CASE ( 'usm_t_window' ) 1884 1948 !-- window temperature for iwl layer of walls and land 1885 DO m = 1, surf_usm_h%ns 1886 surf_usm_h%t_window_av(iwl,m) = & 1887 surf_usm_h%t_window_av(iwl,m) + & 1888 t_window_h(iwl,m) 1889 ENDDO 1890 DO l = 0, 3 1949 IF ( l == -1 ) THEN 1950 DO m = 1, surf_usm_h%ns 1951 surf_usm_h%t_window_av(iwl,m) = & 1952 surf_usm_h%t_window_av(iwl,m) + & 1953 t_window_h(iwl,m) 1954 ENDDO 1955 ELSE 1891 1956 DO m = 1, surf_usm_v(l)%ns 1892 1957 surf_usm_v(l)%t_window_av(iwl,m) = & … … 1894 1959 t_window_v(l)%t(iwl,m) 1895 1960 ENDDO 1896 END DO1961 ENDIF 1897 1962 1898 1963 CASE ( 'usm_t_green' ) 1899 1964 !-- green temperature for iwl layer of walls and land 1900 DO m = 1, surf_usm_h%ns 1901 surf_usm_h%t_green_av(iwl,m) = & 1902 surf_usm_h%t_green_av(iwl,m) + & 1903 t_green_h(iwl,m) 1904 ENDDO 1905 DO l = 0, 3 1965 IF ( l == -1 ) THEN 1966 DO m = 1, surf_usm_h%ns 1967 surf_usm_h%t_green_av(iwl,m) = & 1968 surf_usm_h%t_green_av(iwl,m) + & 1969 t_green_h(iwl,m) 1970 ENDDO 1971 ELSE 1906 1972 DO m = 1, surf_usm_v(l)%ns 1907 1973 surf_usm_v(l)%t_green_av(iwl,m) = & … … 1909 1975 t_green_v(l)%t(iwl,m) 1910 1976 ENDDO 1911 END DO1977 ENDIF 1912 1978 1913 1979 CASE DEFAULT … … 1922 1988 CASE ( 'usm_rad_net' ) 1923 1989 !-- array of complete radiation balance 1924 DO m = 1, surf_usm_h%ns 1925 surf_usm_h%rad_net_av(m) = & 1926 surf_usm_h%rad_net_av(m) / & 1927 REAL( average_count_3d, kind=wp ) 1928 ENDDO 1929 DO l = 0, 3 1990 IF ( l == -1 ) THEN 1991 DO m = 1, surf_usm_h%ns 1992 surf_usm_h%rad_net_av(m) = & 1993 surf_usm_h%rad_net_av(m) / & 1994 REAL( average_count_3d, kind=wp ) 1995 ENDDO 1996 ELSE 1930 1997 DO m = 1, surf_usm_v(l)%ns 1931 1998 surf_usm_v(l)%rad_net_av(m) = & … … 1933 2000 REAL( average_count_3d, kind=wp ) 1934 2001 ENDDO 1935 END DO2002 ENDIF 1936 2003 1937 2004 CASE ( 'usm_rad_insw' ) … … 2023 2090 ENDDO 2024 2091 2092 CASE ( 'usm_rad_pc_inlw' ) 2093 pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp ) 2094 2095 CASE ( 'usm_rad_pc_insw' ) 2096 pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp ) 2097 2098 CASE ( 'usm_rad_pc_inswdir' ) 2099 pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp ) 2100 2101 CASE ( 'usm_rad_pc_inswdif' ) 2102 pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp ) 2103 2104 CASE ( 'usm_rad_pc_inswref' ) 2105 pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp ) 2106 2025 2107 CASE ( 'usm_rad_hf' ) 2026 2108 !-- array of heat flux from radiation for surfaces after i-th reflection 2027 DO m = 1, surf_usm_h%ns 2028 surf_usm_h%surfhf_av(m) = & 2029 surf_usm_h%surfhf_av(m) / & 2030 REAL( average_count_3d, kind=wp ) 2031 ENDDO 2032 DO l = 0, 3 2109 IF ( l == -1 ) THEN 2110 DO m = 1, surf_usm_h%ns 2111 surf_usm_h%surfhf_av(m) = & 2112 surf_usm_h%surfhf_av(m) / & 2113 REAL( average_count_3d, kind=wp ) 2114 ENDDO 2115 ELSE 2033 2116 DO m = 1, surf_usm_v(l)%ns 2034 2117 surf_usm_v(l)%surfhf_av(m) = & … … 2036 2119 REAL( average_count_3d, kind=wp ) 2037 2120 ENDDO 2038 END DO2121 ENDIF 2039 2122 2040 2123 CASE ( 'usm_wshf' ) 2041 2124 !-- array of sensible heat flux from surfaces (land, roof, wall) 2042 DO m = 1, surf_usm_h%ns 2043 surf_usm_h%wshf_eb_av(m) = & 2044 surf_usm_h%wshf_eb_av(m) / & 2045 REAL( average_count_3d, kind=wp ) 2046 ENDDO 2047 DO l = 0, 3 2125 IF ( l == -1 ) THEN 2126 DO m = 1, surf_usm_h%ns 2127 surf_usm_h%wshf_eb_av(m) = & 2128 surf_usm_h%wshf_eb_av(m) / & 2129 REAL( average_count_3d, kind=wp ) 2130 ENDDO 2131 ELSE 2048 2132 DO m = 1, surf_usm_v(l)%ns 2049 2133 surf_usm_v(l)%wshf_eb_av(m) = & … … 2051 2135 REAL( average_count_3d, kind=wp ) 2052 2136 ENDDO 2053 END DO2137 ENDIF 2054 2138 2055 2139 CASE ( 'usm_wghf' ) 2056 2140 !-- array of heat flux from ground (wall, roof, land) 2057 DO m = 1, surf_usm_h%ns 2058 surf_usm_h%wghf_eb_av(m) = & 2059 surf_usm_h%wghf_eb_av(m) / & 2060 REAL( average_count_3d, kind=wp ) 2061 ENDDO 2062 DO l = 0, 3 2141 IF ( l == -1 ) THEN 2142 DO m = 1, surf_usm_h%ns 2143 surf_usm_h%wghf_eb_av(m) = & 2144 surf_usm_h%wghf_eb_av(m) / & 2145 REAL( average_count_3d, kind=wp ) 2146 ENDDO 2147 ELSE 2063 2148 DO m = 1, surf_usm_v(l)%ns 2064 2149 surf_usm_v(l)%wghf_eb_av(m) = & … … 2066 2151 REAL( average_count_3d, kind=wp ) 2067 2152 ENDDO 2068 END DO2153 ENDIF 2069 2154 2070 2155 CASE ( 'usm_wghf_window' ) 2071 2156 !-- array of heat flux from window ground (wall, roof, land) 2072 DO m = 1, surf_usm_h%ns 2073 surf_usm_h%wghf_eb_window_av(m) = & 2074 surf_usm_h%wghf_eb_window_av(m) / & 2075 REAL( average_count_3d, kind=wp ) 2076 ENDDO 2077 DO l = 0, 3 2157 IF ( l == -1 ) THEN 2158 DO m = 1, surf_usm_h%ns 2159 surf_usm_h%wghf_eb_window_av(m) = & 2160 surf_usm_h%wghf_eb_window_av(m) / & 2161 REAL( average_count_3d, kind=wp ) 2162 ENDDO 2163 ELSE 2078 2164 DO m = 1, surf_usm_v(l)%ns 2079 2165 surf_usm_v(l)%wghf_eb_window_av(m) = & … … 2081 2167 REAL( average_count_3d, kind=wp ) 2082 2168 ENDDO 2083 END DO2169 ENDIF 2084 2170 2085 2171 CASE ( 'usm_wghf_green' ) 2086 2172 !-- array of heat flux from green ground (wall, roof, land) 2087 DO m = 1, surf_usm_h%ns 2088 surf_usm_h%wghf_eb_green_av(m) = & 2089 surf_usm_h%wghf_eb_green_av(m) / & 2090 REAL( average_count_3d, kind=wp ) 2091 ENDDO 2092 DO l = 0, 3 2173 IF ( l == -1 ) THEN 2174 DO m = 1, surf_usm_h%ns 2175 surf_usm_h%wghf_eb_green_av(m) = & 2176 surf_usm_h%wghf_eb_green_av(m) / & 2177 REAL( average_count_3d, kind=wp ) 2178 ENDDO 2179 ELSE 2093 2180 DO m = 1, surf_usm_v(l)%ns 2094 2181 surf_usm_v(l)%wghf_eb_green_av(m) = & … … 2096 2183 REAL( average_count_3d, kind=wp ) 2097 2184 ENDDO 2098 END DO2185 ENDIF 2099 2186 2100 2187 CASE ( 'usm_iwghf' ) 2101 2188 !-- array of heat flux from indoor ground (wall, roof, land) 2102 DO m = 1, surf_usm_h%ns 2103 surf_usm_h%iwghf_eb_av(m) = & 2104 surf_usm_h%iwghf_eb_av(m) / & 2105 REAL( average_count_3d, kind=wp ) 2106 ENDDO 2107 DO l = 0, 3 2189 IF ( l == -1 ) THEN 2190 DO m = 1, surf_usm_h%ns 2191 surf_usm_h%iwghf_eb_av(m) = & 2192 surf_usm_h%iwghf_eb_av(m) / & 2193 REAL( average_count_3d, kind=wp ) 2194 ENDDO 2195 ELSE 2108 2196 DO m = 1, surf_usm_v(l)%ns 2109 2197 surf_usm_v(l)%iwghf_eb_av(m) = & … … 2111 2199 REAL( average_count_3d, kind=wp ) 2112 2200 ENDDO 2113 END DO2201 ENDIF 2114 2202 2115 2203 CASE ( 'usm_iwghf_window' ) 2116 2204 !-- array of heat flux from indoor window ground (wall, roof, land) 2117 DO m = 1, surf_usm_h%ns 2118 surf_usm_h%iwghf_eb_window_av(m) = & 2119 surf_usm_h%iwghf_eb_window_av(m) / & 2120 REAL( average_count_3d, kind=wp ) 2121 ENDDO 2122 DO l = 0, 3 2205 IF ( l == -1 ) THEN 2206 DO m = 1, surf_usm_h%ns 2207 surf_usm_h%iwghf_eb_window_av(m) = & 2208 surf_usm_h%iwghf_eb_window_av(m) / & 2209 REAL( average_count_3d, kind=wp ) 2210 ENDDO 2211 ELSE 2123 2212 DO m = 1, surf_usm_v(l)%ns 2124 2213 surf_usm_v(l)%iwghf_eb_window_av(m) = & … … 2126 2215 REAL( average_count_3d, kind=wp ) 2127 2216 ENDDO 2128 END DO2217 ENDIF 2129 2218 2130 2219 CASE ( 'usm_t_surf' ) 2131 2220 !-- surface temperature for surfaces 2132 DO m = 1, surf_usm_h%ns 2133 surf_usm_h%t_surf_av(m) = & 2134 surf_usm_h%t_surf_av(m) / & 2135 REAL( average_count_3d, kind=wp ) 2136 ENDDO 2137 DO l = 0, 3 2221 IF ( l == -1 ) THEN 2222 DO m = 1, surf_usm_h%ns 2223 surf_usm_h%t_surf_av(m) = & 2224 surf_usm_h%t_surf_av(m) / & 2225 REAL( average_count_3d, kind=wp ) 2226 ENDDO 2227 ELSE 2138 2228 DO m = 1, surf_usm_v(l)%ns 2139 2229 surf_usm_v(l)%t_surf_av(m) = & … … 2141 2231 REAL( average_count_3d, kind=wp ) 2142 2232 ENDDO 2143 END DO2233 ENDIF 2144 2234 2145 2235 CASE ( 'usm_t_surf_window' ) 2146 2236 !-- surface temperature for window surfaces 2147 DO m = 1, surf_usm_h%ns 2148 surf_usm_h%t_surf_window_av(m) = & 2149 surf_usm_h%t_surf_window_av(m) / & 2150 REAL( average_count_3d, kind=wp ) 2151 ENDDO 2152 DO l = 0, 3 2237 IF ( l == -1 ) THEN 2238 DO m = 1, surf_usm_h%ns 2239 surf_usm_h%t_surf_window_av(m) = & 2240 surf_usm_h%t_surf_window_av(m) / & 2241 REAL( average_count_3d, kind=wp ) 2242 ENDDO 2243 ELSE 2153 2244 DO m = 1, surf_usm_v(l)%ns 2154 2245 surf_usm_v(l)%t_surf_window_av(m) = & … … 2156 2247 REAL( average_count_3d, kind=wp ) 2157 2248 ENDDO 2158 END DO2249 ENDIF 2159 2250 2160 2251 CASE ( 'usm_t_surf_green' ) 2161 2252 !-- surface temperature for green surfaces 2162 DO m = 1, surf_usm_h%ns 2163 surf_usm_h%t_surf_green_av(m) = & 2164 surf_usm_h%t_surf_green_av(m) / & 2165 REAL( average_count_3d, kind=wp ) 2166 ENDDO 2167 DO l = 0, 3 2253 IF ( l == -1 ) THEN 2254 DO m = 1, surf_usm_h%ns 2255 surf_usm_h%t_surf_green_av(m) = & 2256 surf_usm_h%t_surf_green_av(m) / & 2257 REAL( average_count_3d, kind=wp ) 2258 ENDDO 2259 ELSE 2168 2260 DO m = 1, surf_usm_v(l)%ns 2169 2261 surf_usm_v(l)%t_surf_green_av(m) = & … … 2171 2263 REAL( average_count_3d, kind=wp ) 2172 2264 ENDDO 2173 END DO2265 ENDIF 2174 2266 2175 2267 CASE ( 'usm_t_surf_10cm' ) 2176 2268 !-- near surface temperature for whole surfaces 2177 DO m = 1, surf_usm_h%ns 2178 surf_usm_h%t_surf_10cm_av(m) = & 2179 surf_usm_h%t_surf_10cm_av(m) / & 2180 REAL( average_count_3d, kind=wp ) 2181 ENDDO 2182 DO l = 0, 3 2269 IF ( l == -1 ) THEN 2270 DO m = 1, surf_usm_h%ns 2271 surf_usm_h%t_surf_10cm_av(m) = & 2272 surf_usm_h%t_surf_10cm_av(m) / & 2273 REAL( average_count_3d, kind=wp ) 2274 ENDDO 2275 ELSE 2183 2276 DO m = 1, surf_usm_v(l)%ns 2184 2277 surf_usm_v(l)%t_surf_10cm_av(m) = & … … 2186 2279 REAL( average_count_3d, kind=wp ) 2187 2280 ENDDO 2188 END DO2281 ENDIF 2189 2282 2190 2283 CASE ( 'usm_t_wall' ) 2191 2284 !-- wall temperature for iwl layer of walls and land 2192 DO m = 1, surf_usm_h%ns 2193 surf_usm_h%t_wall_av(iwl,m) = & 2194 surf_usm_h%t_wall_av(iwl,m) / & 2195 REAL( average_count_3d, kind=wp ) 2196 ENDDO 2197 DO l = 0, 3 2285 IF ( l == -1 ) THEN 2286 DO m = 1, surf_usm_h%ns 2287 surf_usm_h%t_wall_av(iwl,m) = & 2288 surf_usm_h%t_wall_av(iwl,m) / & 2289 REAL( average_count_3d, kind=wp ) 2290 ENDDO 2291 ELSE 2198 2292 DO m = 1, surf_usm_v(l)%ns 2199 2293 surf_usm_v(l)%t_wall_av(iwl,m) = & … … 2201 2295 REAL( average_count_3d, kind=wp ) 2202 2296 ENDDO 2203 END DO2297 ENDIF 2204 2298 2205 2299 CASE ( 'usm_t_window' ) 2206 2300 !-- window temperature for iwl layer of walls and land 2207 DO m = 1, surf_usm_h%ns 2208 surf_usm_h%t_window_av(iwl,m) = & 2209 surf_usm_h%t_window_av(iwl,m) / & 2210 REAL( average_count_3d, kind=wp ) 2211 ENDDO 2212 DO l = 0, 3 2301 IF ( l == -1 ) THEN 2302 DO m = 1, surf_usm_h%ns 2303 surf_usm_h%t_window_av(iwl,m) = & 2304 surf_usm_h%t_window_av(iwl,m) / & 2305 REAL( average_count_3d, kind=wp ) 2306 ENDDO 2307 ELSE 2213 2308 DO m = 1, surf_usm_v(l)%ns 2214 2309 surf_usm_v(l)%t_window_av(iwl,m) = & … … 2216 2311 REAL( average_count_3d, kind=wp ) 2217 2312 ENDDO 2218 END DO2313 ENDIF 2219 2314 2220 2315 CASE ( 'usm_t_green' ) 2221 2316 !-- green temperature for iwl layer of walls and land 2222 DO m = 1, surf_usm_h%ns 2223 surf_usm_h%t_green_av(iwl,m) = & 2224 surf_usm_h%t_green_av(iwl,m) / & 2225 REAL( average_count_3d, kind=wp ) 2226 ENDDO 2227 DO l = 0, 3 2317 IF ( l == -1 ) THEN 2318 DO m = 1, surf_usm_h%ns 2319 surf_usm_h%t_green_av(iwl,m) = & 2320 surf_usm_h%t_green_av(iwl,m) / & 2321 REAL( average_count_3d, kind=wp ) 2322 ENDDO 2323 ELSE 2228 2324 DO m = 1, surf_usm_v(l)%ns 2229 2325 surf_usm_v(l)%t_green_av(iwl,m) = & … … 2231 2327 REAL( average_count_3d, kind=wp ) 2232 2328 ENDDO 2233 END DO2329 ENDIF 2234 2330 2235 2331 … … 2309 2405 var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. & 2310 2406 var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. & 2311 var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' ) THEN 2407 var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. & 2408 var(1:17) == 'usm_surfwintrans_' ) THEN 2312 2409 unit = 'W/m2' 2313 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. 2410 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. & 2314 2411 var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. & 2315 2412 var(1:16) == 'usm_t_surf_green' .OR. & 2316 2413 var(1:11) == 'usm_t_green' .OR. & 2317 var(1:15) == 'usm_t_surf_10cm' ) THEN2414 var(1:15) == 'usm_t_surf_10cm' ) THEN 2318 2415 unit = 'K' 2416 ELSE IF ( var == 'usm_rad_pc_inlw' .OR. var == 'usm_rad_pc_insw' .OR. & 2417 var == 'usm_rad_pc_inswdir' .OR. var == 'usm_rad_pc_inswdif' .OR. & 2418 var == 'usm_rad_pc_inswref' ) THEN 2419 unit = 'W' 2319 2420 ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. & 2320 2421 var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. & … … 2417 2518 INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart 2418 2519 INTEGER(iwp), DIMENSION(0:nd-1) :: dirend 2419 INTEGER(iwp) :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt 2520 INTEGER(iwp) :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt,ipcgb 2420 2521 INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m 2421 2522 … … 2431 2532 k = len(TRIM(var)) 2432 2533 j = len(TRIM(dirname(i))) 2433 IF ( var(k-j+1:k) == dirname(i) ) THEN2534 IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN 2434 2535 ids = i 2435 2536 idsint = dirint(ids) … … 2802 2903 ENDIF 2803 2904 ENDIF 2905 ENDDO 2906 2907 CASE ( 'usm_rad_pc_inlw' ) 2908 !-- array of lw radiation absorbed by plant canopy 2909 DO ipcgb = 1, npcbl 2910 IF ( av == 0 ) THEN 2911 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw(ipcgb) 2912 ELSE 2913 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw_av(ipcgb) 2914 ENDIF 2915 ENDDO 2916 2917 CASE ( 'usm_rad_pc_insw' ) 2918 !-- array of sw radiation absorbed by plant canopy 2919 DO ipcgb = 1, npcbl 2920 IF ( av == 0 ) THEN 2921 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb) 2922 ELSE 2923 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw_av(ipcgb) 2924 ENDIF 2925 ENDDO 2926 2927 CASE ( 'usm_rad_pc_inswdir' ) 2928 !-- array of direct sw radiation absorbed by plant canopy 2929 DO ipcgb = 1, npcbl 2930 IF ( av == 0 ) THEN 2931 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir(ipcgb) 2932 ELSE 2933 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir_av(ipcgb) 2934 ENDIF 2935 ENDDO 2936 2937 CASE ( 'usm_rad_pc_inswdif' ) 2938 !-- array of diffuse sw radiation absorbed by plant canopy 2939 DO ipcgb = 1, npcbl 2940 IF ( av == 0 ) THEN 2941 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif(ipcgb) 2942 ELSE 2943 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif_av(ipcgb) 2944 ENDIF 2945 ENDDO 2946 2947 CASE ( 'usm_rad_pc_inswref' ) 2948 !-- array of reflected sw radiation absorbed by plant canopy 2949 DO ipcgb = 1, npcbl 2950 IF ( av == 0 ) THEN 2951 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb) & 2952 - pcbinswdir(ipcgb) & 2953 - pcbinswdif(ipcgb) 2954 ELSE 2955 temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswref_av(ipcgb) 2956 ENDIF 2804 2957 ENDDO 2805 2958 … … 3356 3509 CASE DEFAULT 3357 3510 found = .FALSE. 3358 3511 RETURN 3359 3512 END SELECT 3360 3513 … … 3399 3552 var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. & 3400 3553 var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. & 3401 var(1:11) == 'usm_rad_hf_' .OR. & 3554 var(1:11) == 'usm_rad_hf_' .OR. var == 'usm_rad_pc_inlw' .OR. & 3555 var == 'usm_rad_pc_insw' .OR. var == 'usm_rad_pc_inswdir' .OR. & 3556 var == 'usm_rad_pc_inswdif' .OR. var == 'usm_rad_pc_inswref' .OR. & 3402 3557 var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. & 3403 3558 var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. & … … 4965 5120 * wintend(nzb_wall:nzt_wall) + tsc(3) & 4966 5121 * surf_usm_h%tt_window_m(nzb_wall:nzt_wall,m) ) 4967 5122 4968 5123 ! 4969 5124 !-- calculate t_wall tendencies for the next Runge-Kutta step … … 5168 5323 ENDIF 5169 5324 ENDDO 5325 !!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!! 5326 ! t_window_v_p(l)%t = t_wall_v_p(l)%t 5327 ! surf_usm_v(l)%tt_window_m = surf_usm_v(l)%tt_wall_m 5328 ! t_green_v_p(l)%t = t_wall_v_p(l)%t 5329 ! surf_usm_v(l)%tt_green_m = surf_usm_v(l)%tt_wall_m 5330 !!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!! 5170 5331 ENDDO 5171 5332 … … 5320 5481 naheatlayers, & 5321 5482 pedestrian_category, & 5483 roughness_concrete, & 5322 5484 read_wall_temp_3d, & 5323 5485 roof_category, & … … 5337 5499 naheatlayers, & 5338 5500 pedestrian_category, & 5501 roughness_concrete, & 5339 5502 read_wall_temp_3d, & 5340 5503 roof_category, & … … 6937 7100 ENDDO 6938 7101 IF ( ip == -99999 ) THEN 6939 !-- wall category not found 6940 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, & 6941 ' not found for i,j,k=', iw,jw,kw 6942 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 ) 7102 !-- land/roof category not found 7103 WRITE (9,"(A,I5,A,3I5)") 'land/roof category ', it, & 7104 ' not found for i,j,k=', iw,jw,kw 7105 FLUSH(9) 7106 IF ( surf_usm_h%isroof_surf(m) ) THEN 7107 category = roof_category 7108 ELSE 7109 category = land_category 7110 ENDIF 7111 DO k = 1, n_surface_types 7112 IF ( surface_type_codes(k) == roof_category ) THEN 7113 ip = k 7114 EXIT 7115 ENDIF 7116 ENDDO 7117 IF ( ip == -99999 ) THEN 7118 !-- default land/roof category not found 7119 WRITE (9,"(A,I5,A,3I5)") 'Default land/roof category', category, ' not found!' 7120 FLUSH(9) 7121 ip = 1 7122 ENDIF 6943 7123 ENDIF 6944 7124 ! … … 7036 7216 surf_usm_v(l)%albedo(:,m) = -1.0_wp 7037 7217 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 7218 surf_usm_v(l)%thickness_window(m) = -1.0_wp 7219 surf_usm_v(l)%thickness_green(m) = -1.0_wp 7220 surf_usm_v(l)%transmissivity(m) = -1.0_wp 7038 7221 ELSE IF ( kw <= usm_par(ii,jw,iw) ) THEN 7039 7222 !-- pedestrian zone … … 7089 7272 ELSE 7090 7273 ! 7274 WRITE(9,*) 'Problem reading USM data:' 7275 WRITE(9,*) l,i,j,kw,get_topography_top_index_ji( j, i, 's' ) 7276 WRITE(9,*) ii,iw,jw,kw,get_topography_top_index_ji( jw, iw, 's' ) 7277 WRITE(9,*) usm_par(ii,jw,iw),usm_par(ii+1,jw,iw) 7278 WRITE(9,*) usm_par(ii+2,jw,iw),usm_par(ii+3,jw,iw) 7279 WRITE(9,*) usm_par(ii+4,jw,iw),usm_par(ii+5,jw,iw) 7280 WRITE(9,*) kw,roof_height_limit,wall_category,roof_category 7281 FLUSH(9) 7091 7282 !-- supply the default category 7092 7283 IF ( kw <= roof_height_limit ) THEN … … 7097 7288 surf_usm_v(l)%albedo(:,m) = -1.0_wp 7098 7289 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 7290 surf_usm_v(l)%thickness_window(m) = -1.0_wp 7291 surf_usm_v(l)%thickness_green(m) = -1.0_wp 7292 surf_usm_v(l)%transmissivity(m) = -1.0_wp 7099 7293 ENDIF 7100 7294 ! … … 7110 7304 IF ( ip == -99999 ) THEN 7111 7305 !-- wall category not found 7112 WRITE (message_string, "(A,I7,A,3I5)") 'wall category ', it, & 7113 ' not found for i,j,k=', iw,jw,kw 7114 WRITE(9,*) message_string 7306 WRITE (9, "(A,I7,A,3I5)") 'wall category ', it, & 7307 ' not found for i,j,k=', iw,jw,kw 7308 FLUSH(9) 7309 category = wall_category 7310 DO k = 1, n_surface_types 7311 IF ( surface_type_codes(k) == category ) THEN 7312 ip = k 7313 EXIT 7314 ENDIF 7315 ENDDO 7316 IF ( ip == -99999 ) THEN 7317 !-- default wall category not found 7318 WRITE (9, "(A,I5,A,3I5)") 'Default wall category', category, ' not found!' 7319 FLUSH(9) 7320 ip = 1 7321 ENDIF 7115 7322 ENDIF 7323 7116 7324 ! 7117 7325 !-- Albedo … … 7192 7400 ENDDO 7193 7401 7402 7403 WRITE(9,*) 'Urban surfaces read' 7404 FLUSH(9) 7405 7194 7406 CALL location_message( ' types and parameters of urban surfaces read', .TRUE. ) 7195 7407 … … 7473 7685 IF ( humidity ) surf_usm_h%vpt_surface(m) = & 7474 7686 surf_usm_h%pt_surface(m) 7475 7687 7476 7688 !-- calculate true tendency 7477 7689 stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) * & … … 7508 7720 ! 7509 7721 !-- calculate fluxes 7510 !-- rad_net_l is never used! 7722 !-- rad_net_l is never used! 7511 7723 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) + & 7512 7724 surf_usm_h%frac(ind_veg_wall,m) * & … … 7522 7734 surf_usm_h%wghf_eb(m) = lambda_surface * & 7523 7735 ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) ) 7524 surf_usm_h%wghf_eb_green(m) = lambda_surface_green * &7736 surf_usm_h%wghf_eb_green(m) = lambda_surface_green * & 7525 7737 ( t_surf_green_h_p(m) - t_green_h(nzb_wall,m) ) 7526 surf_usm_h%wghf_eb_window(m) = lambda_surface_window * &7738 surf_usm_h%wghf_eb_window(m) = lambda_surface_window * & 7527 7739 ( t_surf_window_h_p(m) - t_window_h(nzb_wall,m) ) 7528 7740 … … 7539 7751 !-- diffusion_s, surface_layer_fluxes,... 7540 7752 surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p 7541 7753 7542 7754 ENDDO 7543 7755 ! … … 7562 7774 #if ! defined( __nopointer ) 7563 7775 ! 7564 !-- calculate rho * c_p coefficient at surfacelayer7776 !-- calculate rho * c_p coefficient at wall layer 7565 7777 rho_cp = c_p * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exner(k) ) 7566 7778 #endif … … 7589 7801 !-- obtained by simple linear interpolation. ( An alternative would 7590 7802 !-- be an logarithmic interpolation. ) 7591 !-- A roughness lenght of 0.001 is assumed for concrete (the inverse,7592 !-- 1000 is used in the nominator for scaling)7593 surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) * 1000.0_wp&7594 * ( 11.8_wp + 4.2_wp *&7803 !-- Parameter roughness_concrete (default value = 0.001) is used 7804 !-- to calculation of roughness relative to concrete 7805 surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) / & 7806 roughness_concrete * ( 11.8_wp + 4.2_wp * & 7595 7807 SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + & 7596 7808 ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + & … … 7606 7818 f_shf_window = rho_cp / surf_usm_v(l)%r_a(m) 7607 7819 f_shf_green = rho_cp / surf_usm_v(l)%r_a(m) 7608 7609 7610 7820 7611 7821 !-- add LW up so that it can be removed in prognostic equation … … 7657 7867 ( surf_usm_v(l)%c_surface_green(m) + coef_green_2 * dt_3d * tsc(2) ) 7658 7868 7659 7660 7661 7869 !-- add RK3 term 7662 7870 t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) * & … … 7677 7885 IF ( humidity ) surf_usm_v(l)%vpt_surface(m) = & 7678 7886 surf_usm_v(l)%pt_surface(m) 7679 7887 7680 7888 !-- calculate true tendency 7681 7889 stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *& … … 7766 7974 dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp) 7767 7975 dhour = INT(dtime/3600.0_wp) 7768 DO m = 1, naheatlayers 7769 !-- Get indices of respective grid point 7770 i = surf_usm_h%i(m) 7771 j = surf_usm_h%j(m) 7772 k = surf_usm_h%k(m) 7773 IF ( k > get_topography_top_index_ji( j, i, 's' ) .AND. & 7774 k <= naheatlayers ) THEN 7775 !-- increase of pt in box i,j,k in time dt_3d 7776 !-- given to anthropogenic heat aheat*acoef (W*m-2) 7777 !-- linear interpolation of coeficient 7778 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + & 7779 (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1) 7780 IF ( aheat(k,j,i) > 0.0_wp ) THEN 7781 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dzu(k)) 7782 ENDIF 7783 ENDIF 7976 DO i = nxl, nxr 7977 DO j = nys, nyn 7978 DO k = nzub, min(nzut,naheatlayers) 7979 IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN 7980 !-- increase of pt in box i,j,k in time dt_3d 7981 !-- given to anthropogenic heat aheat*acoef (W*m-2) 7982 !-- linear interpolation of coeficient 7983 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + & 7984 (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1) 7985 IF ( aheat(k,j,i) > 0.0_wp ) THEN 7986 !-- calculate rho * c_p coefficient at layer k 7987 rho_cp = c_p * hyp(k) / ( r_d * pt(k+1,j,i) * exner(k) ) 7988 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dz(1)) 7989 ENDIF 7990 ENDIF 7991 ENDDO 7992 ENDDO 7784 7993 ENDDO 7785 7994
Note: See TracChangeset
for help on using the changeset viewer.
