Changeset 2920 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- Mar 22, 2018 11:22:01 AM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/urban_surface_mod.f90
r2906 r2920 16 16 ! 17 17 ! Copyright 2015-2018 Czech Technical University in Prague 18 ! Copyright 2015-2018 Institute of Computer Science of the 19 ! Czech Academy of Sciences, Prague 18 20 ! Copyright 1997-2018 Leibniz Universitaet Hannover 19 21 !------------------------------------------------------------------------------! … … 26 28 ! ----------------- 27 29 ! $Id$ 30 ! Remove unused pcbl, npcbl from ONLY list 31 ! moh.hefny: 32 ! Fixed bugs introduced by new structures and by moving radiation interaction 33 ! into radiation_model_mod.f90. 34 ! Bugfix: usm data output 3D didn't respect directions 35 ! 36 ! 2906 2018-03-19 08:56:40Z Giersch 28 37 ! Local variable ids has to be initialized with a value of -1 in 29 38 ! usm_average_3d_data … … 217 226 !> Further work: 218 227 !> ------------- 219 !> 1. Reduce number of shape view factors by merging factors for distant surfaces 220 !> under shallow angles. Idea: Iteratively select the smallest shape view 221 !> factor by value (among all sources and targets) which has a similarly 222 !> oriented source neighbor (or near enough) SVF and merge them by adding 223 !> value of the smaller SVF to the larger one and deleting the smaller one. 224 !> This will allow for better scaling at higher resolutions. 225 !> 226 !> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity 228 !> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity 227 229 !> from surfaces that are visible from local surfaces (i.e. there is a SVF 228 230 !> where target is local). To do that, radiosity will be exchanged after each 229 231 !> reflection step using MPI_Alltoall instead of current MPI_Allgather. 230 232 !> 231 !> 3. Temporarily large values of surface heat flux can be observed, up to233 !> 2. Temporarily large values of surface heat flux can be observed, up to 232 234 !> 1.2 Km/s, which seem to be not realistic. 233 235 !> … … 248 250 #if ! defined( __nopointer ) 249 251 USE arrays_3d, & 250 ONLY: zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend252 ONLY: hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend 251 253 #endif 252 254 … … 272 274 273 275 USE date_and_time_mod, & 274 ONLY: d_seconds_year, day_of_year_init,time_utc_init276 ONLY: time_utc_init 275 277 276 278 USE grid_variables, & … … 288 290 289 291 USE plant_canopy_model_mod, & 290 ONLY: pc_heating_rate , usm_lad_rma292 ONLY: pc_heating_rate 291 293 292 294 USE radiation_model_mod, & 293 ONLY: albedo_type, radiation , calc_zenith, zenith, &295 ONLY: albedo_type, radiation_interaction, calc_zenith, zenith, & 294 296 rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, & 295 297 sigma_sb, solar_constant, sun_direction, sun_dir_lat, & … … 297 299 force_radiation_call, surfinsw, surfinlw, surfinswdir, & 298 300 surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, & 299 surfinl, surfinlwdif, energy_balance_surf_h, & 300 energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff, & 301 surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff, & 301 302 rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, & 302 surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy, & 303 endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, & 304 inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id, & 305 iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut, & 306 isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf, & 307 iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf, & 308 idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall 303 surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, & 304 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, & 305 inorth_l, isouth_l, ieast_l, iwest_l, id, & 306 iz, iy, ix, idir, jdir, kdir, nsurf_type, nzub, nzut, & 307 nzu, pch, nsurf, idsvf, ndsvf, & 308 iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a, & 309 idcsf, ndcsf, kdcsf, pct, & 310 startland, endland, startwall, endwall, skyvf, skyvft 309 311 310 312 USE statistics, & … … 322 324 LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls 323 325 LOGICAL :: indoor_model = .FALSE. !< whether to use the indoor model 324 325 326 LOGICAL :: read_wall_temp_3d = .FALSE. 327 328 326 329 INTEGER(iwp) :: building_type = 1 !< default building type (preleminary setting) 327 330 INTEGER(iwp) :: land_category = 2 !< default category for land surface 328 331 INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone 329 INTEGER(iwp) :: pedestr ant_category = 2 !< default category for wall surface in pedestrian zone332 INTEGER(iwp) :: pedestrian_category = 2 !< default category for wall surface in pedestrian zone 330 333 INTEGER(iwp) :: roof_category = 2 !< default category for root surface 331 334 REAL(wp) :: roughness_concrete = 0.001_wp !< roughness length of average concrete surface 332 335 ! 333 336 !-- Indices of input attributes for (above) ground floor level … … 380 383 REAL(wp) :: roof_height_limit = 4._wp !< height for distinguish between land surfaces and roofs 381 384 REAL(wp) :: ground_floor_level = 4.0_wp !< default ground floor level 382 REAL(wp) :: ra_horiz_coef = 5.0_wp !< mysterious coefficient for correction of overestimation 383 !< of r_a for horizontal surfaces -> TODO 385 384 386 385 387 CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/ & … … 524 526 !-- anthropogenic heat sources 525 527 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 526 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2) 527 REAL(wp), DIMENSION(:), ALLOCATABLE :: aheatprof !< diurnal profile of anthropogenic heat 528 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2) 529 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheatprof !< diurnal profiles of anthropogenic heat for particular layers 530 INTEGER(wp) :: naheatlayers = 1 !< number of layers of anthropogenic heat 528 531 529 532 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 547 550 ! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /) 548 551 549 550 REAL(wp) :: wall_inner_temperature = 296.0_wp !< temperature of the inner wall surface (~23degrees C) (K)551 REAL(wp) :: roof_inner_temperature = 296.0_wp !< temperature of the inner roof surface (~23degrees C) (K)552 REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10degrees C) (K)553 REAL(wp) :: window_inner_temperature = 296.0_wp !< temperature of the inner window surface (~23degrees C) (K)552 553 REAL(wp) :: wall_inner_temperature = 295.0_wp !< temperature of the inner wall surface (~22 degrees C) (K) 554 REAL(wp) :: roof_inner_temperature = 295.0_wp !< temperature of the inner roof surface (~22 degrees C) (K) 555 REAL(wp) :: soil_inner_temperature = 288.0_wp !< temperature of the deep soil (~15 degrees C) (K) 556 REAL(wp) :: window_inner_temperature = 295.0_wp !< temperature of the inner window surface (~22 degrees C) (K) 554 557 555 558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 674 677 !-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity 675 678 INTEGER(iwp) :: n_surface_types !< number of the wall type categories 676 INTEGER(iwp), PARAMETER :: n_surface_params = 8!< number of parameters for each type of the wall679 INTEGER(iwp), PARAMETER :: n_surface_params = 9 !< number of parameters for each type of the wall 677 680 INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface 678 681 INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface 679 INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity λS between air and surface ( W mâ2 Kâ1 ) 680 INTEGER(iwp), PARAMETER :: irough = 4 !< roughness relative to concrete 681 INTEGER(iwp), PARAMETER :: icsurf = 5 !< Surface skin layer heat capacity (J mâ2 Kâ1 ) 682 INTEGER(iwp), PARAMETER :: ithick = 6 !< thickness of the surface (wall, roof, land) ( m ) 683 INTEGER(iwp), PARAMETER :: irhoC = 7 !< volumetric heat capacity rho*C of the material ( J mâ3 Kâ1 ) 684 INTEGER(iwp), PARAMETER :: ilambdah = 8 !< thermal conductivity λH of the wall (W mâ1 Kâ1 ) 682 INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity lambda S between surface and material ( W m-2 K-1 ) 683 INTEGER(iwp), PARAMETER :: irough = 4 !< roughness length z0 for movements 684 INTEGER(iwp), PARAMETER :: iroughh = 5 !< roughness length z0h for scalars (heat, humidity,...) 685 INTEGER(iwp), PARAMETER :: icsurf = 6 !< Surface skin layer heat capacity (J m-2 K-1 ) 686 INTEGER(iwp), PARAMETER :: ithick = 7 !< thickness of the surface (wall, roof, land) ( m ) 687 INTEGER(iwp), PARAMETER :: irhoC = 8 !< volumetric heat capacity rho*C of the material ( J m-3 K-1 ) 688 INTEGER(iwp), PARAMETER :: ilambdah = 9 !< thermal conductivity lambda H of the wall (W m-1 K-1 ) 685 689 CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports) 686 690 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types … … 760 764 !-- Public functions 761 765 PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,& 762 usm_rrd_local, 763 usm_surface_energy_balance, usm_material_heat_model, 764 usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, 765 usm_data_output_3d, usm_define_netcdf_grid, usm_parin, 766 usm_rrd_local, & 767 usm_surface_energy_balance, usm_material_heat_model, & 768 usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, & 769 usm_data_output_3d, usm_define_netcdf_grid, usm_parin, & 766 770 usm_wrd_local, usm_allocate_surface 767 771 768 772 !-- Public parameters, constants and initial values 769 PUBLIC usm_anthropogenic_heat, usm_material_model, ra_horiz_coef,&773 PUBLIC usm_anthropogenic_heat, usm_material_model, & 770 774 usm_green_heat_model, usm_temperature_near_surface 771 775 … … 1175 1179 CHARACTER (len=*), INTENT(IN) :: variable 1176 1180 1177 INTEGER(iwp) :: i, j, k, l, m, ids, i wl,istat1181 INTEGER(iwp) :: i, j, k, l, m, ids, idsint, iwl, istat 1178 1182 CHARACTER (len=varnamelength) :: var, surfid 1179 1183 INTEGER(iwp), PARAMETER :: nd = 5 1180 1184 CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 1185 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /) 1181 1186 1182 1187 !-- find the real name of the variable … … 1188 1193 IF ( var(k-j+1:k) == dirname(i) ) THEN 1189 1194 ids = i 1195 idsint = dirint(ids) 1190 1196 var = var(:k-j) 1191 1197 EXIT … … 1245 1251 CASE ( 'usm_rad_insw' ) 1246 1252 !-- array of sw radiation falling to surface after i-th reflection 1247 IF ( .NOT. ALLOCATED(surf _usm_h%surfinsw_av) ) THEN1248 ALLOCATE( surf _usm_h%surfinsw_av(1:surf_usm_h%ns) )1249 surf _usm_h%surfinsw_av = 0.0_wp1253 IF ( .NOT. ALLOCATED(surfinsw_av) ) THEN 1254 ALLOCATE( surfinsw_av(nsurfl) ) 1255 surfinsw_av = 0.0_wp 1250 1256 ENDIF 1251 DO l = 0, 3 1252 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinsw_av) ) THEN 1253 ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) ) 1254 surf_usm_v(l)%surfinsw_av = 0.0_wp 1255 ENDIF 1256 ENDDO 1257 1257 1258 1258 CASE ( 'usm_rad_inlw' ) 1259 1259 !-- array of lw radiation falling to surface after i-th reflection 1260 IF ( .NOT. ALLOCATED(surf _usm_h%surfinlw_av) ) THEN1261 ALLOCATE( surf _usm_h%surfinlw_av(1:surf_usm_h%ns) )1262 surf _usm_h%surfinlw_av = 0.0_wp1260 IF ( .NOT. ALLOCATED(surfinlw_av) ) THEN 1261 ALLOCATE( surfinlw_av(nsurfl) ) 1262 surfinlw_av = 0.0_wp 1263 1263 ENDIF 1264 DO l = 0, 31265 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinlw_av) ) THEN1266 ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )1267 surf_usm_v(l)%surfinlw_av = 0.0_wp1268 ENDIF1269 ENDDO1270 1264 1271 1265 CASE ( 'usm_rad_inswdir' ) 1272 1266 !-- array of direct sw radiation falling to surface from sun 1273 1267 IF ( .NOT. ALLOCATED(surfinswdir_av) ) THEN 1274 ALLOCATE( surfinswdir_av( startenergy:endenergy) )1268 ALLOCATE( surfinswdir_av(nsurfl) ) 1275 1269 surfinswdir_av = 0.0_wp 1276 1270 ENDIF … … 1279 1273 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 1280 1274 IF ( .NOT. ALLOCATED(surfinswdif_av) ) THEN 1281 ALLOCATE( surfinswdif_av( startenergy:endenergy) )1275 ALLOCATE( surfinswdif_av(nsurfl) ) 1282 1276 surfinswdif_av = 0.0_wp 1283 1277 ENDIF … … 1286 1280 !-- array of sw radiation falling to surface from reflections 1287 1281 IF ( .NOT. ALLOCATED(surfinswref_av) ) THEN 1288 ALLOCATE( surfinswref_av( startenergy:endenergy) )1282 ALLOCATE( surfinswref_av(nsurfl) ) 1289 1283 surfinswref_av = 0.0_wp 1290 1284 ENDIF … … 1293 1287 !-- array of sw radiation falling to surface after i-th reflection 1294 1288 IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN 1295 ALLOCATE( surfinlwdif_av( startenergy:endenergy) )1289 ALLOCATE( surfinlwdif_av(nsurfl) ) 1296 1290 surfinlwdif_av = 0.0_wp 1297 1291 ENDIF … … 1300 1294 !-- array of lw radiation falling to surface from reflections 1301 1295 IF ( .NOT. ALLOCATED(surfinlwref_av) ) THEN 1302 ALLOCATE( surfinlwref_av( startenergy:endenergy) )1296 ALLOCATE( surfinlwref_av(nsurfl) ) 1303 1297 surfinlwref_av = 0.0_wp 1304 1298 ENDIF … … 1307 1301 !-- array of sw radiation emitted from surface after i-th reflection 1308 1302 IF ( .NOT. ALLOCATED(surfoutsw_av) ) THEN 1309 ALLOCATE( surfoutsw_av( startenergy:endenergy) )1303 ALLOCATE( surfoutsw_av(nsurfl) ) 1310 1304 surfoutsw_av = 0.0_wp 1311 1305 ENDIF … … 1314 1308 !-- array of lw radiation emitted from surface after i-th reflection 1315 1309 IF ( .NOT. ALLOCATED(surfoutlw_av) ) THEN 1316 ALLOCATE( surfoutlw_av( startenergy:endenergy) )1310 ALLOCATE( surfoutlw_av(nsurfl) ) 1317 1311 surfoutlw_av = 0.0_wp 1318 1312 ENDIF … … 1320 1314 !-- array of residua of sw radiation absorbed in surface after last reflection 1321 1315 IF ( .NOT. ALLOCATED(surfins_av) ) THEN 1322 ALLOCATE( surfins_av( startenergy:endenergy) )1316 ALLOCATE( surfins_av(nsurfl) ) 1323 1317 surfins_av = 0.0_wp 1324 1318 ENDIF … … 1327 1321 !-- array of residua of lw radiation absorbed in surface after last reflection 1328 1322 IF ( .NOT. ALLOCATED(surfinl_av) ) THEN 1329 ALLOCATE( surfinl_av( startenergy:endenergy) )1323 ALLOCATE( surfinl_av(nsurfl) ) 1330 1324 surfinl_av = 0.0_wp 1331 1325 ENDIF … … 1544 1538 CASE ( 'usm_rad_insw' ) 1545 1539 !-- array of sw radiation falling to surface after i-th reflection 1546 DO l = startenergy, endenergy1547 IF ( surfl(id,l) == ids ) THEN1540 DO l = 1, nsurfl 1541 IF ( surfl(id,l) == idsint ) THEN 1548 1542 surfinsw_av(l) = surfinsw_av(l) + surfinsw(l) 1549 1543 ENDIF … … 1552 1546 CASE ( 'usm_rad_inlw' ) 1553 1547 !-- array of lw radiation falling to surface after i-th reflection 1554 DO l = startenergy, endenergy1555 IF ( surfl(id,l) == ids ) THEN1548 DO l = 1, nsurfl 1549 IF ( surfl(id,l) == idsint ) THEN 1556 1550 surfinlw_av(l) = surfinlw_av(l) + surfinlw(l) 1557 1551 ENDIF … … 1560 1554 CASE ( 'usm_rad_inswdir' ) 1561 1555 !-- array of direct sw radiation falling to surface from sun 1562 DO l = startenergy, endenergy1563 IF ( surfl(id,l) == ids ) THEN1556 DO l = 1, nsurfl 1557 IF ( surfl(id,l) == idsint ) THEN 1564 1558 surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l) 1565 1559 ENDIF … … 1568 1562 CASE ( 'usm_rad_inswdif' ) 1569 1563 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 1570 DO l = startenergy, endenergy1571 IF ( surfl(id,l) == ids ) THEN1564 DO l = 1, nsurfl 1565 IF ( surfl(id,l) == idsint ) THEN 1572 1566 surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l) 1573 1567 ENDIF … … 1576 1570 CASE ( 'usm_rad_inswref' ) 1577 1571 !-- array of sw radiation falling to surface from reflections 1578 DO l = startenergy, endenergy1579 IF ( surfl(id,l) == ids ) THEN1572 DO l = 1, nsurfl 1573 IF ( surfl(id,l) == idsint ) THEN 1580 1574 surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - & 1581 1575 surfinswdir(l) - surfinswdif(l) … … 1586 1580 CASE ( 'usm_rad_inlwdif' ) 1587 1581 !-- array of sw radiation falling to surface after i-th reflection 1588 DO l = startenergy, endenergy 1589 IF ( surfl(id,l) == ids ) THEN 1590 surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - & 1591 surfinswdir(l) - surfinswdif(l) 1582 DO l = 1, nsurfl 1583 IF ( surfl(id,l) == idsint ) THEN 1584 surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l) 1592 1585 ENDIF 1593 1586 ENDDO … … 1595 1588 CASE ( 'usm_rad_inlwref' ) 1596 1589 !-- array of lw radiation falling to surface from reflections 1597 DO l = startenergy, endenergy 1598 IF ( surfl(id,l) == ids ) THEN 1599 surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l) 1590 DO l = 1, nsurfl 1591 IF ( surfl(id,l) == idsint ) THEN 1592 surfinlwref_av(l) = surfinlwref_av(l) + & 1593 surfinlw(l) - surfinlwdif(l) 1600 1594 ENDIF 1601 1595 ENDDO … … 1603 1597 CASE ( 'usm_rad_outsw' ) 1604 1598 !-- array of sw radiation emitted from surface after i-th reflection 1605 DO l = startenergy, endenergy 1606 IF ( surfl(id,l) == ids ) THEN 1607 surfinlwref_av(l) = surfinlwref_av(l) + & 1608 surfinlw(l) - surfinlwdif(l) 1599 DO l = 1, nsurfl 1600 IF ( surfl(id,l) == idsint ) THEN 1601 surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l) 1609 1602 ENDIF 1610 1603 ENDDO … … 1612 1605 CASE ( 'usm_rad_outlw' ) 1613 1606 !-- array of lw radiation emitted from surface after i-th reflection 1614 DO l = startenergy, endenergy1615 IF ( surfl(id,l) == ids ) THEN1616 surfout sw_av(l) = surfoutsw_av(l) + surfoutsw(l)1607 DO l = 1, nsurfl 1608 IF ( surfl(id,l) == idsint ) THEN 1609 surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l) 1617 1610 ENDIF 1618 1611 ENDDO … … 1620 1613 CASE ( 'usm_rad_ressw' ) 1621 1614 !-- array of residua of sw radiation absorbed in surface after last reflection 1622 DO l = startenergy, endenergy1623 IF ( surfl(id,l) == ids ) THEN1624 surf outlw_av(l) = surfoutlw_av(l) + surfoutlw(l)1615 DO l = 1, nsurfl 1616 IF ( surfl(id,l) == idsint ) THEN 1617 surfins_av(l) = surfins_av(l) + surfins(l) 1625 1618 ENDIF 1626 1619 ENDDO … … 1628 1621 CASE ( 'usm_rad_reslw' ) 1629 1622 !-- array of residua of lw radiation absorbed in surface after last reflection 1630 DO l = startenergy, endenergy1631 IF ( surfl(id,l) == ids ) THEN1632 surfin s_av(l) = surfins_av(l) + surfins(l)1623 DO l = 1, nsurfl 1624 IF ( surfl(id,l) == idsint ) THEN 1625 surfinl_av(l) = surfinl_av(l) + surfinl(l) 1633 1626 ENDIF 1634 1627 ENDDO … … 1871 1864 CASE ( 'usm_rad_insw' ) 1872 1865 !-- array of sw radiation falling to surface after i-th reflection 1873 DO l = startenergy, endenergy1874 IF ( surfl(id,l) == ids ) THEN1866 DO l = 1, nsurfl 1867 IF ( surfl(id,l) == idsint ) THEN 1875 1868 surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp ) 1876 1869 ENDIF … … 1879 1872 CASE ( 'usm_rad_inlw' ) 1880 1873 !-- array of lw radiation falling to surface after i-th reflection 1881 DO l = startenergy, endenergy1882 IF ( surfl(id,l) == ids ) THEN1874 DO l = 1, nsurfl 1875 IF ( surfl(id,l) == idsint ) THEN 1883 1876 surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp ) 1884 1877 ENDIF … … 1887 1880 CASE ( 'usm_rad_inswdir' ) 1888 1881 !-- array of direct sw radiation falling to surface from sun 1889 DO l = startenergy, endenergy1890 IF ( surfl(id,l) == ids ) THEN1882 DO l = 1, nsurfl 1883 IF ( surfl(id,l) == idsint ) THEN 1891 1884 surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp ) 1892 1885 ENDIF … … 1895 1888 CASE ( 'usm_rad_inswdif' ) 1896 1889 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 1897 DO l = startenergy, endenergy1898 IF ( surfl(id,l) == ids ) THEN1890 DO l = 1, nsurfl 1891 IF ( surfl(id,l) == idsint ) THEN 1899 1892 surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp ) 1900 1893 ENDIF … … 1903 1896 CASE ( 'usm_rad_inswref' ) 1904 1897 !-- array of sw radiation falling to surface from reflections 1905 DO l = startenergy, endenergy1906 IF ( surfl(id,l) == ids ) THEN1898 DO l = 1, nsurfl 1899 IF ( surfl(id,l) == idsint ) THEN 1907 1900 surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp ) 1908 1901 ENDIF … … 1911 1904 CASE ( 'usm_rad_inlwdif' ) 1912 1905 !-- array of sw radiation falling to surface after i-th reflection 1913 DO l = startenergy, endenergy1914 IF ( surfl(id,l) == ids ) THEN1906 DO l = 1, nsurfl 1907 IF ( surfl(id,l) == idsint ) THEN 1915 1908 surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp ) 1916 1909 ENDIF … … 1919 1912 CASE ( 'usm_rad_inlwref' ) 1920 1913 !-- array of lw radiation falling to surface from reflections 1921 DO l = startenergy, endenergy1922 IF ( surfl(id,l) == ids ) THEN1914 DO l = 1, nsurfl 1915 IF ( surfl(id,l) == idsint ) THEN 1923 1916 surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp ) 1924 1917 ENDIF … … 1927 1920 CASE ( 'usm_rad_outsw' ) 1928 1921 !-- array of sw radiation emitted from surface after i-th reflection 1929 DO l = startenergy, endenergy1930 IF ( surfl(id,l) == ids ) THEN1922 DO l = 1, nsurfl 1923 IF ( surfl(id,l) == idsint ) THEN 1931 1924 surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp ) 1932 1925 ENDIF … … 1935 1928 CASE ( 'usm_rad_outlw' ) 1936 1929 !-- array of lw radiation emitted from surface after i-th reflection 1937 DO l = startenergy, endenergy1938 IF ( surfl(id,l) == ids ) THEN1930 DO l = 1, nsurfl 1931 IF ( surfl(id,l) == idsint ) THEN 1939 1932 surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp ) 1940 1933 ENDIF … … 1943 1936 CASE ( 'usm_rad_ressw' ) 1944 1937 !-- array of residua of sw radiation absorbed in surface after last reflection 1945 DO l = startenergy, endenergy1946 IF ( surfl(id,l) == ids ) THEN1938 DO l = 1, nsurfl 1939 IF ( surfl(id,l) == idsint ) THEN 1947 1940 surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp ) 1948 1941 ENDIF … … 1951 1944 CASE ( 'usm_rad_reslw' ) 1952 1945 !-- array of residua of lw radiation absorbed in surface after last reflection 1953 DO l = startenergy, endenergy1954 IF ( surfl(id,l) == ids ) THEN1946 DO l = 1, nsurfl 1947 IF ( surfl(id,l) == idsint ) THEN 1955 1948 surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp ) 1956 1949 ENDIF … … 2245 2238 var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' ) THEN 2246 2239 unit = 'W/m2' 2247 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. &2240 ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. & 2248 2241 var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. & 2249 2242 var(1:16) == 'usm_t_surf_green' .OR. & … … 2253 2246 ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. & 2254 2247 var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. & 2255 var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis') THEN 2248 var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' .OR. & 2249 var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' ) THEN 2256 2250 unit = '1' 2257 2251 ELSE … … 2302 2296 CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 ) 2303 2297 ENDIF 2298 ! 2299 !-- naheatlayers 2300 IF ( naheatlayers > nzt ) THEN 2301 message_string = 'number of anthropogenic heat layers '// & 2302 '"naheatlayers" can not be larger than'// & 2303 ' number of domain layers "nzt"' 2304 CALL message( 'check_parameters', 'PA0593', 1, 2, 0, 6, 0 ) 2305 ENDIF 2304 2306 2305 2307 END SUBROUTINE usm_check_parameters … … 2331 2333 INTEGER(iwp), PARAMETER :: nd = 5 2332 2334 CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 2333 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /) 2335 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /) 2336 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: diridx = (/ -1, 1, 0, 3, 2 /) 2337 !< index for surf_*_v: 0:3 = (North, South, East, West) 2334 2338 INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart 2335 2339 INTEGER(iwp), DIMENSION(0:nd-1) :: dirend 2336 INTEGER(iwp) :: ids,i surf,isvf,isurfs,isurflt2340 INTEGER(iwp) :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt 2337 2341 INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m 2338 2342 INTEGER(iwp) :: k_topo !< topography top index … … 2351 2355 IF ( var(k-j+1:k) == dirname(i) ) THEN 2352 2356 ids = i 2357 idsint = dirint(ids) 2358 idsidx = diridx(ids) 2353 2359 var = var(:k-j) 2354 2360 EXIT … … 2400 2406 CASE ( 'usm_surfz' ) 2401 2407 !-- array of lw radiation falling to local surface after i-th reflection 2402 DO m = 1, surf_usm_h%ns 2403 i = surf_usm_h%i(m) 2404 j = surf_usm_h%j(m) 2405 k = surf_usm_h%k(m) 2406 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) ) 2407 ENDDO 2408 DO l = 0, 3 2408 IF ( idsint == iup_u ) THEN 2409 DO m = 1, surf_usm_h%ns 2410 i = surf_usm_h%i(m) 2411 j = surf_usm_h%j(m) 2412 k = surf_usm_h%k(m) 2413 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) ) 2414 ENDDO 2415 ELSE 2416 l = idsidx 2409 2417 DO m = 1, surf_usm_v(l)%ns 2410 2418 i = surf_usm_v(l)%i(m) … … 2413 2421 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp ) 2414 2422 ENDDO 2415 END DO2423 ENDIF 2416 2424 2417 2425 CASE ( 'usm_surfcat' ) 2418 2426 !-- surface category 2419 DO m = 1, surf_usm_h%ns 2420 i = surf_usm_h%i(m) 2421 j = surf_usm_h%j(m) 2422 k = surf_usm_h%k(m) 2423 temp_pf(k,j,i) = surf_usm_h%surface_types(m) 2424 ENDDO 2425 DO l = 0, 3 2427 IF ( idsint == iup_u ) THEN 2428 DO m = 1, surf_usm_h%ns 2429 i = surf_usm_h%i(m) 2430 j = surf_usm_h%j(m) 2431 k = surf_usm_h%k(m) 2432 temp_pf(k,j,i) = surf_usm_h%surface_types(m) 2433 ENDDO 2434 ELSE 2435 l = idsidx 2426 2436 DO m = 1, surf_usm_v(l)%ns 2427 2437 i = surf_usm_v(l)%i(m) … … 2430 2440 temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m) 2431 2441 ENDDO 2432 END DO2442 ENDIF 2433 2443 2434 2444 CASE ( 'usm_surfalb' ) 2435 2445 !-- surface albedo, weighted average 2436 DO m = 1, surf_usm_h%ns 2437 i = surf_usm_h%i(m) 2438 j = surf_usm_h%j(m) 2439 k = surf_usm_h%k(m) 2440 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2441 surf_usm_h%albedo(0,m) + & 2442 surf_usm_h%frac(1,m) * & 2443 surf_usm_h%albedo(1,m) + & 2444 surf_usm_h%frac(2,m) * & 2445 surf_usm_h%albedo(2,m) 2446 ENDDO 2447 DO l = 0, 3 2446 IF ( idsint == iup_u ) THEN 2447 DO m = 1, surf_usm_h%ns 2448 i = surf_usm_h%i(m) 2449 j = surf_usm_h%j(m) 2450 k = surf_usm_h%k(m) 2451 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2452 surf_usm_h%albedo(0,m) + & 2453 surf_usm_h%frac(1,m) * & 2454 surf_usm_h%albedo(1,m) + & 2455 surf_usm_h%frac(2,m) * & 2456 surf_usm_h%albedo(2,m) 2457 ENDDO 2458 ELSE 2459 l = idsidx 2448 2460 DO m = 1, surf_usm_v(l)%ns 2449 2461 i = surf_usm_v(l)%i(m) … … 2457 2469 surf_usm_v(l)%albedo(2,m) 2458 2470 ENDDO 2459 END DO2471 ENDIF 2460 2472 2461 2473 CASE ( 'usm_surfemis' ) 2462 2474 !-- surface emissivity, weighted average 2463 DO m = 1, surf_usm_h%ns 2464 i = surf_usm_h%i(m) 2465 j = surf_usm_h%j(m) 2466 k = surf_usm_h%k(m) 2467 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2468 surf_usm_h%emissivity(0,m) + & 2469 surf_usm_h%frac(1,m) * & 2470 surf_usm_h%emissivity(1,m) + & 2471 surf_usm_h%frac(2,m) * & 2472 surf_usm_h%emissivity(2,m) 2473 ENDDO 2474 DO l = 0, 3 2475 IF ( idsint == iup_u ) THEN 2476 DO m = 1, surf_usm_h%ns 2477 i = surf_usm_h%i(m) 2478 j = surf_usm_h%j(m) 2479 k = surf_usm_h%k(m) 2480 temp_pf(k,j,i) = surf_usm_h%frac(0,m) * & 2481 surf_usm_h%emissivity(0,m) + & 2482 surf_usm_h%frac(1,m) * & 2483 surf_usm_h%emissivity(1,m) + & 2484 surf_usm_h%frac(2,m) * & 2485 surf_usm_h%emissivity(2,m) 2486 ENDDO 2487 ELSE 2488 l = idsidx 2475 2489 DO m = 1, surf_usm_v(l)%ns 2476 2490 i = surf_usm_v(l)%i(m) … … 2484 2498 surf_usm_v(l)%emissivity(2,m) 2485 2499 ENDDO 2486 END DO2500 ENDIF 2487 2501 2488 2502 CASE ( 'usm_surfwintrans' ) 2489 2503 !-- transmissivity window tiles 2490 DO m = 1, surf_usm_h%ns 2491 i = surf_usm_h%i(m) 2492 j = surf_usm_h%j(m) 2493 k = surf_usm_h%k(m) 2494 temp_pf(k,j,i) = surf_usm_h%transmissivity(m) 2495 ENDDO 2496 DO l = 0, 3 2504 IF ( idsint == iup_u ) THEN 2505 DO m = 1, surf_usm_h%ns 2506 i = surf_usm_h%i(m) 2507 j = surf_usm_h%j(m) 2508 k = surf_usm_h%k(m) 2509 temp_pf(k,j,i) = surf_usm_h%transmissivity(m) 2510 ENDDO 2511 ELSE 2512 l = idsidx 2497 2513 DO m = 1, surf_usm_v(l)%ns 2498 2514 i = surf_usm_v(l)%i(m) … … 2501 2517 temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m) 2502 2518 ENDDO 2503 2519 ENDIF 2520 2521 CASE ( 'usm_skyvf' ) 2522 !-- sky view factor 2523 DO isurf = dirstart(ids), dirend(ids) 2524 IF ( surfl(id,isurf) == idsint ) THEN 2525 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvf(isurf) 2526 ENDIF 2527 ENDDO 2528 2529 CASE ( 'usm_skyvft' ) 2530 !-- sky view factor 2531 DO isurf = dirstart(ids), dirend(ids) 2532 IF ( surfl(id,isurf) == ids ) THEN 2533 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvft(isurf) 2534 ENDIF 2504 2535 ENDDO 2505 2536 … … 2518 2549 2519 2550 IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. & 2520 surf(iz,isurfs) == ks .AND. surf(id,isurfs) == ids ) THEN2551 surf(iz,isurfs) == ks .AND. surf(id,isurfs) == idsint ) THEN 2521 2552 !-- correct source surface 2522 2553 temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf) … … 2527 2558 !-- array of complete radiation balance 2528 2559 IF ( av == 0 ) THEN 2529 DO m = 1, surf_usm_h%ns 2530 i = surf_usm_h%i(m) 2531 j = surf_usm_h%j(m) 2532 k = surf_usm_h%k(m) 2533 temp_pf(k,j,i) = surf_usm_h%rad_net_l(m) 2534 ENDDO 2535 DO l = 0, 3 2560 IF ( idsint == iup_u ) THEN 2561 DO m = 1, surf_usm_h%ns 2562 i = surf_usm_h%i(m) 2563 j = surf_usm_h%j(m) 2564 k = surf_usm_h%k(m) 2565 temp_pf(k,j,i) = surf_usm_h%rad_net_l(m) 2566 ENDDO 2567 ELSE 2568 l = idsidx 2536 2569 DO m = 1, surf_usm_v(l)%ns 2537 2570 i = surf_usm_v(l)%i(m) … … 2540 2573 temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m) 2541 2574 ENDDO 2542 END DO2575 ENDIF 2543 2576 ELSE 2544 DO m = 1, surf_usm_h%ns 2545 i = surf_usm_h%i(m) 2546 j = surf_usm_h%j(m) 2547 k = surf_usm_h%k(m) 2548 temp_pf(k,j,i) = surf_usm_h%rad_net_av(m) 2549 ENDDO 2550 DO l = 0, 3 2577 IF ( idsint == iup_u ) THEN 2578 DO m = 1, surf_usm_h%ns 2579 i = surf_usm_h%i(m) 2580 j = surf_usm_h%j(m) 2581 k = surf_usm_h%k(m) 2582 temp_pf(k,j,i) = surf_usm_h%rad_net_av(m) 2583 ENDDO 2584 ELSE 2585 l = idsidx 2551 2586 DO m = 1, surf_usm_v(l)%ns 2552 2587 i = surf_usm_v(l)%i(m) … … 2555 2590 temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m) 2556 2591 ENDDO 2557 END DO2592 ENDIF 2558 2593 ENDIF 2559 2594 … … 2561 2596 !-- array of sw radiation falling to surface after i-th reflection 2562 2597 DO isurf = dirstart(ids), dirend(ids) 2563 IF ( surfl(id,isurf) == ids ) THEN2598 IF ( surfl(id,isurf) == idsint ) THEN 2564 2599 IF ( av == 0 ) THEN 2565 2600 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf) … … 2573 2608 !-- array of lw radiation falling to surface after i-th reflection 2574 2609 DO isurf = dirstart(ids), dirend(ids) 2575 IF ( surfl(id,isurf) == ids ) THEN2610 IF ( surfl(id,isurf) == idsint ) THEN 2576 2611 IF ( av == 0 ) THEN 2577 2612 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) … … 2585 2620 !-- array of direct sw radiation falling to surface from sun 2586 2621 DO isurf = dirstart(ids), dirend(ids) 2587 IF ( surfl(id,isurf) == ids ) THEN2622 IF ( surfl(id,isurf) == idsint ) THEN 2588 2623 IF ( av == 0 ) THEN 2589 2624 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf) … … 2597 2632 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 2598 2633 DO isurf = dirstart(ids), dirend(ids) 2599 IF ( surfl(id,isurf) == ids ) THEN2634 IF ( surfl(id,isurf) == idsint ) THEN 2600 2635 IF ( av == 0 ) THEN 2601 2636 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf) … … 2609 2644 !-- array of sw radiation falling to surface from reflections 2610 2645 DO isurf = dirstart(ids), dirend(ids) 2611 IF ( surfl(id,isurf) == ids ) THEN2646 IF ( surfl(id,isurf) == idsint ) THEN 2612 2647 IF ( av == 0 ) THEN 2613 2648 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = & … … 2619 2654 ENDDO 2620 2655 2656 CASE ( 'usm_rad_inlwdif' ) 2657 !-- array of difusion lw radiation falling to surface from sky and borders of the domain 2658 DO isurf = dirstart(ids), dirend(ids) 2659 IF ( surfl(id,isurf) == idsint ) THEN 2660 IF ( av == 0 ) THEN 2661 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf) 2662 ELSE 2663 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf) 2664 ENDIF 2665 ENDIF 2666 ENDDO 2667 2621 2668 CASE ( 'usm_rad_inlwref' ) 2622 2669 !-- array of lw radiation falling to surface from reflections 2623 2670 DO isurf = dirstart(ids), dirend(ids) 2624 IF ( surfl(id,isurf) == ids ) THEN2671 IF ( surfl(id,isurf) == idsint ) THEN 2625 2672 IF ( av == 0 ) THEN 2626 2673 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf) … … 2634 2681 !-- array of sw radiation emitted from surface after i-th reflection 2635 2682 DO isurf = dirstart(ids), dirend(ids) 2636 IF ( surfl(id,isurf) == ids ) THEN2683 IF ( surfl(id,isurf) == idsint ) THEN 2637 2684 IF ( av == 0 ) THEN 2638 2685 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf) … … 2646 2693 !-- array of lw radiation emitted from surface after i-th reflection 2647 2694 DO isurf = dirstart(ids), dirend(ids) 2648 IF ( surfl(id,isurf) == ids ) THEN2695 IF ( surfl(id,isurf) == idsint ) THEN 2649 2696 IF ( av == 0 ) THEN 2650 2697 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf) … … 2658 2705 !-- average of array of residua of sw radiation absorbed in surface after last reflection 2659 2706 DO isurf = dirstart(ids), dirend(ids) 2660 IF ( surfl(id,isurf) == ids ) THEN2707 IF ( surfl(id,isurf) == idsint ) THEN 2661 2708 IF ( av == 0 ) THEN 2662 2709 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf) … … 2670 2717 !-- average of array of residua of lw radiation absorbed in surface after last reflection 2671 2718 DO isurf = dirstart(ids), dirend(ids) 2672 IF ( surfl(id,isurf) == ids ) THEN2719 IF ( surfl(id,isurf) == idsint ) THEN 2673 2720 IF ( av == 0 ) THEN 2674 2721 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf) … … 2682 2729 !-- array of heat flux from radiation for surfaces after all reflections 2683 2730 IF ( av == 0 ) THEN 2684 DO m = 1, surf_usm_h%ns 2685 i = surf_usm_h%i(m) 2686 j = surf_usm_h%j(m) 2687 k = surf_usm_h%k(m) 2688 temp_pf(k,j,i) = surf_usm_h%surfhf(m) 2689 ENDDO 2690 DO l = 0, 3 2731 IF ( idsint == iup_u ) THEN 2732 DO m = 1, surf_usm_h%ns 2733 i = surf_usm_h%i(m) 2734 j = surf_usm_h%j(m) 2735 k = surf_usm_h%k(m) 2736 temp_pf(k,j,i) = surf_usm_h%surfhf(m) 2737 ENDDO 2738 ELSE 2739 l = idsidx 2691 2740 DO m = 1, surf_usm_v(l)%ns 2692 2741 i = surf_usm_v(l)%i(m) … … 2695 2744 temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m) 2696 2745 ENDDO 2697 END DO2746 ENDIF 2698 2747 ELSE 2699 DO m = 1, surf_usm_h%ns 2700 i = surf_usm_h%i(m) 2701 j = surf_usm_h%j(m) 2702 k = surf_usm_h%k(m) 2703 temp_pf(k,j,i) = surf_usm_h%surfhf_av(m) 2704 ENDDO 2705 DO l = 0, 3 2748 IF ( idsint == iup_u ) THEN 2749 DO m = 1, surf_usm_h%ns 2750 i = surf_usm_h%i(m) 2751 j = surf_usm_h%j(m) 2752 k = surf_usm_h%k(m) 2753 temp_pf(k,j,i) = surf_usm_h%surfhf_av(m) 2754 ENDDO 2755 ELSE 2756 l = idsidx 2706 2757 DO m = 1, surf_usm_v(l)%ns 2707 2758 i = surf_usm_v(l)%i(m) … … 2710 2761 temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m) 2711 2762 ENDDO 2712 END DO2763 ENDIF 2713 2764 ENDIF 2714 2765 … … 2716 2767 !-- array of sensible heat flux from surfaces 2717 2768 IF ( av == 0 ) THEN 2718 DO m = 1, surf_usm_h%ns 2719 i = surf_usm_h%i(m) 2720 j = surf_usm_h%j(m) 2721 k = surf_usm_h%k(m) 2722 temp_pf(k,j,i) = surf_usm_h%wshf_eb(m) 2723 ENDDO 2724 DO l = 0, 3 2769 IF ( idsint == iup_u ) THEN 2770 DO m = 1, surf_usm_h%ns 2771 i = surf_usm_h%i(m) 2772 j = surf_usm_h%j(m) 2773 k = surf_usm_h%k(m) 2774 temp_pf(k,j,i) = surf_usm_h%wshf_eb(m) 2775 ENDDO 2776 ELSE 2777 l = idsidx 2725 2778 DO m = 1, surf_usm_v(l)%ns 2726 2779 i = surf_usm_v(l)%i(m) … … 2729 2782 temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m) 2730 2783 ENDDO 2731 END DO2784 ENDIF 2732 2785 ELSE 2733 DO m = 1, surf_usm_h%ns 2734 i = surf_usm_h%i(m) 2735 j = surf_usm_h%j(m) 2736 k = surf_usm_h%k(m) 2737 temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m) 2738 ENDDO 2739 DO l = 0, 3 2786 IF ( idsint == iup_u ) THEN 2787 DO m = 1, surf_usm_h%ns 2788 i = surf_usm_h%i(m) 2789 j = surf_usm_h%j(m) 2790 k = surf_usm_h%k(m) 2791 temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m) 2792 ENDDO 2793 ELSE 2794 l = idsidx 2740 2795 DO m = 1, surf_usm_v(l)%ns 2741 2796 i = surf_usm_v(l)%i(m) … … 2744 2799 temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m) 2745 2800 ENDDO 2746 END DO2801 ENDIF 2747 2802 ENDIF 2748 2803 … … 2751 2806 !-- array of heat flux from ground (land, wall, roof) 2752 2807 IF ( av == 0 ) THEN 2753 DO m = 1, surf_usm_h%ns 2754 i = surf_usm_h%i(m) 2755 j = surf_usm_h%j(m) 2756 k = surf_usm_h%k(m) 2757 temp_pf(k,j,i) = surf_usm_h%wghf_eb(m) 2758 ENDDO 2759 DO l = 0, 3 2808 IF ( idsint == iup_u ) THEN 2809 DO m = 1, surf_usm_h%ns 2810 i = surf_usm_h%i(m) 2811 j = surf_usm_h%j(m) 2812 k = surf_usm_h%k(m) 2813 temp_pf(k,j,i) = surf_usm_h%wghf_eb(m) 2814 ENDDO 2815 ELSE 2816 l = idsidx 2760 2817 DO m = 1, surf_usm_v(l)%ns 2761 2818 i = surf_usm_v(l)%i(m) … … 2764 2821 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m) 2765 2822 ENDDO 2766 END DO2823 ENDIF 2767 2824 ELSE 2768 DO m = 1, surf_usm_h%ns 2769 i = surf_usm_h%i(m) 2770 j = surf_usm_h%j(m) 2771 k = surf_usm_h%k(m) 2772 temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m) 2773 ENDDO 2774 DO l = 0, 3 2825 IF ( idsint == iup_u ) THEN 2826 DO m = 1, surf_usm_h%ns 2827 i = surf_usm_h%i(m) 2828 j = surf_usm_h%j(m) 2829 k = surf_usm_h%k(m) 2830 temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m) 2831 ENDDO 2832 ELSE 2833 l = idsidx 2775 2834 DO m = 1, surf_usm_v(l)%ns 2776 2835 i = surf_usm_v(l)%i(m) … … 2779 2838 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m) 2780 2839 ENDDO 2781 END DO2840 ENDIF 2782 2841 ENDIF 2783 2842 … … 2786 2845 2787 2846 IF ( av == 0 ) THEN 2788 DO m = 1, surf_usm_h%ns 2789 i = surf_usm_h%i(m) 2790 j = surf_usm_h%j(m) 2791 k = surf_usm_h%k(m) 2792 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m) 2793 ENDDO 2794 DO l = 0, 3 2847 IF ( idsint == iup_u ) THEN 2848 DO m = 1, surf_usm_h%ns 2849 i = surf_usm_h%i(m) 2850 j = surf_usm_h%j(m) 2851 k = surf_usm_h%k(m) 2852 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m) 2853 ENDDO 2854 ELSE 2855 l = idsidx 2795 2856 DO m = 1, surf_usm_v(l)%ns 2796 2857 i = surf_usm_v(l)%i(m) … … 2799 2860 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m) 2800 2861 ENDDO 2801 END DO2862 ENDIF 2802 2863 ELSE 2803 DO m = 1, surf_usm_h%ns 2804 i = surf_usm_h%i(m) 2805 j = surf_usm_h%j(m) 2806 k = surf_usm_h%k(m) 2807 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m) 2808 ENDDO 2809 DO l = 0, 3 2864 IF ( idsint == iup_u ) THEN 2865 DO m = 1, surf_usm_h%ns 2866 i = surf_usm_h%i(m) 2867 j = surf_usm_h%j(m) 2868 k = surf_usm_h%k(m) 2869 temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m) 2870 ENDDO 2871 ELSE 2872 l = idsidx 2810 2873 DO m = 1, surf_usm_v(l)%ns 2811 2874 i = surf_usm_v(l)%i(m) … … 2814 2877 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m) 2815 2878 ENDDO 2816 END DO2879 ENDIF 2817 2880 ENDIF 2818 2881 … … 2821 2884 2822 2885 IF ( av == 0 ) THEN 2823 DO m = 1, surf_usm_h%ns 2824 i = surf_usm_h%i(m) 2825 j = surf_usm_h%j(m) 2826 k = surf_usm_h%k(m) 2827 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m) 2828 ENDDO 2829 DO l = 0, 3 2886 IF ( idsint == iup_u ) THEN 2887 DO m = 1, surf_usm_h%ns 2888 i = surf_usm_h%i(m) 2889 j = surf_usm_h%j(m) 2890 k = surf_usm_h%k(m) 2891 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m) 2892 ENDDO 2893 ELSE 2894 l = idsidx 2830 2895 DO m = 1, surf_usm_v(l)%ns 2831 2896 i = surf_usm_v(l)%i(m) … … 2834 2899 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m) 2835 2900 ENDDO 2836 END DO2901 ENDIF 2837 2902 ELSE 2838 DO m = 1, surf_usm_h%ns 2839 i = surf_usm_h%i(m) 2840 j = surf_usm_h%j(m) 2841 k = surf_usm_h%k(m) 2842 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m) 2843 ENDDO 2844 DO l = 0, 3 2903 IF ( idsint == iup_u ) THEN 2904 DO m = 1, surf_usm_h%ns 2905 i = surf_usm_h%i(m) 2906 j = surf_usm_h%j(m) 2907 k = surf_usm_h%k(m) 2908 temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m) 2909 ENDDO 2910 ELSE 2911 l = idsidx 2845 2912 DO m = 1, surf_usm_v(l)%ns 2846 2913 i = surf_usm_v(l)%i(m) … … 2849 2916 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m) 2850 2917 ENDDO 2851 END DO2918 ENDIF 2852 2919 ENDIF 2853 2920 … … 2855 2922 !-- array of heat flux from indoor ground (land, wall, roof) 2856 2923 IF ( av == 0 ) THEN 2857 DO m = 1, surf_usm_h%ns 2858 i = surf_usm_h%i(m) 2859 j = surf_usm_h%j(m) 2860 k = surf_usm_h%k(m) 2861 temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m) 2862 ENDDO 2863 DO l = 0, 3 2924 IF ( idsint == iup_u ) THEN 2925 DO m = 1, surf_usm_h%ns 2926 i = surf_usm_h%i(m) 2927 j = surf_usm_h%j(m) 2928 k = surf_usm_h%k(m) 2929 temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m) 2930 ENDDO 2931 ELSE 2932 l = idsidx 2864 2933 DO m = 1, surf_usm_v(l)%ns 2865 2934 i = surf_usm_v(l)%i(m) … … 2868 2937 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m) 2869 2938 ENDDO 2870 END DO2939 ENDIF 2871 2940 ELSE 2872 DO m = 1, surf_usm_h%ns 2873 i = surf_usm_h%i(m) 2874 j = surf_usm_h%j(m) 2875 k = surf_usm_h%k(m) 2876 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m) 2877 ENDDO 2878 DO l = 0, 3 2941 IF ( idsint == iup_u ) THEN 2942 DO m = 1, surf_usm_h%ns 2943 i = surf_usm_h%i(m) 2944 j = surf_usm_h%j(m) 2945 k = surf_usm_h%k(m) 2946 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m) 2947 ENDDO 2948 ELSE 2949 l = idsidx 2879 2950 DO m = 1, surf_usm_v(l)%ns 2880 2951 i = surf_usm_v(l)%i(m) … … 2883 2954 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m) 2884 2955 ENDDO 2885 END DO2956 ENDIF 2886 2957 ENDIF 2887 2958 … … 2890 2961 2891 2962 IF ( av == 0 ) THEN 2892 DO m = 1, surf_usm_h%ns 2893 i = surf_usm_h%i(m) 2894 j = surf_usm_h%j(m) 2895 k = surf_usm_h%k(m) 2896 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m) 2897 ENDDO 2898 DO l = 0, 3 2963 IF ( idsint == iup_u ) THEN 2964 DO m = 1, surf_usm_h%ns 2965 i = surf_usm_h%i(m) 2966 j = surf_usm_h%j(m) 2967 k = surf_usm_h%k(m) 2968 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m) 2969 ENDDO 2970 ELSE 2971 l = idsidx 2899 2972 DO m = 1, surf_usm_v(l)%ns 2900 2973 i = surf_usm_v(l)%i(m) … … 2903 2976 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m) 2904 2977 ENDDO 2905 END DO2978 ENDIF 2906 2979 ELSE 2907 DO m = 1, surf_usm_h%ns 2908 i = surf_usm_h%i(m) 2909 j = surf_usm_h%j(m) 2910 k = surf_usm_h%k(m) 2911 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m) 2912 ENDDO 2913 DO l = 0, 3 2980 IF ( idsint == iup_u ) THEN 2981 DO m = 1, surf_usm_h%ns 2982 i = surf_usm_h%i(m) 2983 j = surf_usm_h%j(m) 2984 k = surf_usm_h%k(m) 2985 temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m) 2986 ENDDO 2987 ELSE 2988 l = idsidx 2914 2989 DO m = 1, surf_usm_v(l)%ns 2915 2990 i = surf_usm_v(l)%i(m) … … 2918 2993 temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m) 2919 2994 ENDDO 2920 END DO2995 ENDIF 2921 2996 ENDIF 2922 2997 … … 2924 2999 !-- surface temperature for surfaces 2925 3000 IF ( av == 0 ) THEN 2926 DO m = 1, surf_usm_h%ns 2927 i = surf_usm_h%i(m) 2928 j = surf_usm_h%j(m) 2929 k = surf_usm_h%k(m) 2930 temp_pf(k,j,i) = t_surf_h(m) 2931 ENDDO 2932 DO l = 0, 3 3001 IF ( idsint == iup_u ) THEN 3002 DO m = 1, surf_usm_h%ns 3003 i = surf_usm_h%i(m) 3004 j = surf_usm_h%j(m) 3005 k = surf_usm_h%k(m) 3006 temp_pf(k,j,i) = t_surf_h(m) 3007 ENDDO 3008 ELSE 3009 l = idsidx 2933 3010 DO m = 1, surf_usm_v(l)%ns 2934 3011 i = surf_usm_v(l)%i(m) … … 2937 3014 temp_pf(k,j,i) = t_surf_v(l)%t(m) 2938 3015 ENDDO 2939 END DO3016 ENDIF 2940 3017 ELSE 2941 DO m = 1, surf_usm_h%ns 2942 i = surf_usm_h%i(m) 2943 j = surf_usm_h%j(m) 2944 k = surf_usm_h%k(m) 2945 temp_pf(k,j,i) = surf_usm_h%t_surf_av(m) 2946 ENDDO 2947 DO l = 0, 3 3018 IF ( idsint == iup_u ) THEN 3019 DO m = 1, surf_usm_h%ns 3020 i = surf_usm_h%i(m) 3021 j = surf_usm_h%j(m) 3022 k = surf_usm_h%k(m) 3023 temp_pf(k,j,i) = surf_usm_h%t_surf_av(m) 3024 ENDDO 3025 ELSE 3026 l = idsidx 2948 3027 DO m = 1, surf_usm_v(l)%ns 2949 3028 i = surf_usm_v(l)%i(m) … … 2952 3031 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m) 2953 3032 ENDDO 2954 END DO3033 ENDIF 2955 3034 ENDIF 2956 3035 … … 2959 3038 2960 3039 IF ( av == 0 ) THEN 2961 DO m = 1, surf_usm_h%ns 2962 i = surf_usm_h%i(m) 2963 j = surf_usm_h%j(m) 2964 k = surf_usm_h%k(m) 2965 temp_pf(k,j,i) = t_surf_window_h(m) 2966 ENDDO 2967 DO l = 0, 3 3040 IF ( idsint == iup_u ) THEN 3041 DO m = 1, surf_usm_h%ns 3042 i = surf_usm_h%i(m) 3043 j = surf_usm_h%j(m) 3044 k = surf_usm_h%k(m) 3045 temp_pf(k,j,i) = t_surf_window_h(m) 3046 ENDDO 3047 ELSE 3048 l = idsidx 2968 3049 DO m = 1, surf_usm_v(l)%ns 2969 3050 i = surf_usm_v(l)%i(m) … … 2972 3053 temp_pf(k,j,i) = t_surf_window_v(l)%t(m) 2973 3054 ENDDO 2974 END DO3055 ENDIF 2975 3056 2976 3057 ELSE 2977 DO m = 1, surf_usm_h%ns 2978 i = surf_usm_h%i(m) 2979 j = surf_usm_h%j(m) 2980 k = surf_usm_h%k(m) 2981 temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m) 2982 ENDDO 2983 DO l = 0, 3 3058 IF ( idsint == iup_u ) THEN 3059 DO m = 1, surf_usm_h%ns 3060 i = surf_usm_h%i(m) 3061 j = surf_usm_h%j(m) 3062 k = surf_usm_h%k(m) 3063 temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m) 3064 ENDDO 3065 ELSE 3066 l = idsidx 2984 3067 DO m = 1, surf_usm_v(l)%ns 2985 3068 i = surf_usm_v(l)%i(m) … … 2989 3072 ENDDO 2990 3073 2991 END DO3074 ENDIF 2992 3075 2993 3076 ENDIF … … 2997 3080 2998 3081 IF ( av == 0 ) THEN 2999 DO m = 1, surf_usm_h%ns 3000 i = surf_usm_h%i(m) 3001 j = surf_usm_h%j(m) 3002 k = surf_usm_h%k(m) 3003 temp_pf(k,j,i) = t_surf_green_h(m) 3004 ENDDO 3005 DO l = 0, 3 3082 IF ( idsint == iup_u ) THEN 3083 DO m = 1, surf_usm_h%ns 3084 i = surf_usm_h%i(m) 3085 j = surf_usm_h%j(m) 3086 k = surf_usm_h%k(m) 3087 temp_pf(k,j,i) = t_surf_green_h(m) 3088 ENDDO 3089 ELSE 3090 l = idsidx 3006 3091 DO m = 1, surf_usm_v(l)%ns 3007 3092 i = surf_usm_v(l)%i(m) … … 3010 3095 temp_pf(k,j,i) = t_surf_green_v(l)%t(m) 3011 3096 ENDDO 3012 END DO3097 ENDIF 3013 3098 3014 3099 ELSE 3015 DO m = 1, surf_usm_h%ns 3016 i = surf_usm_h%i(m) 3017 j = surf_usm_h%j(m) 3018 k = surf_usm_h%k(m) 3019 temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m) 3020 ENDDO 3021 DO l = 0, 3 3100 IF ( idsint == iup_u ) THEN 3101 DO m = 1, surf_usm_h%ns 3102 i = surf_usm_h%i(m) 3103 j = surf_usm_h%j(m) 3104 k = surf_usm_h%k(m) 3105 temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m) 3106 ENDDO 3107 ELSE 3108 l = idsidx 3022 3109 DO m = 1, surf_usm_v(l)%ns 3023 3110 i = surf_usm_v(l)%i(m) … … 3027 3114 ENDDO 3028 3115 3029 END DO3116 ENDIF 3030 3117 3031 3118 ENDIF … … 3035 3122 3036 3123 IF ( av == 0 ) THEN 3037 DO m = 1, surf_usm_h%ns 3038 i = surf_usm_h%i(m) 3039 j = surf_usm_h%j(m) 3040 k = surf_usm_h%k(m) 3041 temp_pf(k,j,i) = t_surf_10cm_h(m) 3042 ENDDO 3043 DO l = 0, 3 3124 IF ( idsint == iup_u ) THEN 3125 DO m = 1, surf_usm_h%ns 3126 i = surf_usm_h%i(m) 3127 j = surf_usm_h%j(m) 3128 k = surf_usm_h%k(m) 3129 temp_pf(k,j,i) = t_surf_10cm_h(m) 3130 ENDDO 3131 ELSE 3132 l = idsidx 3044 3133 DO m = 1, surf_usm_v(l)%ns 3045 3134 i = surf_usm_v(l)%i(m) … … 3048 3137 temp_pf(k,j,i) = t_surf_10cm_v(l)%t(m) 3049 3138 ENDDO 3050 END DO3139 ENDIF 3051 3140 3052 3141 ELSE 3053 DO m = 1, surf_usm_h%ns 3054 i = surf_usm_h%i(m) 3055 j = surf_usm_h%j(m) 3056 k = surf_usm_h%k(m) 3057 temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m) 3058 ENDDO 3059 DO l = 0, 3 3142 IF ( idsint == iup_u ) THEN 3143 DO m = 1, surf_usm_h%ns 3144 i = surf_usm_h%i(m) 3145 j = surf_usm_h%j(m) 3146 k = surf_usm_h%k(m) 3147 temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m) 3148 ENDDO 3149 ELSE 3150 l = idsidx 3060 3151 DO m = 1, surf_usm_v(l)%ns 3061 3152 i = surf_usm_v(l)%i(m) … … 3065 3156 ENDDO 3066 3157 3067 END DO3158 ENDIF 3068 3159 3069 3160 ENDIF … … 3073 3164 !-- wall temperature for iwl layer of walls and land 3074 3165 IF ( av == 0 ) THEN 3075 DO m = 1, surf_usm_h%ns 3076 i = surf_usm_h%i(m) 3077 j = surf_usm_h%j(m) 3078 k = surf_usm_h%k(m) 3079 temp_pf(k,j,i) = t_wall_h(iwl,m) 3080 ENDDO 3081 DO l = 0, 3 3166 IF ( idsint == iup_u ) THEN 3167 DO m = 1, surf_usm_h%ns 3168 i = surf_usm_h%i(m) 3169 j = surf_usm_h%j(m) 3170 k = surf_usm_h%k(m) 3171 temp_pf(k,j,i) = t_wall_h(iwl,m) 3172 ENDDO 3173 ELSE 3174 l = idsidx 3082 3175 DO m = 1, surf_usm_v(l)%ns 3083 3176 i = surf_usm_v(l)%i(m) … … 3086 3179 temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m) 3087 3180 ENDDO 3088 END DO3181 ENDIF 3089 3182 ELSE 3090 DO m = 1, surf_usm_h%ns 3091 i = surf_usm_h%i(m) 3092 j = surf_usm_h%j(m) 3093 k = surf_usm_h%k(m) 3094 temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m) 3095 ENDDO 3096 DO l = 0, 3 3183 IF ( idsint == iup_u ) THEN 3184 DO m = 1, surf_usm_h%ns 3185 i = surf_usm_h%i(m) 3186 j = surf_usm_h%j(m) 3187 k = surf_usm_h%k(m) 3188 temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m) 3189 ENDDO 3190 ELSE 3191 l = idsidx 3097 3192 DO m = 1, surf_usm_v(l)%ns 3098 3193 i = surf_usm_v(l)%i(m) … … 3101 3196 temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m) 3102 3197 ENDDO 3103 END DO3198 ENDIF 3104 3199 ENDIF 3105 3200 … … 3107 3202 !-- window temperature for iwl layer of walls and land 3108 3203 IF ( av == 0 ) THEN 3109 DO m = 1, surf_usm_h%ns 3110 i = surf_usm_h%i(m) 3111 j = surf_usm_h%j(m) 3112 k = surf_usm_h%k(m) 3113 temp_pf(k,j,i) = t_window_h(iwl,m) 3114 ENDDO 3115 DO l = 0, 3 3204 IF ( idsint == iup_u ) THEN 3205 DO m = 1, surf_usm_h%ns 3206 i = surf_usm_h%i(m) 3207 j = surf_usm_h%j(m) 3208 k = surf_usm_h%k(m) 3209 temp_pf(k,j,i) = t_window_h(iwl,m) 3210 ENDDO 3211 ELSE 3212 l = idsidx 3116 3213 DO m = 1, surf_usm_v(l)%ns 3117 3214 i = surf_usm_v(l)%i(m) … … 3120 3217 temp_pf(k,j,i) = t_window_v(l)%t(iwl,m) 3121 3218 ENDDO 3122 END DO3219 ENDIF 3123 3220 ELSE 3124 DO m = 1, surf_usm_h%ns 3125 i = surf_usm_h%i(m) 3126 j = surf_usm_h%j(m) 3127 k = surf_usm_h%k(m) 3128 temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m) 3129 ENDDO 3130 DO l = 0, 3 3221 IF ( idsint == iup_u ) THEN 3222 DO m = 1, surf_usm_h%ns 3223 i = surf_usm_h%i(m) 3224 j = surf_usm_h%j(m) 3225 k = surf_usm_h%k(m) 3226 temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m) 3227 ENDDO 3228 ELSE 3229 l = idsidx 3131 3230 DO m = 1, surf_usm_v(l)%ns 3132 3231 i = surf_usm_v(l)%i(m) … … 3135 3234 temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m) 3136 3235 ENDDO 3137 END DO3236 ENDIF 3138 3237 ENDIF 3139 3238 … … 3141 3240 !-- green temperature for iwl layer of walls and land 3142 3241 IF ( av == 0 ) THEN 3143 DO m = 1, surf_usm_h%ns 3144 i = surf_usm_h%i(m) 3145 j = surf_usm_h%j(m) 3146 k = surf_usm_h%k(m) 3147 temp_pf(k,j,i) = t_green_h(iwl,m) 3148 ENDDO 3149 DO l = 0, 3 3242 IF ( idsint == iup_u ) THEN 3243 DO m = 1, surf_usm_h%ns 3244 i = surf_usm_h%i(m) 3245 j = surf_usm_h%j(m) 3246 k = surf_usm_h%k(m) 3247 temp_pf(k,j,i) = t_green_h(iwl,m) 3248 ENDDO 3249 ELSE 3250 l = idsidx 3150 3251 DO m = 1, surf_usm_v(l)%ns 3151 3252 i = surf_usm_v(l)%i(m) … … 3154 3255 temp_pf(k,j,i) = t_green_v(l)%t(iwl,m) 3155 3256 ENDDO 3156 END DO3257 ENDIF 3157 3258 ELSE 3158 DO m = 1, surf_usm_h%ns 3159 i = surf_usm_h%i(m) 3160 j = surf_usm_h%j(m) 3161 k = surf_usm_h%k(m) 3162 temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m) 3163 ENDDO 3164 DO l = 0, 3 3259 IF ( idsint == iup_u ) THEN 3260 DO m = 1, surf_usm_h%ns 3261 i = surf_usm_h%i(m) 3262 j = surf_usm_h%j(m) 3263 k = surf_usm_h%k(m) 3264 temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m) 3265 ENDDO 3266 ELSE 3267 l = idsidx 3165 3268 DO m = 1, surf_usm_v(l)%ns 3166 3269 i = surf_usm_v(l)%i(m) … … 3169 3272 temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m) 3170 3273 ENDDO 3171 END DO3274 ENDIF 3172 3275 ENDIF 3173 3276 … … 3180 3283 ! 3181 3284 !-- Rearrange dimensions for NetCDF output 3285 !-- FIXME: this may generate FPE overflow upon conversion from DP to SP 3182 3286 DO j = nys, nyn 3183 3287 DO i = nxl, nxr 3184 3288 DO k = nzb_do, nzt_do 3289 ! print*, j,nys,nyn,i,nxl,nxr,k,nzb_do,nzt_do,local_pf(i,j,k),temp_pf(k,j,i) 3185 3290 local_pf(i,j,k) = temp_pf(k,j,i) 3186 3291 ENDDO … … 3228 3333 var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. & 3229 3334 var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' .OR. & 3230 var(1:16) == 'usm_surfwintrans' ) THEN 3335 var(1:16) == 'usm_surfwintrans' .OR. & 3336 var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' ) THEN 3231 3337 3232 3338 found = .TRUE. … … 3416 3522 INTEGER(iwp) :: st !< dummy 3417 3523 3418 REAL(wp) :: c, d, tin, twin , exn3524 REAL(wp) :: c, d, tin, twin 3419 3525 REAL(wp) :: ground_floor_level_l !< local height of ground floor level 3420 3526 REAL(wp) :: z_agl !< height above ground 3527 REAL(wp), DIMENSION(nzb:nzt) :: exn !< value of the Exner function in layers 3421 3528 3422 3529 ! … … 4402 4509 CALL usm_read_anthropogenic_heat() 4403 4510 ENDIF 4404 4511 4405 4512 IF ( plant_canopy ) THEN 4406 4513 … … 4419 4526 4420 4527 !-- Calculate initial surface temperature from pt of adjacent gridbox 4421 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 4528 #if ! defined( __nopointer ) 4529 exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp !< Exner function 4530 #endif 4422 4531 4423 4532 ! … … 4430 4539 k = surf_usm_h%k(m) 4431 4540 4432 t_surf_h(m) = pt(k,j,i) * exn 4433 t_surf_window_h(m) = pt(k,j,i) * exn 4434 t_surf_green_h(m) = pt(k,j,i) * exn 4435 surf_usm_h%pt_surface(m) = pt(k,j,i) * exn 4541 t_surf_h(m) = pt(k,j,i) * exn(k) 4542 t_surf_window_h(m) = pt(k,j,i) * exn(k) 4543 t_surf_green_h(m) = pt(k,j,i) * exn(k) 4544 surf_usm_h%pt_surface(m) = pt(k,j,i) * exn(k) 4436 4545 ENDDO 4437 4546 ! … … 4443 4552 k = surf_usm_v(l)%k(m) 4444 4553 4445 t_surf_v(l)%t(m) = pt(k,j,i) * exn 4446 t_surf_window_v(l)%t(m) = pt(k,j,i) * exn 4447 t_surf_green_v(l)%t(m) = pt(k,j,i) * exn 4448 surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exn 4554 t_surf_v(l)%t(m) = pt(k,j,i) * exn(k) 4555 t_surf_window_v(l)%t(m) = pt(k,j,i) * exn(k) 4556 t_surf_green_v(l)%t(m) = pt(k,j,i) * exn(k) 4557 surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exn(k) 4449 4558 ENDDO 4450 4559 ENDDO … … 4496 4605 ENDDO 4497 4606 ENDDO 4498 4607 ELSE 4608 !-- If specified, replace constant wall temperatures with fully 3D values from file 4609 IF ( read_wall_temp_3d ) CALL usm_read_wall_temperature() 4610 ! 4499 4611 ENDIF 4500 4612 … … 4519 4631 t_green_h_p = t_green_h 4520 4632 t_green_v_p = t_green_v 4633 4634 !-- Adjust radiative fluxes for urban surface at model start 4635 !CALL radiation_interaction 4636 !-- TODO: interaction should be called once before first output, 4637 !-- that is not yet possible. 4521 4638 4522 4639 CALL cpu_log( log_point_s(78), 'usm_init', 'stop' ) … … 5029 5146 building_type, & 5030 5147 land_category, & 5031 pedestrant_category, & 5032 ra_horiz_coef, & 5148 naheatlayers, & 5149 pedestrian_category, & 5150 roughness_concrete, & 5151 read_wall_temp_3d, & 5033 5152 roof_category, & 5034 5153 urban_surface, & 5035 5154 usm_anthropogenic_heat, & 5036 5155 usm_material_model, & 5037 usm_lad_rma, &5038 5156 wall_category, & 5039 indoor_model 5040 5041 line = ' ' 5157 indoor_model, & 5158 wall_inner_temperature, & 5159 roof_inner_temperature, & 5160 soil_inner_temperature, & 5161 window_inner_temperature 5042 5162 5043 5163 ! … … 5134 5254 SUBROUTINE usm_read_anthropogenic_heat 5135 5255 5136 INTEGER(iwp) :: i,j, ii5256 INTEGER(iwp) :: i,j,k,ii 5137 5257 REAL(wp) :: heat 5138 5258 5139 5259 !-- allocation of array of sources of anthropogenic heat and their diural profile 5140 ALLOCATE( aheat(n ys:nyn,nxl:nxr) )5141 ALLOCATE( aheatprof( 0:24) )5260 ALLOCATE( aheat(naheatlayers,nys:nyn,nxl:nxr) ) 5261 ALLOCATE( aheatprof(naheatlayers,0:24) ) 5142 5262 5143 5263 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 5154 5274 j = 0 5155 5275 DO 5156 READ( 151, *, err=12, end=13 ) i, j, heat5276 READ( 151, *, err=12, end=13 ) i, j, k, heat 5157 5277 IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN 5158 !-- write heat into the array 5159 aheat(j,i) = heat 5278 IF ( k <= naheatlayers .AND. k > get_topography_top_index_ji( j, i, 's' ) ) THEN 5279 !-- write heat into the array 5280 aheat(k,j,i) = heat 5281 ENDIF 5160 5282 ENDIF 5161 5283 CYCLE … … 5186 5308 i = 0 5187 5309 DO 5188 READ( 151, *, err=22, end=23 ) i, heat5189 IF ( i >= 0 .AND. i <= 24 ) THEN5310 READ( 151, *, err=22, end=23 ) i, k, heat 5311 IF ( i >= 0 .AND. i <= 24 .AND. k <= naheatlayers ) THEN 5190 5312 !-- write heat into the array 5191 aheatprof( i) = heat5313 aheatprof(k,i) = heat 5192 5314 ENDIF 5193 5315 CYCLE … … 5196 5318 CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 ) 5197 5319 ENDDO 5198 aheatprof( 24) = aheatprof(0)5320 aheatprof(:,24) = aheatprof(:,0) 5199 5321 23 CLOSE(151) 5200 5322 CYCLE … … 5212 5334 5213 5335 !------------------------------------------------------------------------------! 5214 !5215 5336 ! Description: 5216 5337 ! ------------ 5217 !> Soubroutine reads t_surf and t_wall data from restart file(s) 5218 !kanani: Renamed this routine according to corresponging routines in PALM 5219 !kanani: Modified the routine to match rrd_global, from where usm_rrd_local 5220 ! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod) 5221 ! Also, I had some trouble with the allocation of t_surf, since this is a pointer. 5222 ! So, I added some directives here. 5338 !> Soubroutine reads t_surf and t_wall data from restart files 5223 5339 !------------------------------------------------------------------------------! 5224 5340 SUBROUTINE usm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & … … 6623 6739 surf_usm_h%albedo(:,m) = surface_params(ialbedo,ip) 6624 6740 ENDIF 6741 !-- Albedo type is 0 (custom), others are replaced later 6742 surf_usm_h%albedo_type(:,m) = 0 6625 6743 !-- Transmissivity 6626 6744 IF ( surf_usm_h%transmissivity(m) < 0.0_wp ) THEN … … 6636 6754 surf_usm_h%lambda_surf_green(m) = surface_params(ilambdas,ip) 6637 6755 ! 6638 !-- roughness relative to concrete6756 !-- roughness length for momentum, heat and humidity 6639 6757 surf_usm_h%z0(m) = surface_params(irough,ip) 6640 ! 6758 surf_usm_h%z0h(m) = surface_params(iroughh,ip) 6759 surf_usm_h%z0q(m) = surface_params(iroughh,ip) 6760 ! 6641 6761 !-- Surface skin layer heat capacity (J mâ2 Kâ1 ) 6642 6762 surf_usm_h%c_surface(m) = surface_params(icsurf,ip) … … 6676 6796 j = surf_usm_v(l)%j(m) 6677 6797 kw = surf_usm_v(l)%k(m) 6678 6798 6679 6799 IF ( l == 3 ) THEN ! westward facing 6680 6800 iw = i … … 6699 6819 ENDIF 6700 6820 6701 IF ( kw <= usm_par(ii,jw,iw) ) THEN 6702 !-- pedestrant zone 6821 IF ( iw < 0 .OR. jw < 0 ) THEN 6822 !-- wall on west or south border of the domain - assign default category 6823 IF ( kw <= roof_height_limit ) THEN 6824 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface in wall zone 6825 ELSE 6826 surf_usm_v(l)%surface_types(m) = roof_category !< default category for wall surface in roof zone 6827 END IF 6828 surf_usm_v(l)%albedo(:,m) = -1.0_wp 6829 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 6830 ELSE IF ( kw <= usm_par(ii,jw,iw) ) THEN 6831 !-- pedestrian zone 6703 6832 IF ( usm_par(ii+1,jw,iw) == 0 ) THEN 6704 surf_usm_v(l)%surface_types(m) = pedestr ant_category !< default category for wall surface in pedestrantzone6833 surf_usm_v(l)%surface_types(m) = pedestrian_category !< default category for wall surface in pedestrian zone 6705 6834 surf_usm_v(l)%albedo(:,m) = -1.0_wp 6706 6835 surf_usm_v(l)%thickness_wall(m) = -1.0_wp … … 6751 6880 ENDIF 6752 6881 ELSE 6753 !-- something wrong 6754 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 ) 6882 ! 6883 !-- supply the default category 6884 IF ( kw <= roof_height_limit ) THEN 6885 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface in wall zone 6886 ELSE 6887 surf_usm_v(l)%surface_types(m) = roof_category !< default category for wall surface in roof zone 6888 END IF 6889 surf_usm_v(l)%albedo(:,m) = -1.0_wp 6890 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 6755 6891 ENDIF 6756 6757 6892 ! 6758 6893 !-- Find the type position … … 6767 6902 IF ( ip == -99999 ) THEN 6768 6903 !-- wall category not found 6769 WRITE (message_string, "(A,I 5,A,3I5)") 'wall category ', it, &6904 WRITE (message_string, "(A,I7,A,3I5)") 'wall category ', it, & 6770 6905 ' not found for i,j,k=', iw,jw,kw 6771 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )6906 WRITE(9,*) message_string 6772 6907 ENDIF 6773 6908 ! … … 6776 6911 surf_usm_v(l)%albedo(:,m) = surface_params(ialbedo,ip) 6777 6912 ENDIF 6913 !-- Albedo type is 0 (custom), others are replaced later 6914 surf_usm_v(l)%albedo_type(:,m) = 0 6778 6915 !-- Transmissivity of the windows 6779 6916 IF ( surf_usm_v(l)%transmissivity(m) < 0.0_wp ) THEN … … 6784 6921 surf_usm_v(l)%emissivity(:,m) = surface_params(iemiss,ip) 6785 6922 ! 6786 !-- heat conductivity λS between air and wall ( W mâ2 Kâ1 )6923 !-- heat conductivity lambda S between air and wall ( W m-2 K-1 ) 6787 6924 surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip) 6788 6925 surf_usm_v(l)%lambda_surf_window(m) = surface_params(ilambdas,ip) 6789 6926 surf_usm_v(l)%lambda_surf_green(m) = surface_params(ilambdas,ip) 6790 6927 ! 6791 !-- roughness relative to concrete6928 !-- roughness length 6792 6929 surf_usm_v(l)%z0(m) = surface_params(irough,ip) 6930 surf_usm_v(l)%z0h(m) = surface_params(iroughh,ip) 6931 surf_usm_v(l)%z0q(m) = surface_params(iroughh,ip) 6793 6932 ! 6794 !-- Surface skin layer heat capacity (J m â2 Kâ1 )6933 !-- Surface skin layer heat capacity (J m-2 K-1 ) 6795 6934 surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip) 6796 6935 surf_usm_v(l)%c_surface_window(m) = surface_params(icsurf,ip) … … 6809 6948 surf_usm_v(l)%thickness_green(m) = surface_params(ithick,ip) 6810 6949 ENDIF 6811 ! 6812 !-- volumetric heat capacity rho*C of the wall ( J m â3 Kâ1 )6950 ! 6951 !-- volumetric heat capacity rho*C of the wall ( J m-3 K-1 ) 6813 6952 surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip) 6814 6953 surf_usm_v(l)%rho_c_window(:,m) = surface_params(irhoC,ip) 6815 6954 surf_usm_v(l)%rho_c_green(:,m) = surface_params(irhoC,ip) 6816 6955 ! 6817 !-- thermal conductivity λH of the wall (W mâ1 Kâ1 )6956 !-- thermal conductivity lambda H of the wall (W m-1 K-1 ) 6818 6957 surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip) 6819 6958 surf_usm_v(l)%lambda_h_window(:,m) = surface_params(ilambdah,ip) … … 6848 6987 6849 6988 END SUBROUTINE usm_read_urban_surface_types 6989 6990 6991 !------------------------------------------------------------------------------! 6992 ! Description: 6993 ! ------------ 6994 ! 6995 !> This function advances through the list of local surfaces to find given 6996 !> x, y, d, z coordinates 6997 !------------------------------------------------------------------------------! 6998 PURE FUNCTION advance_surface(isurfl_start, isurfl_stop, x, y, z, d) & 6999 result(isurfl) 7000 7001 INTEGER(iwp), INTENT(in) :: isurfl_start, isurfl_stop 7002 INTEGER(iwp), INTENT(in) :: x, y, z, d 7003 INTEGER(iwp) :: isx, isy, isz, isd 7004 INTEGER(iwp) :: isurfl 7005 7006 DO isurfl = isurfl_start, isurfl_stop 7007 isx = surfl(ix, isurfl) 7008 isy = surfl(iy, isurfl) 7009 isz = surfl(iz, isurfl) 7010 isd = surfl(id, isurfl) 7011 IF ( isx==x .and. isy==y .and. isz==z .and. isd==d ) RETURN 7012 ENDDO 7013 7014 !-- coordinate not found 7015 isurfl = -1 7016 7017 END FUNCTION 7018 7019 7020 !------------------------------------------------------------------------------! 7021 ! Description: 7022 ! ------------ 7023 ! 7024 !> This subroutine reads temperatures of respective material layers in walls, 7025 !> roofs and ground from input files. Data in the input file must be in 7026 !> standard order, i.e. horizontal surfaces first ordered by x, y and then 7027 !> vertical surfaces ordered by x, y, direction, z 7028 !------------------------------------------------------------------------------! 7029 SUBROUTINE usm_read_wall_temperature 7030 7031 INTEGER(iwp) :: i, j, k, d, ii, iline 7032 INTEGER(iwp) :: isurfl 7033 REAL(wp) :: rtsurf 7034 REAL(wp), DIMENSION(nzb_wall:nzt_wall+1) :: rtwall 7035 7036 7037 7038 7039 DO ii = 0, io_blocks-1 7040 IF ( ii == io_group ) THEN 7041 7042 !-- open wall temperature file 7043 OPEN( 152, file='WALL_TEMPERATURE'//coupling_char, action='read', & 7044 status='old', form='formatted', err=15 ) 7045 7046 isurfl = 0 7047 iline = 1 7048 DO 7049 rtwall = -9999.0_wp !< for incomplete lines 7050 READ( 152, *, err=13, end=14 ) i, j, k, d, rtsurf, rtwall 7051 7052 IF ( nxl <= i .and. i <= nxr .and. & 7053 nys <= j .and. j <= nyn) THEN !< local processor 7054 !-- identify surface id 7055 isurfl = advance_surface(isurfl+1, nsurfl, i, j, k, d) 7056 IF ( isurfl == -1 ) THEN 7057 WRITE(message_string, '(a,4i5,a,i5,a)') 'Coordinates (xyzd) ', i, j, k, d, & 7058 ' on line ', iline, & 7059 ' in file WALL_TEMPERATURE are either not present or out of standard order of surfaces.' 7060 CALL message( 'usm_read_wall_temperature', 'PA0521', 1, 2, 0, 6, 0 ) 7061 ENDIF 7062 7063 !-- assign temperatures 7064 IF ( d == 0 ) THEN 7065 t_surf_h(isurfl) = rtsurf 7066 t_wall_h(:,isurfl) = rtwall(:) 7067 ELSE 7068 t_surf_v(d-1)%t(isurfl) = rtsurf 7069 t_wall_v(d-1)%t(:,isurfl) = rtwall(:) 7070 ENDIF 7071 ENDIF 7072 7073 iline = iline + 1 7074 CYCLE 7075 13 WRITE(message_string, '(a,i5,a)') 'Error reading line ', iline, & 7076 ' in file WALL_TEMPERATURE.' 7077 CALL message( 'usm_read_wall_temperature', 'PA0522', 1, 2, 0, 6, 0 ) 7078 ENDDO 7079 14 CLOSE(152) 7080 CYCLE 7081 15 message_string = 'file WALL_TEMPERATURE'//TRIM(coupling_char)//' does not exist' 7082 CALL message( 'usm_read_wall_temperature', 'PA0523', 1, 2, 0, 6, 0 ) 7083 ENDIF 7084 #if defined( __parallel ) && ! defined ( __check ) 7085 CALL MPI_BARRIER( comm2d, ierr ) 7086 #endif 7087 ENDDO 7088 7089 CALL location_message( ' wall layer temperatures read', .TRUE. ) 7090 7091 END SUBROUTINE usm_read_wall_temperature 7092 6850 7093 6851 7094 … … 6945 7188 IF ( surf_usm_h%r_a_window(m) < 1.0_wp ) & 6946 7189 surf_usm_h%r_a_window(m) = 1.0_wp 6947 6948 6949 !-- the parameterization is developed originally for larger scales6950 !-- (compare with remark in TUF-3D)6951 !-- our first experiences show that the parameterization underestimates6952 !-- r_a in meter resolution.6953 !-- A temporary solution would be multiplication by magic constant :-(.6954 !-- For the moment this is comment out.6955 surf_usm_h%r_a(m) = surf_usm_h%r_a(m) !* ra_horiz_coef6956 surf_usm_h%r_a_window(m) = surf_usm_h%r_a_window(m) !* ra_horiz_coef6957 surf_usm_h%r_a_green(m) = surf_usm_h%r_a_green(m) !* ra_horiz_coef6958 7190 6959 7191 !-- factor for shf_eb … … 7070 7302 !-- calculate fluxes 7071 7303 !-- rad_net_l is never used! 7072 surf_usm_h%rad_net_l(m) = surf_usm_h%frac(0,m) * & 7073 ( surf_usm_h%rad_net_l(m) + & 7074 3.0_wp * sigma_sb * & 7075 t_surf_h(m)**4 - 4.0_wp * sigma_sb * & 7076 t_surf_h(m)**3 * t_surf_h_p(m) ) & 7077 + surf_usm_h%frac(2,m) * & 7078 ( surf_usm_h%rad_net_l(m) + & 7079 3.0_wp * sigma_sb * & 7080 t_surf_window_h(m)**4 - 4.0_wp * sigma_sb * & 7081 t_surf_window_h(m)**3 * t_surf_window_h_p(m) ) & 7082 + surf_usm_h%frac(1,m) * & 7083 ( surf_usm_h%rad_net_l(m) + & 7084 3.0_wp * sigma_sb * & 7085 t_surf_green_h(m)**4 - 4.0_wp * sigma_sb * & 7086 t_surf_green_h(m)**3 * t_surf_green_h_p(m) ) 7304 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) + & 7305 surf_usm_h%frac(0,m) * & 7306 sigma_sb * surf_usm_h%emissivity(0,m) * & 7307 ( t_surf_h_p(m)**4 - t_surf_h(m)**4 ) & 7308 + surf_usm_h%frac(2,m) * & 7309 sigma_sb * surf_usm_h%emissivity(2,m) * & 7310 ( t_surf_window_h_p(m)**4 - t_surf_window_h(m)**4 ) & 7311 + surf_usm_h%frac(1,m) * & 7312 sigma_sb * surf_usm_h%emissivity(1,m) * & 7313 ( t_surf_green_h_p(m)**4 - t_surf_green_h(m)**4 ) 7314 7087 7315 surf_usm_h%wghf_eb(m) = lambda_surface * & 7088 7316 ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) ) … … 7269 7497 7270 7498 !-- calculate fluxes 7271 !-- rad_net_l is never used!7499 !-- prognostic rad_net_l is used just for output! 7272 7500 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%frac(0,m) * & 7273 7501 ( surf_usm_v(l)%rad_net_l(m) + & … … 7322 7550 dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp) 7323 7551 dhour = INT(dtime/3600.0_wp) 7324 !-- linear interpolation of coeficient 7325 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1) 7326 7327 DO m = 1, surf_usm_h%ns 7328 ! 7552 DO m = 1, naheatlayers 7329 7553 !-- Get indices of respective grid point 7330 7554 i = surf_usm_h%i(m) 7331 7555 j = surf_usm_h%j(m) 7332 7556 k = surf_usm_h%k(m) 7333 7334 IF ( aheat(j,i) > 0.0_wp) THEN7335 !-- TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1in time dt_3d7557 IF ( k > get_topography_top_index_ji( j, i, 's' ) .AND. & 7558 k <= naheatlayers ) THEN 7559 !-- increase of pt in box i,j,k in time dt_3d 7336 7560 !-- given to anthropogenic heat aheat*acoef (W*m-2) 7337 !-- k = nzb_s_inner(j,i)+1 7338 !-- pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz) 7339 !-- Instead of this, we can adjust shf in case AH only at surface 7340 surf_usm_h%shf(m) = surf_usm_h%shf(m) + & 7341 aheat(j,i) * acoef * ddx * ddy / cp 7561 !-- linear interpolation of coeficient 7562 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + & 7563 (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1) 7564 IF ( aheat(k,j,i) > 0.0_wp ) THEN 7565 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz) 7566 ENDIF 7342 7567 ENDIF 7343 7568 ENDDO … … 7437 7662 7438 7663 !------------------------------------------------------------------------------! 7439 !7440 7664 ! Description: 7441 7665 ! ------------ 7442 !> Subroutine writes t_surf and t_wall data into restart files. It is necessary 7443 !> to write start_index and end_index several times 7444 !kanani: Renamed this routine according to corresponging routines in PALM 7445 !kanani: Modified the routine to match wrd_global, from where usm_wrd_local 7446 ! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod) 7447 ! Also, I had some trouble with the allocation of t_surf, since this is a pointer. 7448 ! So, I added some directives here. 7666 !> Subroutine writes t_surf and t_wall data into restart files 7449 7667 !------------------------------------------------------------------------------! 7450 7668 SUBROUTINE usm_wrd_local
Note: See TracChangeset
for help on using the changeset viewer.