Changeset 2011 for palm/trunk/SOURCE/urban_surface_mod.f90
 Timestamp:
 Sep 19, 2016 5:29:57 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/urban_surface_mod.f90
</r2008 r2011 1 MODULE urban_surface_mod2 3 1 !> @file urban_surface_mod.f90 4 2 !! … … 23 21 ! Current revisions: 24 22 !  25 ! 23 ! Major reformatting according to PALM coding standard (comments, blanks, 24 ! alphabetical ordering, etc.), 25 ! removed debug_prints, 26 ! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is 27 ! defined in MODULE control_parameters (modules.f90) to avoid circular 28 ! dependencies, 29 ! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed. 26 30 ! 27 31 ! Former revisions: … … 67 71 !> 68 72 !! 73 MODULE urban_surface_mod 74 75 USE arrays_3d, & 76 ONLY: zu, pt, pt_1, pt_2, p, ol, shf, ts, us, u, v, w, hyp, tend 77 78 USE cloud_parameters, & 79 ONLY: cp, r_d 69 80 70 81 USE constants, & 71 only: pi 82 ONLY: pi 83 84 USE control_parameters, & 85 ONLY: dz, topography, dt_3d, intermediate_timestep_count, & 86 initializing_actions, intermediate_timestep_count_max, & 87 simulated_time, end_time, timestep_scheme, tsc, & 88 coupling_char, io_blocks, io_group, message_string, & 89 time_since_reference_point, surface_pressure, & 90 g, pt_surface, large_scale_forcing, lsf_surf, & 91 time_do3d, dt_do3d, average_count_3d, urban_surface 92 93 USE cpulog, & 94 ONLY: cpu_log, log_point, log_point_s 95 96 USE grid_variables, & 97 ONLY: dx, dy, ddx, ddy, ddx2, ddy2 98 99 USE indices, & 100 ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, & 101 nysg, nzb_s_inner, nzb_s_outer, nzb, nzt, nbgp 102 103 USE, INTRINSIC :: iso_c_binding 72 104 73 105 USE kinds … … 75 107 USE pegrid 76 108 77 USE indices, & 78 only: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, & 79 nysg, nzb_s_inner, nzb_s_outer, nzb, nzt, nbgp 80 81 USE control_parameters, & 82 ONLY: dz, topography, dt_3d, intermediate_timestep_count, & 83 initializing_actions, intermediate_timestep_count_max, & 84 simulated_time, end_time, timestep_scheme, tsc, & 85 coupling_char, io_blocks, io_group, message_string, & 86 time_since_reference_point, surface_pressure, & 87 g, pt_surface, large_scale_forcing, lsf_surf, & 88 time_do3d, dt_do3d, average_count_3d 89 90 USE grid_variables, & 91 ONLY: dx, dy, ddx, ddy, ddx2, ddy2 92 93 USE arrays_3d, & 94 ONLY: zu, pt, pt_1, pt_2, p, ol, shf, ts, us, u, v, w, hyp, tend 95 96 USE cloud_parameters, & 97 ONLY: cp, r_d 109 USE plant_canopy_model_mod, & 110 ONLY: plant_canopy, pch_index, & 111 pc_heating_rate, lad_s 98 112 99 113 USE radiation_model_mod, & … … 102 116 sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon, & 103 117 force_radiation_call 104 105 USE plant_canopy_model_mod, &106 ONLY: plant_canopy, pch_index, &107 canopy_heat_flux, lad_s108 118 109 119 USE statistics, & 110 ONLY: hom, statistic_regions 111 112 USE cpulog, & 113 ONLY: cpu_log, log_point, log_point_s 114 115 ! USE ieee_arithmetic 116 USE, INTRINSIC :: iso_c_binding 120 ONLY: hom, statistic_regions 117 121 118 122 … … 120 124 IMPLICIT NONE 121 125 122 ! configuration parameters (they can be setup in PALM config) 123 LOGICAL :: urban_surface = .FALSE. !< switch for use of urban surface model 126 ! configuration parameters (they can be setup in PALM config) 124 127 LOGICAL :: split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation 125 128 !< (.F. in case the radiation model already does it) 126 129 LOGICAL :: usm_energy_balance_land = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs 127 130 LOGICAL :: usm_energy_balance_wall = .TRUE. !< flag parameter indicating wheather the energy balance is calculated for land and roofs 131 LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used 128 132 LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used 129 LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the wsm is used130 133 LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls 131 134 LOGICAL :: mrt_factors = .FALSE. !< whether to generate MRT factor files during init … … 133 136 LOGICAL :: read_svf_on_init = .FALSE. 134 137 LOGICAL :: usm_lad_rma = .TRUE. !< use MPI RMA to access LAD for raytracing (instead of global array) 135 LOGICAL :: debug_prints = .FALSE. !< print debug messages into process debug files136 138 137 139 INTEGER(iwp) :: nrefsteps = 0 !< number of reflection steps to perform … … 147 149 !< of r_a for horizontal surfaces > TODO 148 150 149 151 ! parameters of urban surface model 150 152 INTEGER(iwp), PARAMETER :: usm_version_len = 10 !< length of identification string of usm version 151 153 CHARACTER(usm_version_len), PARAMETER :: usm_version = 'USM v. 1.0' !< identification of version of binary svf and restart files … … 180 182 !< parameter but set in the code 181 183 182 ! indices and sizes of urban surface model 183 !ketelsen: INTEGER(iwp), DIMENSION(:,:,:,:), ALLOCATABLE :: gridsurf !< array of indices of the surfaces in local block in global vector of indices 184 !< of surfaces for grid coordinates (d,z,y,x) 185 !< d = 0 groud/roof, 1 south wall, 2 north wall, 3 west wall, 4 east wall 186 !< 5 free north border, 6 free south b., 7 free east b., 8 free west b. 184 ! indices and sizes of urban surface model 187 185 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surfl !< coordinates of ith local surface in local grid  surfl[:,k] = [d, z, y, x] 188 186 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: surf !< coordinates of ith surface in grid  surf[:,k] = [d, z, y, x] … … 211 209 212 210 213 211 ! type for calculation of svf 214 212 TYPE t_svf 215 213 INTEGER(iwp) :: isurflt !< … … 219 217 END TYPE 220 218 221 219 ! type for calculation of csf 222 220 TYPE t_csf 223 221 INTEGER(iwp) :: ip !< … … 230 228 END TYPE 231 229 232 230 ! arrays for calculation of svf and csf 233 231 TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array 234 232 TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array … … 240 238 INTEGER(iwp), PARAMETER :: gasize = 10000 !< initial size of growing arrays 241 239 242 240 ! arrays storing the values of USM 243 241 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of source and target surface for svf[isvf] 244 242 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors … … 264 262 REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_l !< local copy of rad_net (net radiation at surface) 265 263 266 264 ! arrays for time averages 267 265 REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l 268 266 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections … … 277 275 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface 278 276 279 277 ! block variables needed for calculation of the plant canopy model inside the urban surface model 280 278 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< z,y,x coordinates of ith local plant canopy box pcbl[:,i] = [z, y, x] 281 279 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< index of local pcb[z,y,x] … … 290 288 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing 291 289 292 290 ! radiation related arrays (it should be better in interface of radiation module of PALM 293 291 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation 294 292 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation 295 293 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation 296 294 297 298 299 295 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 296 ! anthropogenic heat sources 297 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 300 298 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2) 301 299 REAL(wp), DIMENSION(:), ALLOCATABLE :: aheatprof !< diurnal profile of anthropogenic heat 302 300 303 304 305 306 301 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 302 ! wall surface model 303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 304 ! wall surface model constants 307 305 INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched) 308 306 INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched) … … 316 314 REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10 degrees C) (K) 317 315 318 319 320 316 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 317 ! surface and material model variables for walls, ground, roofs 318 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 321 319 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_types !< array of types of wall parameters 322 320 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m) … … 339 337 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K) 340 338 341 339 ! Temporal tendencies for time stepping 342 340 REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency (K) 343 341 344 345 346 347 342 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 343 ! Energy balance variables 344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 345 ! parameters of the land, roof and wall surfaces 348 346 LOGICAL, DIMENSION(:), ALLOCATABLE :: isroof_surf !< is the surface the part of a roof 349 347 REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface 350 348 ! parameters of the wall surfaces 351 349 REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface !< heat capacity of the wall surface skin ( J mâ2 Kâ1 ) 352 350 REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface 353 351 REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 ) 354 352 355 353 ! parameters of the walls material 356 354 REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_wall !< thickness of the wall, roof and soil layers 357 355 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_wall !< volumetric heat capacity of the material ( J m3 K1 ) (= 2.19E6) … … 359 357 REAL(wp), DIMENSION(:), ALLOCATABLE :: roughness_wall !< roughness relative to concrete 360 358 361 ! 362 ! output wall heat flux arrays 359 ! output wall heat flux arrays 363 360 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf !< kinematic wall heat flux of sensible heat (needed for diffusion_s!<) 364 361 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb !< wall heat flux of sensible heat in wall normal direction … … 376 373 #endif 377 374 378 375 ! Wall temporal tendencies for time stepping 379 376 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array 380 377 381 382 378 ! Surface and material parameters classes (surface_type) 379 ! albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity 383 380 INTEGER(iwp) :: n_surface_types !< number of the wall type categories 384 381 INTEGER(iwp), PARAMETER :: n_surface_params = 8 !< number of parameters for each type of the wall … … 397 394 CHARACTER(len=*), PARAMETER :: svf_file_name='usm_svf' 398 395 399 ! interfaces of subroutines accessed from outside of this module 396 ! interfaces of subroutines accessed from outside of this module 397 INTERFACE usm_check_data_output 398 MODULE PROCEDURE usm_check_data_output 399 END INTERFACE usm_check_data_output 400 401 INTERFACE usm_check_parameters 402 MODULE PROCEDURE usm_check_parameters 403 END INTERFACE usm_check_parameters 404 405 INTERFACE usm_data_output_3d 406 MODULE PROCEDURE usm_data_output_3d 407 END INTERFACE usm_data_output_3d 408 409 INTERFACE usm_define_netcdf_grid 410 MODULE PROCEDURE usm_define_netcdf_grid 411 END INTERFACE usm_define_netcdf_grid 412 400 413 INTERFACE usm_init_urban_surface 401 414 MODULE PROCEDURE usm_init_urban_surface 402 415 END INTERFACE usm_init_urban_surface 403 416 417 INTERFACE usm_material_heat_model 418 MODULE PROCEDURE usm_material_heat_model 419 END INTERFACE usm_material_heat_model 420 421 INTERFACE usm_parin 422 MODULE PROCEDURE usm_parin 423 END INTERFACE usm_parin 424 404 425 INTERFACE usm_radiation 405 426 MODULE PROCEDURE usm_radiation 406 427 END INTERFACE usm_radiation 428 429 INTERFACE usm_read_restart_data 430 MODULE PROCEDURE usm_read_restart_data 431 END INTERFACE usm_read_restart_data 407 432 408 433 INTERFACE usm_surface_energy_balance 409 434 MODULE PROCEDURE usm_surface_energy_balance 410 435 END INTERFACE usm_surface_energy_balance 411 412 INTERFACE usm_ material_heat_model413 MODULE PROCEDURE usm_ material_heat_model414 END INTERFACE usm_ material_heat_model436 437 INTERFACE usm_swap_timelevel 438 MODULE PROCEDURE usm_swap_timelevel 439 END INTERFACE usm_swap_timelevel 415 440 416 441 INTERFACE usm_wall_heat_flux … … 419 444 END INTERFACE usm_wall_heat_flux 420 445 421 INTERFACE usm_swap_timelevel422 MODULE PROCEDURE usm_swap_timelevel423 END INTERFACE usm_swap_timelevel424 425 INTERFACE usm_check_data_output426 MODULE PROCEDURE usm_check_data_output427 END INTERFACE usm_check_data_output428 429 INTERFACE usm_check_parameters430 MODULE PROCEDURE usm_check_parameters431 END INTERFACE usm_check_parameters432 433 INTERFACE usm_data_output_3d434 MODULE PROCEDURE usm_data_output_3d435 END INTERFACE usm_data_output_3d436 437 INTERFACE usm_define_netcdf_grid438 MODULE PROCEDURE usm_define_netcdf_grid439 END INTERFACE usm_define_netcdf_grid440 441 INTERFACE usm_parin442 MODULE PROCEDURE usm_parin443 END INTERFACE usm_parin444 445 INTERFACE usm_read_restart_data446 MODULE PROCEDURE usm_read_restart_data447 END INTERFACE usm_read_restart_data448 449 446 INTERFACE usm_write_restart_data 450 447 MODULE PROCEDURE usm_write_restart_data … … 454 451 455 452 PRIVATE 456 !457 458 PUBLIC urban_surface, split_diffusion_radiation,&453 454 ! Public parameters, constants and initial values 455 PUBLIC split_diffusion_radiation, & 459 456 usm_anthropogenic_heat, usm_material_model, mrt_factors, & 460 457 usm_check_parameters, & … … 466 463 usm_data_output_3d, usm_define_netcdf_grid, usm_parin, & 467 464 usm_write_restart_data, & 468 nzub, nzut, ra_horiz_coef, usm_lad_rma, debug_prints,&465 nzub, nzut, ra_horiz_coef, usm_lad_rma, & 469 466 land_category, pedestrant_category, wall_category, roof_category, & 470 467 write_svf_on_init, read_svf_on_init … … 473 470 CONTAINS 474 471 475 476 477 !!478 ! Description:479 ! 480 !> Initialization of the urban surface model481 !!482 SUBROUTINE usm_init_urban_surface483 484 IMPLICIT NONE485 486 INTEGER(iwp) :: i, j, k, l ! running indices487 REAL(wp) :: c, d, tin, exn488 489 CALL cpu_log( log_point_s(78), 'usm_init', 'start' )490 ! surface forcing have to be disabled for LSF491 ! in case of enabled urban surface module492 IF ( large_scale_forcing ) THEN493 lsf_surf = .FALSE.494 ENDIF495 496 ! init anthropogenic sources of heat497 CALL usm_allocate_urban_surface()498 499 ! read the surface_types array somewhere500 CALL usm_read_urban_surface_types()501 502 ! init material heat model503 CALL usm_init_material_model()504 505 IF ( usm_anthropogenic_heat ) THEN506 ! init anthropogenic sources of heat (from transportation for now)507 CALL usm_read_anthropogenic_heat()508 ENDIF509 510 IF ( read_svf_on_init ) THEN511 ! read svf and svfsurf data from file512 WRITE(6,*) myid, 'Before read svf from file'513 FLUSH(6)514 CALL usm_read_svf_from_file()515 WRITE(6,*) myid, 'After read svf from file'516 FLUSH(6)517 ELSE518 ! calculate SFV and CSF519 WRITE(6,*) myid, 'Before calc svf'520 FLUSH(6)521 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'start' )522 CALL usm_calc_svf()523 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'stop' )524 WRITE(6,*) myid, 'After calc svf'525 FLUSH(6)526 ENDIF527 528 IF ( write_svf_on_init ) THEN529 ! write svf and svfsurf data to file530 WRITE(6,*) myid, 'Before write svf to file'531 FLUSH(6)532 CALL usm_write_svf_to_file()533 WRITE(6,*) myid, 'After write svf to file'534 FLUSH(6)535 ENDIF536 537 IF ( plant_canopy ) THEN538 ! gridpcbl was only necessary for initialization539 DEALLOCATE( gridpcbl )540 IF ( .NOT. ALLOCATED(canopy_heat_flux) ) THEN541 ! then canopy_heat_flux is allocated in init_plant_canopy542 ! in case of cthf /= 0 => we need to allocate it for our use here543 ALLOCATE( canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )544 ENDIF545 ENDIF546 547 ! Intitialization of the surface and wall/ground/roof temperature548 549 ! Initialization for restart runs550 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN551 552 ! restore data from restart file553 CALL usm_read_restart_data()554 ELSE555 556 ! Calculate initial surface temperature557 exn = ( surface_pressure / 1000.0_wp )**0.286_wp558 559 DO l = startenergy, endenergy560 k = surfl(iz,l)561 j = surfl(iy,l)562 i = surfl(ix,l)563 564 ! Initial surface temperature set from pt of adjacent gridbox565 t_surf(l) = pt(k,j,i) * exn566 ENDDO567 568 ! initial values for t_wall569 ! outer value is set to surface temperature570 ! inner value is set to wall_inner_temperature571 ! and profile is logaritmic (linear in nz)572 DO l = startenergy, endenergy573 IF ( isroof_surf(l) ) THEN574 tin = roof_inner_temperature575 ELSE IF ( surf(id,l)==iroof ) THEN576 tin = soil_inner_temperature577 ELSE578 tin = wall_inner_temperature579 ENDIF580 DO k = nzb_wall, nzt_wall+1581 c = REAL(knzb_wall,wp)/REAL(nzt_wall+1nzb_wall,wp)582 t_wall(k,:) = (1.0_wpc)*t_surf(:) + c*tin583 ENDDO584 ENDDO585 ENDIF586 587 !588 ! Possibly DO userdefined actions (e.g. define heterogeneous wall surface)589 CALL user_init_urban_surface590 591 ! initialize prognostic values for the first timestep592 t_surf_p = t_surf593 t_wall_p = t_wall594 595 ! Adjust radiative fluxes for urban surface at model start596 CALL usm_radiation597 598 CALL cpu_log( log_point_s(78), 'usm_init', 'stop' )599 600 601 END SUBROUTINE usm_init_urban_surface602 603 604 472 605 473 !! … … 620 488 621 489 622 !auxiliary vars490 ! auxiliary vars 623 491 ddxy2 = (/ddy2,ddy2,ddx2,ddx2/) !< 1/dx^2 or 1/dy^2 (in surface normal direction) 624 492 625 493 CALL location_message( '', .TRUE. ) 626 494 CALL location_message( ' allocation of needed arrays', .TRUE. ) 627 !find nzub, nzut, nzu495 ! find nzub, nzut, nzu 628 496 nzubl = minval(nzb_s_inner(nys:nyn,nxl:nxr)) 629 497 nzutl = maxval(nzb_s_inner(nys:nyn,nxl:nxr)) 630 498 nzubl = max(nzubl,nzb) 631 499 632 IF ( plant_canopy)THEN633 !allocate needed arrays500 IF ( plant_canopy ) THEN 501 ! allocate needed arrays 634 502 ALLOCATE( pct(nys:nyn,nxl:nxr) ) 635 503 ALLOCATE( pch(nys:nyn,nxl:nxr) ) 636 504 637 !calculate plant canopy height505 ! calculate plant canopy height 638 506 npcbl = 0 639 507 pct = 0.0_wp … … 642 510 DO j = nys, nyn 643 511 DO k = nzt+1, 0, 1 644 IF ( lad_s(k,j,i) /= 0.0_wp ) THEN645 !we are at the top of the pcs512 IF ( lad_s(k,j,i) /= 0.0_wp ) THEN 513 ! we are at the top of the pcs 646 514 pct(j,i) = k + nzb_s_inner(j,i) 647 515 pch(j,i) = k … … 654 522 655 523 nzutl = max(nzutl, maxval(pct)) 656 !code of plant canopy model uses parameter pch_index657 !we need to setup it here to right value658 !(pch_index, lad_s and other arrays in PCM are defined flat)524 ! code of plant canopy model uses parameter pch_index 525 ! we need to setup it here to right value 526 ! (pch_index, lad_s and other arrays in PCM are defined flat) 659 527 pch_index = maxval(pch) 660 528 661 529 prototype_lad = maxval(lad_s) * .9_wp !< better be *1.0 if lad is either 0 or maxval(lad) everywhere 662 530 IF ( prototype_lad <= 0._wp ) prototype_lad = .3_wp 663 WRITE(message_string, '(a,f6.3)') 'Precomputing effective box optical ' &664 // 'depth using prototype leaf area density = ', prototype_lad665 CALL message('usm_init_urban_surface', 'PA0520', 0, 0, 1, 6, 0)531 !WRITE(message_string, '(a,f6.3)') 'Precomputing effective box optical ' & 532 ! // 'depth using prototype leaf area density = ', prototype_lad 533 !CALL message('usm_init_urban_surface', 'PA0520', 0, 0, 1, 6, 0) 666 534 ENDIF 667 535 668 536 nzutl = min(nzutl+nzut_free, nzt) 669 537 … … 676 544 #endif 677 545 678 !global number of urban layers546 ! global number of urban layers 679 547 nzu = nzut  nzub + 1 680 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 681 WRITE(9,*) 'nzub= ', nzub, ' nzut= ', nzut, ' nzu= ', nzu 682 FLUSH(9) 683 ENDIF 684 685 ! allocate urban surfaces grid 686 !ketelsen: ALLOCATE(gridsurf(0:9,nzub:nzut,nys:nyn,nxl:nxr)) 687 !ketelsen: gridsurf = 0 688 689 ! calc number of surfaces in local proc 548 549 ! allocate urban surfaces grid 550 ! calc number of surfaces in local proc 690 551 CALL location_message( ' calculation of indices for surfaces', .TRUE. ) 691 552 nsurfl = 0 692 !calculate land surface and roof553 ! calculate land surface and roof 693 554 startland = nsurfl+1 694 555 nsurfl = nsurfl+(nxrnxl+1)*(nynnys+1) … … 696 557 nlands = endlandstartland+1 697 558 698 !calculation of the walls559 ! calculation of the walls 699 560 startwall = nsurfl+1 700 561 DO i = nxl, nxr 701 562 DO j = nys, nyn 702 !test for walls703 !(we don't use array flags because it isn't calculated in case of masking_method=.T.)563 ! test for walls 564 ! (we don't use array flags because it isn't calculated in case of masking_method=.T.) 704 565 DO ids = 1, 4 ! four wall directions 705 566 jr = min(max(jjdir(ids),0),ny) … … 712 573 nwalls = endwallstartwall+1 713 574 714 !range of energy balance surfaces575 ! range of energy balance surfaces 715 576 nenergy = 0 716 IF ( usm_energy_balance_land ) THEN577 IF ( usm_energy_balance_land ) THEN 717 578 startenergy = startland 718 579 nenergy = nenergy + nlands … … 720 581 startenergy = startwall 721 582 ENDIF 722 IF ( usm_energy_balance_wall ) THEN583 IF ( usm_energy_balance_wall ) THEN 723 584 endenergy = endwall 724 585 nenergy = nenergy + nwalls … … 727 588 ENDIF 728 589 729 730 !block of virtual surfaces731 732 !calculate sky surfaces590 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 591 ! block of virtual surfaces 592 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 593 ! calculate sky surfaces 733 594 startsky = nsurfl+1 734 595 nsurfl = nsurfl+(nxrnxl+1)*(nynnys+1) … … 736 597 nskys = endskystartsky+1 737 598 738 !border flags599 ! border flags 739 600 #if defined( __parallel ) 740 601 isborder = (/ north_border_pe, south_border_pe, right_border_pe, left_border_pe /) … … 742 603 isborder = (/.TRUE.,.TRUE.,.TRUE.,.TRUE./) 743 604 #endif 744 !fill array of the limits of the local domain borders605 ! fill array of the limits of the local domain borders 745 606 ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) ) 746 !calulation of the free borders of the domain607 ! calulation of the free borders of the domain 747 608 DO ids = 6,9 748 IF ( isborder(ids) ) THEN749 !free border of the domain in direction ids609 IF ( isborder(ids) ) THEN 610 ! free border of the domain in direction ids 750 611 DO i = ijdb(1,ids), ijdb(2,ids) 751 612 DO j = ijdb(3,ids), ijdb(4,ids) … … 757 618 ENDDO 758 619 759 !fill gridpcbl and pcbl760 IF ( plant_canopy ) THEN620 ! fill gridpcbl and pcbl 621 IF ( plant_canopy ) THEN 761 622 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 762 623 ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) ) … … 777 638 ENDIF 778 639 779 ! fill gridsurf andsurfl640 ! fill surfl 780 641 ALLOCATE(surfl(4,nsurfl)) 781 642 isurf = 0 782 643 783 !add land surfaces or roofs644 ! add land surfaces or roofs 784 645 DO i = nxl, nxr 785 646 DO j = nys, nyn 786 647 isurf = isurf + 1 787 648 k = nzb_s_inner(j,i)+1 788 !ketelsen: gridsurf(iroof,k,j,i) = isurf789 649 surfl(:,isurf) = (/iroof,k,j,i/) 790 650 ENDDO 791 651 ENDDO 792 652 793 !add walls653 ! add walls 794 654 DO i = nxl, nxr 795 655 DO j = nys, nyn … … 799 659 DO k = nzb_s_inner(j,i)+1, nzb_s_inner(jr,ir) 800 660 isurf = isurf + 1 801 !ketelsen: gridsurf(isouth,k,j,i) = isurf802 661 surfl(:,isurf) = (/ids,k,j,i/) 803 662 ENDDO … … 806 665 ENDDO 807 666 808 !add sky667 ! add sky 809 668 DO i = nxl, nxr 810 669 DO j = nys, nyn 811 670 isurf = isurf + 1 812 671 k = nzut 813 !ketelsen: gridsurf(isky,k,j,i) = isurf814 672 surfl(:,isurf) = (/isky,k,j,i/) 815 673 ENDDO 816 674 ENDDO 817 675 818 !calulation of the free borders of the domain676 ! calulation of the free borders of the domain 819 677 DO ids = 6,9 820 IF ( isborder(ids) ) THEN821 !free border of the domain in direction ids678 IF ( isborder(ids) ) THEN 679 ! free border of the domain in direction ids 822 680 DO i = ijdb(1,ids), ijdb(2,ids) 823 681 DO j = ijdb(3,ids), ijdb(4,ids) 824 682 DO k = max(nzb_s_inner(j,i),nzb_s_inner(jjdir(ids),iidir(ids)))+1, nzut 825 683 isurf = isurf + 1 826 !ketelsen: gridsurf(ids,k,j,i) = isurf827 684 surfl(:,isurf) = (/ids,k,j,i/) 828 685 ENDDO … … 832 689 ENDDO 833 690 834 !global array surf of indices of surfaces and displacement index array surfstart691 ! global array surf of indices of surfaces and displacement index array surfstart 835 692 ALLOCATE(nsurfs(0:numprocs1)) 836 693 … … 856 713 #endif 857 714 858 859 !allocation of the arrays for direct and diffusion radiation715 ! 716 ! allocation of the arrays for direct and diffusion radiation 860 717 CALL location_message( ' allocation of radiation arrays', .TRUE. ) 861 !rad_sw_in, rad_lw_in are computed in radiation model,862 !splitting of direct and diffusion part is done863 !in usm_calc_diffusion_radiation for now718 ! rad_sw_in, rad_lw_in are computed in radiation model, 719 ! splitting of direct and diffusion part is done 720 ! in usm_calc_diffusion_radiation for now 864 721 ALLOCATE( rad_sw_in_dir(nysg:nyng,nxlg:nxrg) ) 865 722 ALLOCATE( rad_sw_in_diff(nysg:nyng,nxlg:nxrg) ) 866 723 ALLOCATE( rad_lw_in_diff(nysg:nyng,nxlg:nxrg) ) 867 724 868 !allocate radiation arrays725 ! allocate radiation arrays 869 726 ALLOCATE( surfins(nsurfl) ) 870 727 ALLOCATE( surfinl(nsurfl) ) … … 883 740 ALLOCATE( rad_net_l(startenergy:endenergy) ) 884 741 885 !Wall surface model886 !allocate arrays for wall surface model and define pointers742 ! Wall surface model 743 ! allocate arrays for wall surface model and define pointers 887 744 888 !allocate array of wall types and wall parameters745 ! allocate array of wall types and wall parameters 889 746 ALLOCATE ( surface_types(startenergy:endenergy) ) 890 747 891 !broadband albedo of the land, roof and wall surface892 !for domain border and sky set artifically to 1.0893 !what allows us to calculate heat flux leaving over894 !side and top borders of the domain748 ! broadband albedo of the land, roof and wall surface 749 ! for domain border and sky set artifically to 1.0 750 ! what allows us to calculate heat flux leaving over 751 ! side and top borders of the domain 895 752 ALLOCATE ( albedo_surf(nsurfl) ) 896 753 albedo_surf = 1.0_wp 897 754 898 !wall and roof surface parameters755 ! wall and roof surface parameters 899 756 ALLOCATE ( isroof_surf(startenergy:endenergy) ) 900 757 ALLOCATE ( emiss_surf(startenergy:endenergy) ) … … 903 760 ALLOCATE ( roughness_wall(startenergy:endenergy) ) 904 761 905 !allocate wall and roof material parameters762 ! allocate wall and roof material parameters 906 763 ALLOCATE ( thickness_wall(startenergy:endenergy) ) 907 764 ALLOCATE ( lambda_h(nzb_wall:nzt_wall,startenergy:endenergy) ) 908 765 ALLOCATE ( rho_c_wall(nzb_wall:nzt_wall,startenergy:endenergy) ) 909 766 910 !allocate wall and roof layers sizes767 ! allocate wall and roof layers sizes 911 768 ALLOCATE ( zwn(nzb_wall:nzt_wall) ) 912 769 ALLOCATE ( dz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) ) … … 916 773 ALLOCATE ( zw(nzb_wall:nzt_wall, startenergy:endenergy) ) 917 774 918 !allocate wall and roof temperature arrays775 ! allocate wall and roof temperature arrays 919 776 #if defined( __nopointer ) 920 777 ALLOCATE ( t_surf(startenergy:endenergy) ) … … 928 785 ALLOCATE ( t_wall_2(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 929 786 930 !initial assignment of the pointers787 ! initial assignment of the pointers 931 788 t_wall => t_wall_1; t_wall_p => t_wall_2 932 789 t_surf => t_surf_1; t_surf_p => t_surf_2 933 790 #endif 934 791 935 !allocate intermediate timestep arrays792 ! allocate intermediate timestep arrays 936 793 ALLOCATE ( tt_surface_m(startenergy:endenergy) ) 937 794 ALLOCATE ( tt_wall_m(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 938 795 939 !allocate wall heat flux output array796 ! allocate wall heat flux output array 940 797 ALLOCATE ( wshf(startwall:endwall) ) 941 798 ALLOCATE ( wshf_eb(startenergy:endenergy) ) 942 799 ALLOCATE ( wghf_eb(startenergy:endenergy) ) 943 800 944 !set inital values for prognostic quantities801 ! set inital values for prognostic quantities 945 802 tt_surface_m = 0.0_wp 946 803 tt_wall_m = 0.0_wp … … 953 810 954 811 812 955 813 !! 956 814 ! Description: 957 815 !  958 !> Initialization of the wall surface model 959 !! 960 SUBROUTINE usm_init_material_model 816 !> Sum up and timeaverage urban surface output quantities as well as allocate 817 !> the array necessary for storing the average. 818 !! 819 SUBROUTINE usm_average_3d_data( mode, variable ) 961 820 962 821 IMPLICIT NONE 963 822 964 INTEGER(iwp) :: k, l ! running indices 965 966 CALL location_message( ' initialization of wall surface model', .TRUE. ) 967 968 ! Calculate wall grid spacings. 969 ! Temperature is defined at the center of the wall layers, 970 ! whereas gradients/fluxes are defined at the edges (_stag) 971 DO l = nzb_wall, nzt_wall 972 zwn(l) = zwn_default(l) 823 CHARACTER (len=*), INTENT(IN) :: mode 824 CHARACTER (len=*), INTENT(IN) :: variable 825 826 INTEGER(iwp) :: i, j, k, l, ids, iwl,istat 827 CHARACTER (len=20) :: var, surfid 828 INTEGER(iwp), PARAMETER :: nd = 5 829 CHARACTER(len=6), DIMENSION(0:nd1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 830 831 ! find the real name of the variable 832 var = TRIM(variable) 833 DO i = 0, nd1 834 k = len(TRIM(var)) 835 j = len(TRIM(dirname(i))) 836 IF ( var(kj+1:k) == dirname(i) ) THEN 837 ids = i 838 var = var(:kj) 839 EXIT 840 ENDIF 973 841 ENDDO 974 975 ! apply for all particular wall grids 976 DO l = startenergy, endenergy 977 zw(:,l) = zwn(:) * thickness_wall(l) 978 dz_wall(nzb_wall,l) = zw(nzb_wall,l) 979 DO k = nzb_wall+1, nzt_wall 980 dz_wall(k,l) = zw(k,l)  zw(k1,l) 981 ENDDO 842 IF ( ids == 1 ) THEN 843 var = TRIM(variable) 844 ENDIF 845 IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN 846 ! wall layers 847 READ(var(12:12), '(I1)', iostat=istat ) iwl 848 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 849 var = var(1:10) 850 ELSE 851 ! wrong wall layer index 852 RETURN 853 ENDIF 854 ENDIF 855 856 IF ( mode == 'allocate' ) THEN 982 857 983 dz_wall(nzt_wall+1,l) = dz_wall(nzt_wall,l) 984 985 DO k = nzb_wall, nzt_wall1 986 dz_wall_stag(k,l) = 0.5 * (dz_wall(k+1,l) + dz_wall(k,l)) 987 ENDDO 988 dz_wall_stag(nzt_wall,l) = dz_wall(nzt_wall,l) 858 SELECT CASE ( TRIM( var ) ) 859 860 CASE ( 'usm_radnet' ) 861 ! array of complete radiation balance 862 IF ( .NOT. ALLOCATED(rad_net_av) ) THEN 863 ALLOCATE( rad_net_av(startenergy:endenergy) ) 864 rad_net_av = 0.0_wp 865 ENDIF 866 867 CASE ( 'usm_rad_insw' ) 868 ! array of sw radiation falling to surface after ith reflection 869 IF ( .NOT. ALLOCATED(surfinsw_av) ) THEN 870 ALLOCATE( surfinsw_av(startenergy:endenergy) ) 871 surfinsw_av = 0.0_wp 872 ENDIF 873 874 CASE ( 'usm_rad_inlw' ) 875 ! array of lw radiation falling to surface after ith reflection 876 IF ( .NOT. ALLOCATED(surfinlw_av) ) THEN 877 ALLOCATE( surfinlw_av(startenergy:endenergy) ) 878 surfinlw_av = 0.0_wp 879 ENDIF 880 881 CASE ( 'usm_rad_inswdir' ) 882 ! array of direct sw radiation falling to surface from sun 883 IF ( .NOT. ALLOCATED(surfinswdir_av) ) THEN 884 ALLOCATE( surfinswdir_av(startenergy:endenergy) ) 885 surfinswdir_av = 0.0_wp 886 ENDIF 887 888 CASE ( 'usm_rad_inswdif' ) 889 ! array of difusion sw radiation falling to surface from sky and borders of the domain 890 IF ( .NOT. ALLOCATED(surfinswdif_av) ) THEN 891 ALLOCATE( surfinswdif_av(startenergy:endenergy) ) 892 surfinswdif_av = 0.0_wp 893 ENDIF 894 895 CASE ( 'usm_rad_inswref' ) 896 ! array of sw radiation falling to surface from reflections 897 IF ( .NOT. ALLOCATED(surfinswref_av) ) THEN 898 ALLOCATE( surfinswref_av(startenergy:endenergy) ) 899 surfinswref_av = 0.0_wp 900 ENDIF 901 902 CASE ( 'usm_rad_inlwdif' ) 903 ! array of sw radiation falling to surface after ith reflection 904 IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN 905 ALLOCATE( surfinlwdif_av(startenergy:endenergy) ) 906 surfinlwdif_av = 0.0_wp 907 ENDIF 908 909 CASE ( 'usm_rad_inlwref' ) 910 ! array of lw radiation falling to surface from reflections 911 IF ( .NOT. ALLOCATED(surfinlwref_av) ) THEN 912 ALLOCATE( surfinlwref_av(startenergy:endenergy) ) 913 surfinlwref_av = 0.0_wp 914 ENDIF 915 916 CASE ( 'usm_rad_outsw' ) 917 ! array of sw radiation emitted from surface after ith reflection 918 IF ( .NOT. ALLOCATED(surfoutsw_av) ) THEN 919 ALLOCATE( surfoutsw_av(startenergy:endenergy) ) 920 surfoutsw_av = 0.0_wp 921 ENDIF 922 923 CASE ( 'usm_rad_outlw' ) 924 ! array of lw radiation emitted from surface after ith reflection 925 IF ( .NOT. ALLOCATED(surfoutlw_av) ) THEN 926 ALLOCATE( surfoutlw_av(startenergy:endenergy) ) 927 surfoutlw_av = 0.0_wp 928 ENDIF 929 930 CASE ( 'usm_rad_hf' ) 931 ! array of heat flux from radiation for surfaces after ith reflection 932 IF ( .NOT. ALLOCATED(surfhf_av) ) THEN 933 ALLOCATE( surfhf_av(startenergy:endenergy) ) 934 surfhf_av = 0.0_wp 935 ENDIF 936 937 CASE ( 'usm_wshf' ) 938 ! array of sensible heat flux from surfaces 939 ! land surfaces 940 IF ( .NOT. ALLOCATED(wshf_eb_av) ) THEN 941 ALLOCATE( wshf_eb_av(startenergy:endenergy) ) 942 wshf_eb_av = 0.0_wp 943 ENDIF 944 945 CASE ( 'usm_wghf' ) 946 ! array of heat flux from ground (wall, roof, land) 947 IF ( .NOT. ALLOCATED(wghf_eb_av) ) THEN 948 ALLOCATE( wghf_eb_av(startenergy:endenergy) ) 949 wghf_eb_av = 0.0_wp 950 ENDIF 951 952 CASE ( 'usm_t_surf' ) 953 ! surface temperature for surfaces 954 IF ( .NOT. ALLOCATED(t_surf_av) ) THEN 955 ALLOCATE( t_surf_av(startenergy:endenergy) ) 956 t_surf_av = 0.0_wp 957 ENDIF 958 959 CASE ( 'usm_t_wall' ) 960 ! wall temperature for iwl layer of walls and land 961 IF ( .NOT. ALLOCATED(t_wall_av) ) THEN 962 ALLOCATE( t_wall_av(nzb_wall:nzt_wall,startenergy:endenergy) ) 963 t_wall_av = 0.0_wp 964 ENDIF 965 966 CASE DEFAULT 967 CONTINUE 968 969 END SELECT 970 971 ELSEIF ( mode == 'sum' ) THEN 972 973 SELECT CASE ( TRIM( var ) ) 974 975 CASE ( 'usm_radnet' ) 976 ! array of complete radiation balance 977 DO l = startenergy, endenergy 978 IF ( surfl(id,l) == ids ) THEN 979 rad_net_av(l) = rad_net_av(l) + rad_net_l(l) 980 ENDIF 981 ENDDO 982 983 CASE ( 'usm_rad_insw' ) 984 ! array of sw radiation falling to surface after ith reflection 985 DO l = startenergy, endenergy 986 IF ( surfl(id,l) == ids ) THEN 987 surfinsw_av(l) = surfinsw_av(l) + surfinsw(l) 988 ENDIF 989 ENDDO 990 991 CASE ( 'usm_rad_inlw' ) 992 ! array of lw radiation falling to surface after ith reflection 993 DO l = startenergy, endenergy 994 IF ( surfl(id,l) == ids ) THEN 995 surfinlw_av(l) = surfinlw_av(l) + surfinlw(l) 996 ENDIF 997 ENDDO 998 999 CASE ( 'usm_rad_inswdir' ) 1000 ! array of direct sw radiation falling to surface from sun 1001 DO l = startenergy, endenergy 1002 IF ( surfl(id,l) == ids ) THEN 1003 surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l) 1004 ENDIF 1005 ENDDO 1006 1007 CASE ( 'usm_rad_inswdif' ) 1008 ! array of difusion sw radiation falling to surface from sky and borders of the domain 1009 DO l = startenergy, endenergy 1010 IF ( surfl(id,l) == ids ) THEN 1011 surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l) 1012 ENDIF 1013 ENDDO 1014 1015 CASE ( 'usm_rad_inswref' ) 1016 ! array of sw radiation falling to surface from reflections 1017 DO l = startenergy, endenergy 1018 IF ( surfl(id,l) == ids ) THEN 1019 surfinswref_av(l) = surfinswref_av(l) + surfinsw(l)  & 1020 surfinswdir(l)  surfinswdif(l) 1021 ENDIF 1022 ENDDO 1023 1024 CASE ( 'usm_rad_inlwdif' ) 1025 ! array of sw radiation falling to surface after ith reflection 1026 DO l = startenergy, endenergy 1027 IF ( surfl(id,l) == ids ) THEN 1028 surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l) 1029 ENDIF 1030 ENDDO 1031 1032 CASE ( 'usm_rad_inlwref' ) 1033 ! array of lw radiation falling to surface from reflections 1034 DO l = startenergy, endenergy 1035 IF ( surfl(id,l) == ids ) THEN 1036 surfinlwref_av(l) = surfinlwref_av(l) + & 1037 surfinlw(l)  surfinlwdif(l) 1038 ENDIF 1039 ENDDO 1040 1041 CASE ( 'usm_rad_outsw' ) 1042 ! array of sw radiation emitted from surface after ith reflection 1043 DO l = startenergy, endenergy 1044 IF ( surfl(id,l) == ids ) THEN 1045 surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l) 1046 ENDIF 1047 ENDDO 1048 1049 CASE ( 'usm_rad_outlw' ) 1050 ! array of lw radiation emitted from surface after ith reflection 1051 DO l = startenergy, endenergy 1052 IF ( surfl(id,l) == ids ) THEN 1053 surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l) 1054 ENDIF 1055 ENDDO 1056 1057 CASE ( 'usm_rad_hf' ) 1058 ! array of heat flux from radiation for surfaces after ith reflection 1059 DO l = startenergy, endenergy 1060 IF ( surfl(id,l) == ids ) THEN 1061 surfhf_av(l) = surfhf_av(l) + surfhf(l) 1062 ENDIF 1063 ENDDO 1064 1065 CASE ( 'usm_wshf' ) 1066 ! array of sensible heat flux from surfaces (land, roof, wall) 1067 DO l = startenergy, endenergy 1068 IF ( surfl(id,l) == ids ) THEN 1069 wshf_eb_av(l) = wshf_eb_av(l) + wshf_eb(l) 1070 ENDIF 1071 ENDDO 1072 1073 CASE ( 'usm_wghf' ) 1074 ! array of heat flux from ground (wall, roof, land) 1075 DO l = startenergy, endenergy 1076 IF ( surfl(id,l) == ids ) THEN 1077 wghf_eb_av(l) = wghf_eb_av(l) + wghf_eb(l) 1078 ENDIF 1079 ENDDO 1080 1081 CASE ( 'usm_t_surf' ) 1082 ! surface temperature for surfaces 1083 DO l = startenergy, endenergy 1084 IF ( surfl(id,l) == ids ) THEN 1085 t_surf_av(l) = t_surf_av(l) + t_surf(l) 1086 ENDIF 1087 ENDDO 1088 1089 CASE ( 'usm_t_wall' ) 1090 ! wall temperature for iwl layer of walls and land 1091 DO l = startenergy, endenergy 1092 IF ( surfl(id,l) == ids ) THEN 1093 t_wall_av(iwl, l) = t_wall_av(iwl,l) + t_wall(iwl, l) 1094 ENDIF 1095 ENDDO 1096 1097 CASE DEFAULT 1098 CONTINUE 1099 1100 END SELECT 1101 1102 ELSEIF ( mode == 'average' ) THEN 1103 1104 SELECT CASE ( TRIM( var ) ) 1105 1106 CASE ( 'usm_radnet' ) 1107 ! array of complete radiation balance 1108 DO l = startenergy, endenergy 1109 IF ( surfl(id,l) == ids ) THEN 1110 rad_net_av(l) = rad_net_av(l) / REAL( average_count_3d, kind=wp ) 1111 ENDIF 1112 ENDDO 1113 1114 CASE ( 'usm_rad_insw' ) 1115 ! array of sw radiation falling to surface after ith reflection 1116 DO l = startenergy, endenergy 1117 IF ( surfl(id,l) == ids ) THEN 1118 surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp ) 1119 ENDIF 1120 ENDDO 1121 1122 CASE ( 'usm_rad_inlw' ) 1123 ! array of lw radiation falling to surface after ith reflection 1124 DO l = startenergy, endenergy 1125 IF ( surfl(id,l) == ids ) THEN 1126 surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp ) 1127 ENDIF 1128 ENDDO 1129 1130 CASE ( 'usm_rad_inswdir' ) 1131 ! array of direct sw radiation falling to surface from sun 1132 DO l = startenergy, endenergy 1133 IF ( surfl(id,l) == ids ) THEN 1134 surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp ) 1135 ENDIF 1136 ENDDO 1137 1138 CASE ( 'usm_rad_inswdif' ) 1139 ! array of difusion sw radiation falling to surface from sky and borders of the domain 1140 DO l = startenergy, endenergy 1141 IF ( surfl(id,l) == ids ) THEN 1142 surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp ) 1143 ENDIF 1144 ENDDO 1145 1146 CASE ( 'usm_rad_inswref' ) 1147 ! array of sw radiation falling to surface from reflections 1148 DO l = startenergy, endenergy 1149 IF ( surfl(id,l) == ids ) THEN 1150 surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp ) 1151 ENDIF 1152 ENDDO 1153 1154 CASE ( 'usm_rad_inlwdif' ) 1155 ! array of sw radiation falling to surface after ith reflection 1156 DO l = startenergy, endenergy 1157 IF ( surfl(id,l) == ids ) THEN 1158 surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp ) 1159 ENDIF 1160 ENDDO 1161 1162 CASE ( 'usm_rad_inlwref' ) 1163 ! array of lw radiation falling to surface from reflections 1164 DO l = startenergy, endenergy 1165 IF ( surfl(id,l) == ids ) THEN 1166 surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp ) 1167 ENDIF 1168 ENDDO 1169 1170 CASE ( 'usm_rad_outsw' ) 1171 ! array of sw radiation emitted from surface after ith reflection 1172 DO l = startenergy, endenergy 1173 IF ( surfl(id,l) == ids ) THEN 1174 surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp ) 1175 ENDIF 1176 ENDDO 1177 1178 CASE ( 'usm_rad_outlw' ) 1179 ! array of lw radiation emitted from surface after ith reflection 1180 DO l = startenergy, endenergy 1181 IF ( surfl(id,l) == ids ) THEN 1182 surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp ) 1183 ENDIF 1184 ENDDO 1185 1186 CASE ( 'usm_rad_hf' ) 1187 ! array of heat flux from radiation for surfaces after ith reflection 1188 DO l = startenergy, endenergy 1189 IF ( surfl(id,l) == ids ) THEN 1190 surfhf_av(l) = surfhf_av(l) / REAL( average_count_3d, kind=wp ) 1191 ENDIF 1192 ENDDO 1193 1194 CASE ( 'usm_wshf' ) 1195 ! array of sensible heat flux from surfaces (land, roof, wall) 1196 DO l = startenergy, endenergy 1197 IF ( surfl(id,l) == ids ) THEN 1198 wshf_eb_av(l) = wshf_eb_av(l) / REAL( average_count_3d, kind=wp ) 1199 ENDIF 1200 ENDDO 1201 1202 CASE ( 'usm_wghf' ) 1203 ! array of heat flux from ground (wall, roof, land) 1204 DO l = startenergy, endenergy 1205 IF ( surfl(id,l) == ids ) THEN 1206 wghf_eb_av(l) = wghf_eb_av(l) / REAL( average_count_3d, kind=wp ) 1207 ENDIF 1208 ENDDO 1209 1210 CASE ( 'usm_t_surf' ) 1211 ! surface temperature for surfaces 1212 DO l = startenergy, endenergy 1213 IF ( surfl(id,l) == ids ) THEN 1214 t_surf_av(l) = t_surf_av(l) / REAL( average_count_3d, kind=wp ) 1215 ENDIF 1216 ENDDO 1217 1218 CASE ( 'usm_t_wall' ) 1219 ! wall temperature for iwl layer of walls and land 1220 DO l = startenergy, endenergy 1221 IF ( surfl(id,l) == ids ) THEN 1222 t_wall_av(iwl, l) = t_wall_av(iwl,l) / REAL( average_count_3d, kind=wp ) 1223 ENDIF 1224 ENDDO 1225 1226 END SELECT 1227 1228 ENDIF 1229 1230 END SUBROUTINE usm_average_3d_data 1231 1232 1233 !! 1234 !> Calculates radiation absorbed by box with given size and LAD. 1235 !> 1236 !> Simulates resol**2 rays (by equally spacing a bounding horizontal square 1237 !> conatining all possible rays that would cross the box) and calculates 1238 !> average transparency per ray. Returns fraction of absorbed radiation flux 1239 !> and area for which this fraction is effective. 1240 !! 1241 PURE SUBROUTINE usm_box_absorb(boxsize, resol, dens, uvec, area, absorb) 1242 IMPLICIT NONE 1243 1244 REAL(wp), DIMENSION(3), INTENT(in) :: & 1245 boxsize, & !< z, y, x size of box in m 1246 uvec !< z, y, x unit vector of incoming flux 1247 INTEGER(iwp), INTENT(in) :: & 1248 resol !< No. of rays in x and y dimensions 1249 REAL(wp), INTENT(in) :: & 1250 dens !< box density (e.g. Leaf Area Density) 1251 REAL(wp), INTENT(out) :: & 1252 area, & !< horizontal area for flux absorbtion 1253 absorb !< fraction of absorbed flux 1254 REAL(wp) :: & 1255 xshift, yshift, & 1256 xmin, xmax, ymin, ymax, & 1257 xorig, yorig, & 1258 dx1, dy1, dz1, dx2, dy2, dz2, & 1259 crdist, & 1260 transp 1261 INTEGER(iwp) :: & 1262 i, j 1263 1264 xshift = uvec(3) / uvec(1) * boxsize(1) 1265 xmin = min(0._wp, xshift) 1266 xmax = boxsize(3) + max(0._wp, xshift) 1267 yshift = uvec(2) / uvec(1) * boxsize(1) 1268 ymin = min(0._wp, yshift) 1269 ymax = boxsize(2) + max(0._wp, yshift) 1270 1271 transp = 0._wp 1272 DO i = 1, resol 1273 xorig = xmin + (xmaxxmin) * (i.5_wp) / resol 1274 DO j = 1, resol 1275 yorig = ymin + (ymaxymin) * (j.5_wp) / resol 1276 1277 dz1 = 0._wp 1278 dz2 = boxsize(1)/uvec(1) 1279 1280 IF ( uvec(2) > 0._wp ) THEN 1281 dy1 = yorig / uvec(2) !< crossing with y=0 1282 dy2 = (boxsize(2)yorig) / uvec(2) !< crossing with y=boxsize(2) 1283 ELSE IF ( uvec(2) < 0._wp ) THEN 1284 dy1 = (boxsize(2)yorig) / uvec(2) !< crossing with y=boxsize(2) 1285 dy2 = yorig / uvec(2) !< crossing with y=0 1286 ELSE !uvec(2)==0 1287 dy1 = huge(1._wp) 1288 dy2 = huge(1._wp) 1289 ENDIF 1290 1291 IF ( uvec(3) > 0._wp ) THEN 1292 dx1 = xorig / uvec(3) !< crossing with x=0 1293 dx2 = (boxsize(3)xorig) / uvec(3) !< crossing with x=boxsize(3) 1294 ELSE IF ( uvec(3) < 0._wp ) THEN 1295 dx1 = (boxsize(3)xorig) / uvec(3) !< crossing with x=boxsize(3) 1296 dx2 = xorig / uvec(3) !< crossing with x=0 1297 ELSE !uvec(1)==0 1298 dx1 = huge(1._wp) 1299 dx2 = huge(1._wp) 1300 ENDIF 1301 1302 crdist = max(0._wp, (min(dz2, dy2, dx2)  max(dz1, dy1, dx1))) 1303 transp = transp + exp(ext_coef * dens * crdist) 1304 ENDDO 989 1305 ENDDO 990 991 ddz_wall = 1.0_wp / dz_wall 992 ddz_wall_stag = 1.0_wp / dz_wall_stag 993 994 CALL location_message( ' wall structures filed out', .TRUE. ) 995 996 CALL location_message( ' initialization of wall surface model finished', .TRUE. ) 997 998 END SUBROUTINE usm_init_material_model 999 1306 transp = transp / resol**2 1307 area = (boxsize(3)+xshift)*(boxsize(2)+yshift) 1308 absorb = 1._wp  transp 1309 1310 END SUBROUTINE usm_box_absorb 1311 1312 1313 !! 1314 ! Description: 1315 !  1316 !> This subroutine splits direct and diffusion dw radiation 1317 !> It sould not be called in case the radiation model already does it 1318 !> It follows <CITATION> 1319 !! 1320 SUBROUTINE usm_calc_diffusion_radiation 1321 1322 REAL(wp), PARAMETER :: sol_const = 1367.0_wp !< solar conbstant 1323 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 1324 INTEGER(iwp) :: i, j 1325 REAL(wp), PARAMETER :: year_seconds = 86400._wp * 365._wp 1326 REAL(wp) :: year_angle !< angle 1327 REAL(wp) :: etr !< extraterestrial radiation 1328 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 1329 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 1330 REAL(wp) :: clearnessIndex !< clearness index 1331 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 1332 1333 1334 ! Calculate current day and time based on the initial values and simulation time 1335 year_angle = ((day_init*86400) + time_utc_init+time_since_reference_point) & 1336 / year_seconds * 2.0_wp * pi 1337 1338 etr = sol_const * (1.00011_wp + & 1339 0.034221_wp * cos(year_angle) + & 1340 0.001280_wp * sin(year_angle) + & 1341 0.000719_wp * cos(2.0_wp * year_angle) + & 1342 0.000077_wp * sin(2.0_wp * year_angle)) 1343 1344 ! 1345 ! Under a very low angle, we keep extraterestrial radiation at 1346 ! the last small value, therefore the clearness index will be pushed 1347 ! towards 0 while keeping full continuity. 1348 ! 1349 IF ( zenith(0) <= lowest_solarUp ) THEN 1350 corrected_solarUp = lowest_solarUp 1351 ELSE 1352 corrected_solarUp = zenith(0) 1353 ENDIF 1354 1355 horizontalETR = etr * corrected_solarUp 1356 1357 DO i = nxlg, nxrg 1358 DO j = nysg, nyng 1359 clearnessIndex = rad_sw_in(0,j,i) / horizontalETR 1360 diff_frac = 1.0_wp / (1.0_wp + exp(5.0033_wp + 8.6025_wp * clearnessIndex)) 1361 rad_sw_in_diff(j,i) = rad_sw_in(0,j,i) * diff_frac 1362 rad_sw_in_dir(j,i) = rad_sw_in(0,j,i) * (1.0_wp  diff_frac) 1363 rad_lw_in_diff(j,i) = rad_lw_in(0,j,i) 1364 ENDDO 1365 ENDDO 1366 1367 END SUBROUTINE usm_calc_diffusion_radiation 1368 1000 1369 1001 1370 !! … … 1029 1398 INTEGER(kind=MPI_ADDRESS_KIND) :: size_lad_rma 1030 1399 1031 !calculation of the SVF1400 ! calculation of the SVF 1032 1401 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 1033 1402 1034 !precalculate face areas for different face directions using normal vector1403 ! precalculate face areas for different face directions using normal vector 1035 1404 DO d = 0, 9 1036 1405 facearea(d) = 1._wp 1037 IF ( idir(d) == 0) facearea(d) = facearea(d) * dx1038 IF ( jdir(d) == 0) facearea(d) = facearea(d) * dy1039 IF ( kdir(d) == 0) facearea(d) = facearea(d) * dz1406 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 1407 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 1408 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 1040 1409 ENDDO 1041 1410 1042 !initialize variables and temporary arrays for calculation of svf and csf1411 ! initialize variables and temporary arrays for calculation of svf and csf 1043 1412 nsvfl = 0 1044 1413 ncsfl = 0 … … 1047 1416 ALLOCATE( asvf1(nsvfla) ) 1048 1417 asvf => asvf1 1049 IF ( plant_canopy ) THEN1418 IF ( plant_canopy ) THEN 1050 1419 ncsfla = gasize 1051 1420 mcsf = 1 … … 1054 1423 ENDIF 1055 1424 1056 !initialize temporary terrain and plant canopy height arrays (global 2D array!)1425 ! initialize temporary terrain and plant canopy height arrays (global 2D array!) 1057 1426 ALLOCATE( nzterr(0:(nx+1)*(ny+1)1) ) 1058 1427 #if defined( __parallel ) … … 1065 1434 nzterr = RESHAPE( nzb_s_inner(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) ) 1066 1435 #endif 1067 IF ( plant_canopy)THEN1436 IF ( plant_canopy ) THEN 1068 1437 ALLOCATE( plantt(0:(nx+1)*(ny+1)1) ) 1069 1438 #if defined( __parallel ) … … 1075 1444 DEALLOCATE(planthl) 1076 1445 1077 IF ( usm_lad_rma ) THEN1446 IF ( usm_lad_rma ) THEN 1078 1447 CALL MPI_Info_create(minfo, ierr) 1079 1448 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr) … … 1082 1451 CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr) 1083 1452 1084 !Allocate and initialize the MPI RMA window1085 !must be in accordance with allocation of lad_s in plant_canopy_model1086 !optimization of memory should be done1087 !Argument X of function c_sizeof(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now1453 ! Allocate and initialize the MPI RMA window 1454 ! must be in accordance with allocation of lad_s in plant_canopy_model 1455 ! optimization of memory should be done 1456 ! Argument X of function c_sizeof(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now 1088 1457 size_lad_rma = c_sizeof(1.0_wp)*nnx*nny*nzu 1089 1458 CALL MPI_Win_allocate(size_lad_rma, c_sizeof(1.0_wp), minfo, comm2d, & 1090 1459 lad_s_rma_p, win_lad, ierr) 1091 IF ( debug_prints ) THEN1092 WRITE(9,*) 'MPI_Win_allocate', myid, ierr1093 FLUSH(9)1094 ENDIF1095 1460 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzu, nny, nnx /)) 1096 !IF ( debug_prints ) THEN1097 ! WRITE(9,*) 'lad_s', 'tsize=',c_sizeof(lad_s(0,0,0)), 'size=',c_sizeof(lad_s), &1098 ! 'nnz=', nnz, lbound(lad_s, 1), ubound(lad_s, 1), nzb, nzt, &1099 ! 'nny=', nny, lbound(lad_s, 2), ubound(lad_s, 2), nys, nyn, &1100 ! 'nnx=', nnx, lbound(lad_s, 3), ubound(lad_s, 3), nxl, nxr1101 !ENDIF1102 1461 usm_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:) 1103 1462 ELSE … … 1117 1476 1118 1477 #if defined( __parallel ) 1119 IF ( usm_lad_rma ) THEN1478 IF ( usm_lad_rma ) THEN 1120 1479 CALL MPI_Info_free(minfo, ierr) 1121 1480 CALL MPI_Win_lock_all(0, win_lad, ierr) 1122 IF ( debug_prints ) THEN1123 WRITE(9,*) 'MPI_Win_lock_all', myid, ierr1124 FLUSH(9)1125 ENDIF1126 1481 ELSE 1127 1482 ALLOCATE( usm_lad_g(0:(nx+1)*(ny+1)*nzu1) ) … … 1132 1487 ENDIF 1133 1488 1134 IF ( mrt_factors ) THEN 1135 IF ( debug_prints ) THEN 1136 WRITE(9, *) myid, 'mrt factors start' 1137 FLUSH(9) 1138 ENDIF 1489 IF ( mrt_factors ) THEN 1139 1490 OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', & 1140 1491 action='READ', status='OLD', form='FORMATTED', err=524) … … 1144 1495 DO 1145 1496 READ(153, *, end=526, err=524) imrtt, i, j, k 1146 IF ( i < nxl .OR.i > nxr &1147 .OR. j < nys .OR.j > nyn ) CYCLE1497 IF ( i < nxl .OR. i > nxr & 1498 .OR. j < nys .OR. j > nyn ) CYCLE 1148 1499 tx = REAL(i) 1149 1500 ty = REAL(j) … … 1151 1502 1152 1503 DO isurfs = 1, nsurf 1153 IF ( .NOT. usm_facing(i, j, k, 1, &1504 IF ( .NOT. usm_facing(i, j, k, 1, & 1154 1505 surf(ix, isurfs), surf(iy, isurfs), & 1155 surf(iz, isurfs), surf(id, isurfs)) ) THEN1506 surf(iz, isurfs), surf(id, isurfs)) ) THEN 1156 1507 CYCLE 1157 1508 ENDIF … … 1162 1513 sx = REAL(surf(ix, isurfs), wp)  0.5_wp * idir(sd) 1163 1514 1164 !unit vector source > target1515 ! unit vector source > target 1165 1516 uv = (/ (tzsz)*dz, (tysy)*dy, (txsx)*dx /) 1166 1517 sqdist = SUM(uv(:)**2) 1167 1518 uv = uv / SQRT(sqdist) 1168 1519 1169 !irradiance factor  see svf. Here we consider that target face is always normal,1170 !i.e. the second dot product equals 11520 ! irradiance factor  see svf. Here we consider that target face is always normal, 1521 ! i.e. the second dot product equals 1 1171 1522 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & 1172 1523 / (pi * sqdist) * facearea(sd) 1173 1524 1174 !raytrace while not creating any canopy sink factors1525 ! raytrace while not creating any canopy sink factors 1175 1526 CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, 1._wp, .FALSE., & 1176 1527 visible, transparency, win_lad) 1177 IF ( .NOT. visible ) CYCLE1528 IF ( .NOT. visible ) CYCLE 1178 1529 1179 1530 !rsvf = rirrf * transparency … … 1197 1548 526 CLOSE(153) 1198 1549 CLOSE(154) 1199 IF ( debug_prints ) THEN1200 WRITE(9, *) myid, 'mrt factors finished', ncsfl !should be always 01201 FLUSH(9)1202 ENDIF1203 1550 ENDIF !< mrt_factors 1204 1551 1205 1552 1206 1553 DO isurflt = 1, nsurfl 1207 !determine face centers1554 ! determine face centers 1208 1555 td = surfl(id, isurflt) 1209 IF ( td >= isky .AND. .NOT. plant_canopy) CYCLE1556 IF ( td >= isky .AND. .NOT. plant_canopy ) CYCLE 1210 1557 tz = REAL(surfl(iz, isurflt), wp)  0.5_wp * kdir(td) 1211 1558 ty = REAL(surfl(iy, isurflt), wp)  0.5_wp * jdir(td) 1212 1559 tx = REAL(surfl(ix, isurflt), wp)  0.5_wp * idir(td) 1213 IF ( debug_prints ) THEN1214 WRITE(9,'(i3,a,i2,2i4, 4i6,3f7.2,4i6,3f7.2)') myid, 'svft1', &1215 td, isurflt, nsurfl, &1216 surfl(ix, isurflt), surfl(iy, isurflt), &1217 surfl(iz, isurflt), surfl(id, isurflt), tx,ty,tz1218 FLUSH(9)1219 ENDIF1220 1560 DO isurfs = 1, nsurf 1221 IF ( .NOT. usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), &1561 IF ( .NOT. usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 1222 1562 surfl(iz, isurflt), surfl(id, isurflt), & 1223 1563 surf(ix, isurfs), surf(iy, isurfs), & 1224 surf(iz, isurfs), surf(id, isurfs)) ) THEN1564 surf(iz, isurfs), surf(id, isurfs)) ) THEN 1225 1565 CYCLE 1226 1566 ENDIF … … 1231 1571 sx = REAL(surf(ix, isurfs), wp)  0.5_wp * idir(sd) 1232 1572 1233 !unit vector source > target1573 ! unit vector source > target 1234 1574 uv = (/ (tzsz)*dz, (tysy)*dy, (txsx)*dx /) 1235 1575 sqdist = SUM(uv(:)**2) 1236 1576 uv = uv / SQRT(sqdist) 1237 1577 1238 !irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area1578 ! irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 1239 1579 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 1240 * dot_product((/ kdir(td), jdir(td), idir(td) /), uv) & ! cosine of target normal and reverse direction1580 * dot_product((/ kdir(td), jdir(td), idir(td) /), uv) & ! cosine of target normal and reverse direction 1241 1581 / (pi * sqdist) & ! square of distance between centers 1242 1582 * facearea(sd) 1243 1583 1244 !raytrace + process plant canopy sinks within1584 ! raytrace + process plant canopy sinks within 1245 1585 CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, facearea(td), .TRUE., & 1246 1586 visible, transparency, win_lad) 1247 1587 1248 IF ( .NOT. visible ) CYCLE1588 IF ( .NOT. visible ) CYCLE 1249 1589 IF ( td >= isky ) CYCLE !< we calculated these only for raytracing 1250 1590 !< to find plant canopy sinks, we don't need svf for them 1251 1252 !rsvf = rirrf * transparency 1253 1254 ! write to the svf array 1591 ! rsvf = rirrf * transparency 1592 1593 ! write to the svf array 1255 1594 nsvfl = nsvfl + 1 1256 !check dimmension of asvf array and enlarge it if needed1257 IF ( nsvfla < nsvfl ) THEN1595 ! check dimmension of asvf array and enlarge it if needed 1596 IF ( nsvfla < nsvfl ) THEN 1258 1597 k = nsvfla * 2 1259 IF ( debug_prints ) THEN 1260 WRITE(9,*) 'New dimmension of asvf array set to ', k 1261 FLUSH(9) 1262 ENDIF 1263 IF ( msvf == 0 ) THEN 1598 IF ( msvf == 0 ) THEN 1264 1599 msvf = 1 1265 1600 ALLOCATE( asvf1(k) ) … … 1276 1611 nsvfla = k 1277 1612 ENDIF 1278 !write svf values into the array1613 ! write svf values into the array 1279 1614 asvf(nsvfl)%isurflt = isurflt 1280 1615 asvf(nsvfl)%isurfs = isurfs … … 1284 1619 ENDDO 1285 1620 1286 WRITE(6,*) myid, 'waiting for completion of svf and csf calculation in all processes' 1287 FLUSH(6) 1621 CALL location_message( ' waiting for completion of SVF and CSF calculation in all processes', .TRUE. ) 1288 1622 1289 1623 #if defined( __parallel ) 1290 IF ( plant_canopy)THEN1291 IF ( usm_lad_rma ) THEN1624 IF ( plant_canopy ) THEN 1625 IF ( usm_lad_rma ) THEN 1292 1626 CALL MPI_Win_flush_all(win_lad, ierr) 1293 IF ( debug_prints ) THEN 1294 WRITE(9,*) 'MPI_Win_flush_all', myid, ierr 1295 FLUSH(9) 1296 ENDIF 1297 ! unlock MPI window 1627 ! unlock MPI window 1298 1628 CALL MPI_Win_unlock_all(win_lad, ierr) 1299 IF ( debug_prints ) THEN 1300 WRITE(9,*) 'MPI_Win_unlock_all', myid, ierr 1301 FLUSH(9) 1302 ENDIF 1303 ! free MPI window 1629 ! free MPI window 1304 1630 CALL MPI_Win_free(win_lad, ierr) 1305 IF ( debug_prints ) THEN1306 WRITE(9,*) 'MPI_Win_free', myid, ierr1307 FLUSH(9)1308 ENDIF1309 1631 ELSE 1310 1632 DEALLOCATE(usm_lad) … … 1316 1638 #endif 1317 1639 1318 !deallocate temporary global arrays1640 ! deallocate temporary global arrays 1319 1641 IF ( ALLOCATED(nzterr) ) DEALLOCATE(nzterr) 1320 1642 IF ( ALLOCATED(plantt) ) DEALLOCATE(plantt) 1321 1643 1322 IF ( debug_prints ) THEN 1323 WRITE(9,*) myid, 'writing svf finished, nsvfl', nsvfl, ', ncsfl', ncsfl 1324 FLUSH(9) 1325 ENDIF 1326 ! sort svf ( a version of quicksort ) 1327 IF ( debug_prints ) THEN 1328 WRITE(9,*) myid, 'start svf sort' 1329 FLUSH(9) 1330 ENDIF 1644 ! sort svf ( a version of quicksort ) 1331 1645 CALL quicksort_svf(asvf,1,nsvfl) 1332 IF ( debug_prints ) THEN1333 WRITE(9,*) myid, 'sort svf finished'1334 FLUSH(9)1335 ENDIF1336 1646 1337 1647 npcsfl = 0 1338 IF ( plant_canopy)THEN1339 !sort and merge csf for the last time, keeping the array size to minimum1340 CALL merge_and_grow_csf(1)1648 IF ( plant_canopy ) THEN 1649 ! sort and merge csf for the last time, keeping the array size to minimum 1650 CALL usm_merge_and_grow_csf(1) 1341 1651 1342 IF ( debug_prints ) THEN 1343 WRITE(9,*) myid, 'distribute CSF into processors' 1344 FLUSH(9) 1345 ENDIF 1346 1347 ! aggregate csb among processors 1348 ! allocate necessary arrays 1652 ! aggregate csb among processors 1653 ! allocate necessary arrays 1349 1654 ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) ) 1350 1655 ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) ) … … 1354 1659 ALLOCATE( dpcsflt(0:numprocs1) ) 1355 1660 1356 !fill out arrays of csf values and1357 !arrays of number of elements and displacements1358 !for particular precessors1661 ! fill out arrays of csf values and 1662 ! arrays of number of elements and displacements 1663 ! for particular precessors 1359 1664 icsflt = 0 1360 1665 dcsflt = 0 … … 1364 1669 DO kcsf = 1, ncsfl 1365 1670 j = j+1 1366 IF ( acsf(kcsf)%ip /= ip ) THEN1367 !new block of the processor1368 !number of elements of previous block1369 IF ( ip>=0) icsflt(ip) = j1671 IF ( acsf(kcsf)%ip /= ip ) THEN 1672 ! new block of the processor 1673 ! number of elements of previous block 1674 IF ( ip>=0) icsflt(ip) = j 1370 1675 d = d+j 1371 !blank blocks1676 ! blank blocks 1372 1677 DO jp = ip+1, acsf(kcsf)%ip1 1373 !number of elements is zero, displacement is equal to previous1678 ! number of elements is zero, displacement is equal to previous 1374 1679 icsflt(jp) = 0 1375 1680 dcsflt(jp) = d 1376 1681 ENDDO 1377 !the actual block1682 ! the actual block 1378 1683 ip = acsf(kcsf)%ip 1379 1684 dcsflt(ip) = d 1380 1685 j = 0 1381 1686 ENDIF 1382 !fill out real values of rsvf, rtransp1687 ! fill out real values of rsvf, rtransp 1383 1688 csflt(1,kcsf) = acsf(kcsf)%rsvf 1384 1689 csflt(2,kcsf) = acsf(kcsf)%rtransp 1385 !fill out integer values of itz,ity,itx,isurfs1690 ! fill out integer values of itz,ity,itx,isurfs 1386 1691 kcsflt(1,kcsf) = acsf(kcsf)%itz 1387 1692 kcsflt(2,kcsf) = acsf(kcsf)%ity … … 1389 1694 kcsflt(4,kcsf) = acsf(kcsf)%isurfs 1390 1695 ENDDO 1391 !last blank blocks at the end of array1696 ! last blank blocks at the end of array 1392 1697 j = j+1 1393 IF ( ip>=0) icsflt(ip) = j1698 IF ( ip>=0 ) icsflt(ip) = j 1394 1699 d = d+j 1395 1700 DO jp = ip+1, numprocs1 1396 !number of elements is zero, displacement is equal to previous1701 ! number of elements is zero, displacement is equal to previous 1397 1702 icsflt(jp) = 0 1398 1703 dcsflt(jp) = d 1399 1704 ENDDO 1400 1705 1401 !deallocate temporary acsf array1402 !DEALLOCATE(acsf)  ifort has a problem with deallocation of allocatable target1403 !via pointing pointer  we need to test original targets1404 IF ( ALLOCATED(acsf1) ) THEN1706 ! deallocate temporary acsf array 1707 ! DEALLOCATE(acsf)  ifort has a problem with deallocation of allocatable target 1708 ! via pointing pointer  we need to test original targets 1709 IF ( ALLOCATED(acsf1) ) THEN 1405 1710 DEALLOCATE(acsf1) 1406 1711 ENDIF 1407 IF ( ALLOCATED(acsf2) ) THEN1712 IF ( ALLOCATED(acsf2) ) THEN 1408 1713 DEALLOCATE(acsf2) 1409 1714 ENDIF 1410 1715 1411 1716 #if defined( __parallel ) 1412 !scatter and gather the number of elements to and from all processor1413 !and calculate displacements1717 ! scatter and gather the number of elements to and from all processor 1718 ! and calculate displacements 1414 1719 CALL mpi_alltoall(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr) 1415 1720 … … 1421 1726 ENDDO 1422 1727 1423 !exchange csf fields between processors1728 ! exchange csf fields between processors 1424 1729 ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) ) 1425 1730 ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) ) … … 1429 1734 kpcsflt, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr) 1430 1735 1431 IF ( debug_prints ) THEN1432 WRITE(9,*) myid, 'distribution CSF into processors finished'1433 FLUSH(9)1434 ENDIF1435 1736 #else 1436 1737 npcsfl = ncsfl … … 1441 1742 #endif 1442 1743 1443 !deallocate temporary arrays1744 ! deallocate temporary arrays 1444 1745 DEALLOCATE( csflt ) 1445 1746 DEALLOCATE( kcsflt ) … … 1449 1750 DEALLOCATE( dpcsflt ) 1450 1751 1451 ! sort csf ( a version of quicksort ) 1452 IF ( debug_prints ) THEN 1453 WRITE(9,*) myid, 'start csf sort2' 1454 FLUSH(9) 1455 ENDIF 1752 ! sort csf ( a version of quicksort ) 1456 1753 CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl) 1457 IF ( debug_prints ) THEN 1458 WRITE(9,*) myid, 'sort2 csf finished' 1459 FLUSH(9) 1460 ENDIF 1461 1462 ! aggregate canopy sink factor records with identical box & source 1463 ! againg across all values from all processors 1464 IF ( npcsfl > 0 ) THEN 1465 IF ( debug_prints ) THEN 1466 WRITE(9,*) myid, 'csf merge of all values with', npcsfl, 'boxes' 1467 FLUSH(9) 1468 ENDIF 1754 1755 ! aggregate canopy sink factor records with identical box & source 1756 ! againg across all values from all processors 1757 IF ( npcsfl > 0 ) THEN 1469 1758 icsf = 1 !< reading index 1470 1759 kcsf = 1 !< writing index 1471 1760 DO while (icsf < npcsfl) 1472 !here kpcsf(kcsf) already has values from kpcsf(icsf)1473 IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1) .AND.&1474 kpcsflt(2,icsf) == kpcsflt(2,icsf+1) .AND.&1475 kpcsflt(1,icsf) == kpcsflt(1,icsf+1) .AND.&1476 kpcsflt(4,icsf) == kpcsflt(4,icsf+1) ) THEN1477 !We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray1478 !probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that1479 !might mean that the traced beam passes longer through the canopy box.1480 IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) ) THEN1761 ! here kpcsf(kcsf) already has values from kpcsf(icsf) 1762 IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1) .AND. & 1763 kpcsflt(2,icsf) == kpcsflt(2,icsf+1) .AND. & 1764 kpcsflt(1,icsf) == kpcsflt(1,icsf+1) .AND. & 1765 kpcsflt(4,icsf) == kpcsflt(4,icsf+1) ) THEN 1766 ! We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray 1767 ! probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that 1768 ! might mean that the traced beam passes longer through the canopy box. 1769 IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) ) THEN 1481 1770 pcsflt(2,kcsf) = pcsflt(2,icsf+1) 1482 1771 ENDIF 1483 1772 pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1) 1484 1773 1485 !advance reading index, keep writing index1774 ! advance reading index, keep writing index 1486 1775 icsf = icsf + 1 1487 1776 ELSE 1488 !not identical, just advance and copy1777 ! not identical, just advance and copy 1489 1778 icsf = icsf + 1 1490 1779 kcsf = kcsf + 1 … … 1493 1782 ENDIF 1494 1783 ENDDO 1495 !last written item is now also the last item in valid part of array1784 ! last written item is now also the last item in valid part of array 1496 1785 npcsfl = kcsf 1497 IF ( debug_prints ) THEN1498 WRITE(9,*) myid, 'csf merge of all values completed,', npcsfl, 'boxes remaining'1499 FLUSH(9)1500 ENDIF1501 1786 ENDIF 1502 1787 … … 1514 1799 svfsum = 0._wp 1515 1800 DO isvf = 1, nsvfl 1516 !normalize svf per target face1517 IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN1518 IF ( isurflt_prev /= 1 .AND. svfsum /= 0._wp )THEN1519 !TODO detect and log when normalization differs too much from 11801 ! normalize svf per target face 1802 IF ( asvf(ksvf)%isurflt /= isurflt_prev ) THEN 1803 IF ( isurflt_prev /= 1 .AND. svfsum /= 0._wp ) THEN 1804 ! TODO detect and log when normalization differs too much from 1 1520 1805 svf(1, isvf_surflt:isvf1) = svf(1, isvf_surflt:isvf1) / svfsum 1521 1806 ENDIF … … 1530 1815 svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /) 1531 1816 1532 !next element1817 ! next element 1533 1818 ksvf = ksvf + 1 1534 1819 ENDDO 1535 1820 1536 IF ( isurflt_prev /= 1 .AND. svfsum /= 0._wp )THEN1537 !TODO detect and log when normalization differs too much from 11821 IF ( isurflt_prev /= 1 .AND. svfsum /= 0._wp ) THEN 1822 ! TODO detect and log when normalization differs too much from 1 1538 1823 svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum 1539 1824 ENDIF 1540 1825 1541 !deallocate temporary asvf array1542 !DEALLOCATE(asvf)  ifort has a problem with deallocation of allocatable target1543 !via pointing pointer  we need to test original targets1544 IF ( ALLOCATED(asvf1) ) THEN1826 ! deallocate temporary asvf array 1827 ! DEALLOCATE(asvf)  ifort has a problem with deallocation of allocatable target 1828 ! via pointing pointer  we need to test original targets 1829 IF ( ALLOCATED(asvf1) ) THEN 1545 1830 DEALLOCATE(asvf1) 1546 1831 ENDIF 1547 IF ( ALLOCATED(asvf2) ) THEN1832 IF ( ALLOCATED(asvf2) ) THEN 1548 1833 DEALLOCATE(asvf2) 1549 1834 ENDIF 1550 1835 1551 IF ( plant_canopy)THEN1836 IF ( plant_canopy ) THEN 1552 1837 CALL location_message( ' calculation of the complete CSF part of the array', .TRUE. ) 1553 IF ( npcsfl > 0 ) THEN1838 IF ( npcsfl > 0 ) THEN 1554 1839 DO isvf = 1, npcsfl 1555 1840 svf(:,nsvfl+isvf) = pcsflt(:,isvf) … … 1559 1844 ENDIF 1560 1845 1561 !deallocation of temporary arrays1846 ! deallocation of temporary arrays 1562 1847 DEALLOCATE( pcsflt ) 1563 1848 DEALLOCATE( kpcsflt ) … … 1573 1858 1574 1859 END SUBROUTINE usm_calc_svf 1575 1860 1861 1862 !! 1863 ! 1864 ! Description: 1865 !  1866 !> Subroutine checks variables and assigns units. 1867 !> It is caaled out from subroutine check_parameters. 1868 !! 1869 SUBROUTINE usm_check_data_output( variable, unit ) 1870 1871 IMPLICIT NONE 1872 1873 CHARACTER (len=*),INTENT(IN) :: variable !: 1874 CHARACTER (len=*),INTENT(OUT) :: unit !: 1875 1876 CHARACTER (len=20) :: var 1877 1878 var = TRIM(variable) 1879 IF ( var(1:11)=='usm_radnet_' .OR. var(1:13)=='usm_rad_insw_' .OR. & 1880 var(1:13)=='usm_rad_inlw_' .OR. var(1:16)=='usm_rad_inswdir_' .OR. & 1881 var(1:16)=='usm_rad_inswdif_' .OR. var(1:16)=='usm_rad_inswref_' .OR. & 1882 var(1:16)=='usm_rad_inlwdif_' .OR. var(1:16)=='usm_rad_inlwref_' .OR. & 1883 var(1:14)=='usm_rad_outsw_' .OR. var(1:14)=='usm_rad_outlw_' .OR. & 1884 var(1:11)=='usm_rad_hf_' .OR. & 1885 var(1:9) =='usm_wshf_' .OR. var(1:9)=='usm_wghf_' ) THEN 1886 unit = 'W/m2' 1887 ELSE IF ( var(1:10) =='usm_t_surf' .OR. var(1:10) =='usm_t_wall' ) THEN 1888 unit = 'K' 1889 ELSE IF ( var(1:9) =='usm_surfz' .OR. var(1:7) =='usm_svf' .OR. & 1890 var(1:7) =='usm_dif' .OR. var(1:11) =='usm_surfcat' .OR. & 1891 var(1:11) =='usm_surfalb' .OR. var(1:12) =='usm_surfemis') THEN 1892 unit = '1' 1893 ELSE IF ( plant_canopy .AND. var(1:7) =='usm_lad' ) THEN 1894 unit = 'm2/m3' 1895 ELSE IF ( plant_canopy .AND. var(1:14) == 'usm_canopy_khf' ) THEN 1896 unit = 'K/s' 1897 ELSE 1898 unit = 'illegal' 1899 ENDIF 1900 1901 END SUBROUTINE usm_check_data_output 1902 1903 1904 !! 1905 ! Description: 1906 !  1907 !> Check parameters routine for urban surface model 1908 !! 1909 SUBROUTINE usm_check_parameters 1910 1911 USE control_parameters, & 1912 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, & 1913 lsf_surf, topography 1914 1915 ! 1916 ! Dirichlet boundary conditions are required as the surface fluxes are 1917 ! calculated from the temperature/humidity gradients in the urban surface 1918 ! model 1919 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 1920 message_string = 'urban surface model requires setting of '// & 1921 'bc_pt_b = "dirichlet" and '// & 1922 'bc_q_b = "dirichlet"' 1923 CALL message( 'check_parameters', 'PA0590', 1, 2, 0, 6, 0 ) 1924 ENDIF 1925 1926 IF ( .NOT. constant_flux_layer ) THEN 1927 message_string = 'urban surface model requires '// & 1928 'constant_flux_layer = .T.' 1929 CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 1930 ENDIF 1931 ! 1932 ! Surface forcing has to be disabled for LSF in case of enabled 1933 ! urban surface module 1934 IF ( large_scale_forcing ) THEN 1935 lsf_surf = .FALSE. 1936 ENDIF 1937 ! 1938 ! Topography 1939 IF ( topography == 'flat' ) THEN 1940 message_string = 'topography /= "flat" is required '// & 1941 'when using the urban surface model' 1942 CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 ) 1943 ENDIF 1944 1945 1946 END SUBROUTINE usm_check_parameters 1947 1948 1949 !! 1950 ! 1951 ! Description: 1952 !  1953 !> Output of the 3Darrays in netCDF and/or AVS format 1954 !> for variables of urban_surface model. 1955 !> It resorts the urban surface module output quantities from surf style 1956 !> indexing into temporary 3D array with indices (i,j,k). 1957 !> It is called from subroutine data_output_3d. 1958 !! 1959 SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 1960 1961 IMPLICIT NONE 1962 1963 INTEGER(iwp), INTENT(IN) :: av !< 1964 CHARACTER (len=*), INTENT(IN) :: variable !< 1965 INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0) 1966 INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 1967 LOGICAL, INTENT(OUT) :: found !< 1968 REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) :: local_pf !< sp  it has to correspond to module data_output_3d 1969 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: temp_pf !< temp array for urban surface output procedure 1970 1971 CHARACTER (len=20) :: var, surfid 1972 INTEGER(iwp), PARAMETER :: nd = 5 1973 CHARACTER(len=6), DIMENSION(0:nd1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 1974 INTEGER(iwp), DIMENSION(0:nd1), PARAMETER :: dirint = (/ iroof, isouth, inorth, iwest, ieast /) 1975 INTEGER(iwp), DIMENSION(0:nd1) :: dirstart 1976 INTEGER(iwp), DIMENSION(0:nd1) :: dirend 1977 INTEGER(iwp) :: ids,isurf,isvf,isurfs,isurflt 1978 INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat 1979 1980 dirstart = (/ startland, startwall, startwall, startwall, startwall /) 1981 dirend = (/ endland, endwall, endwall, endwall, endwall /) 1982 1983 found = .TRUE. 1984 temp_pf = 1._wp 1985 1986 ids = 1 1987 var = TRIM(variable) 1988 DO i = 0, nd1 1989 k = len(TRIM(var)) 1990 j = len(TRIM(dirname(i))) 1991 IF ( var(kj+1:k) == dirname(i) ) THEN 1992 ids = i 1993 var = var(:kj) 1994 EXIT 1995 ENDIF 1996 ENDDO 1997 IF ( ids == 1 ) THEN 1998 var = TRIM(variable) 1999 ENDIF 2000 IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN 2001 ! wall layers 2002 READ(var(12:12), '(I1)', iostat=istat ) iwl 2003 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN 2004 var = var(1:10) 2005 ENDIF 2006 ENDIF 2007 IF ( (var(1:8) == 'usm_svf_' .OR. var(1:8) == 'usm_dif_') .AND. len(TRIM(var)) >= 13 ) THEN 2008 ! svf values to particular surface 2009 surfid = var(9:) 2010 i = index(surfid,'_') 2011 j = index(surfid(i+1:),'_') 2012 READ(surfid(1:i1),*, iostat=istat ) is 2013 IF ( istat == 0 ) THEN 2014 READ(surfid(i+1:i+j1),*, iostat=istat ) js 2015 ENDIF 2016 IF ( istat == 0 ) THEN 2017 READ(surfid(i+j+1:),*, iostat=istat ) ks 2018 ENDIF 2019 IF ( istat == 0 ) THEN 2020 var = var(1:7) 2021 ENDIF 2022 ENDIF 2023 2024 SELECT CASE ( TRIM(var) ) 2025 2026 CASE ( 'usm_surfz' ) 2027 ! array of lw radiation falling to local surface after ith reflection 2028 DO isurf = dirstart(ids), dirend(ids) 2029 IF ( surfl(id,isurf) == ids ) THEN 2030 IF ( surfl(id,isurf) == iroof ) THEN 2031 temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) = & 2032 max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)), & 2033 REAL(surfl(iz,isurf),wp)) 2034 ELSE 2035 temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) = & 2036 max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)), & 2037 REAL(surfl(iz,isurf),wp)+1.0_wp) 2038 ENDIF 2039 ENDIF 2040 ENDDO 2041 2042 CASE ( 'usm_surfcat' ) 2043 ! surface category 2044 DO isurf = dirstart(ids), dirend(ids) 2045 IF ( surfl(id,isurf) == ids ) THEN 2046 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surface_types(isurf) 2047 ENDIF 2048 ENDDO 2049 2050 CASE ( 'usm_surfalb' ) 2051 ! surface albedo 2052 DO isurf = dirstart(ids), dirend(ids) 2053 IF ( surfl(id,isurf) == ids ) THEN 2054 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = albedo_surf(isurf) 2055 ENDIF 2056 ENDDO 2057 2058 CASE ( 'usm_surfemis' ) 2059 ! surface albedo 2060 DO isurf = dirstart(ids), dirend(ids) 2061 IF ( surfl(id,isurf) == ids ) THEN 2062 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = emiss_surf(isurf) 2063 ENDIF 2064 ENDDO 2065 2066 CASE ( 'usm_svf', 'usm_dif' ) 2067 ! shape view factors or iradiance factors to selected surface 2068 IF ( TRIM(var)=='usm_svf' ) THEN 2069 k = 1 2070 ELSE 2071 k = 2 2072 ENDIF 2073 DO isvf = 1, nsvfl 2074 isurflt = svfsurf(1, isvf) 2075 isurfs = svfsurf(2, isvf) 2076 2077 IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. & 2078 surf(iz,isurfs) == ks .AND. surf(id,isurfs) == ids ) THEN 2079 ! correct source surface 2080 temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf) 2081 ENDIF 2082 ENDDO 2083 2084 CASE ( 'usm_radnet' ) 2085 ! array of complete radiation balance 2086 DO isurf = dirstart(ids), dirend(ids) 2087 IF ( surfl(id,isurf) == ids ) THEN 2088 IF ( av == 0 ) THEN 2089 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_l(isurf) 2090 ELSE 2091 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_av(isurf) 2092 ENDIF 2093 ENDIF 2094 ENDDO 2095 2096 CASE ( 'usm_rad_insw' ) 2097 ! array of sw radiation falling to surface after ith reflection 2098 DO isurf = dirstart(ids), dirend(ids) 2099 IF ( surfl(id,isurf) == ids ) THEN 2100 IF ( av == 0 ) THEN 2101 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf) 2102 ELSE 2103 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf) 2104 ENDIF 2105 ENDIF 2106 ENDDO 2107 2108 CASE ( 'usm_rad_inlw' ) 2109 ! array of lw radiation falling to surface after ith reflection 2110 DO isurf = dirstart(ids), dirend(ids) 2111 IF ( surfl(id,isurf) == ids ) THEN 2112 IF ( av == 0 ) THEN 2113 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) 2114 ELSE 2115 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf) 2116 ENDIF 2117 ENDIF 2118 ENDDO 2119 2120 CASE ( 'usm_rad_inswdir' ) 2121 ! array of direct sw radiation falling to surface from sun 2122 DO isurf = dirstart(ids), dirend(ids) 2123 IF ( surfl(id,isurf) == ids ) THEN 2124 IF ( av == 0 ) THEN 2125 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf) 2126 ELSE 2127 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf) 2128 ENDIF 2129 ENDIF 2130 ENDDO 2131 2132 CASE ( 'usm_rad_inswdif' ) 2133 ! array of difusion sw radiation falling to surface from sky and borders of the domain 2134 DO isurf = dirstart(ids), dirend(ids) 2135 IF ( surfl(id,isurf) == ids ) THEN 2136 IF ( av == 0 ) THEN 2137 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf) 2138 ELSE 2139 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf) 2140 ENDIF 2141 ENDIF 2142 ENDDO 2143 2144 CASE ( 'usm_rad_inswref' ) 2145 ! array of sw radiation falling to surface from reflections 2146 DO isurf = dirstart(ids), dirend(ids) 2147 IF ( surfl(id,isurf) == ids ) THEN 2148 IF ( av == 0 ) THEN 2149 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = & 2150 surfinsw(isurf)  surfinswdir(isurf)  surfinswdif(isurf) 2151 ELSE 2152 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf) 2153 ENDIF 2154 ENDIF 2155 ENDDO 2156 2157 CASE ( 'usm_rad_inlwdif' ) 2158 ! array of sw radiation falling to surface after ith reflection 2159 DO isurf = dirstart(ids), dirend(ids) 2160 IF ( surfl(id,isurf) == ids ) THEN 2161 IF ( av == 0 ) THEN 2162 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf) 2163 ELSE 2164 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf) 2165 ENDIF 2166 ENDIF 2167 ENDDO 2168 2169 CASE ( 'usm_rad_inlwref' ) 2170 ! array of lw radiation falling to surface from reflections 2171 DO isurf = dirstart(ids), dirend(ids) 2172 IF ( surfl(id,isurf) == ids ) THEN 2173 IF ( av == 0 ) THEN 2174 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)  surfinlwdif(isurf) 2175 ELSE 2176 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf) 2177 ENDIF 2178 ENDIF 2179 ENDDO 2180 2181 CASE ( 'usm_rad_outsw' ) 2182 ! array of sw radiation emitted from surface after ith reflection 2183 DO isurf = dirstart(ids), dirend(ids) 2184 IF ( surfl(id,isurf) == ids ) THEN 2185 IF ( av == 0 ) THEN 2186 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf) 2187 ELSE 2188 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf) 2189 ENDIF 2190 ENDIF 2191 ENDDO 2192 2193 CASE ( 'usm_rad_outlw' ) 2194 ! array of lw radiation emitted from surface after ith reflection 2195 DO isurf = dirstart(ids), dirend(ids) 2196 IF ( surfl(id,isurf) == ids ) THEN 2197 IF ( av == 0 ) THEN 2198 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf) 2199 ELSE 2200 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf) 2201 ENDIF 2202 ENDIF 2203 ENDDO 2204 2205 CASE ( 'usm_rad_hf' ) 2206 ! array of heat flux from radiation for surfaces after all reflections 2207 DO isurf = dirstart(ids), dirend(ids) 2208 IF ( surfl(id,isurf) == ids ) THEN 2209 IF ( av == 0 ) THEN 2210 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf(isurf) 2211 ELSE 2212 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf_av(isurf) 2213 ENDIF 2214 ENDIF 2215 ENDDO 2216 2217 CASE ( 'usm_wshf' ) 2218 ! array of sensible heat flux from surfaces 2219 ! horizontal surfaces 2220 DO isurf = dirstart(ids), dirend(ids) 2221 IF ( surfl(id,isurf) == ids ) THEN 2222 IF ( av == 0 ) THEN 2223 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb(isurf) 2224 ELSE 2225 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb_av(isurf) 2226 ENDIF 2227 ENDIF 2228 ENDDO 2229 2230 CASE ( 'usm_wghf' ) 2231 ! array of heat flux from ground (land, wall, roof) 2232 DO isurf = dirstart(ids), dirend(ids) 2233 IF ( surfl(id,isurf) == ids ) THEN 2234 IF ( av == 0 ) THEN 2235 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb(isurf) 2236 ELSE 2237 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb_av(isurf) 2238 ENDIF 2239 ENDIF 2240 ENDDO 2241 2242 CASE ( 'usm_t_surf' ) 2243 ! surface temperature for surfaces 2244 DO isurf = max(startenergy,dirstart(ids)), min(endenergy,dirend(ids)) 2245 IF ( surfl(id,isurf) == ids ) THEN 2246 IF ( av == 0 ) THEN 2247 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf(isurf) 2248 ELSE 2249 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf_av(isurf) 2250 ENDIF 2251 ENDIF 2252 ENDDO 2253 2254 CASE ( 'usm_t_wall' ) 2255 ! wall temperature for iwl layer of walls and land 2256 DO isurf = dirstart(ids), dirend(ids) 2257 IF ( surfl(id,isurf) == ids ) THEN 2258 IF ( av == 0 ) THEN 2259 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall(iwl,isurf) 2260 ELSE 2261 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall_av(iwl,isurf) 2262 ENDIF 2263 ENDIF 2264 ENDDO 2265 2266 CASE ( 'usm_lad' ) 2267 ! leaf area density 2268 DO i = nxl, nxr 2269 DO j = nys, nyn 2270 DO k = nzb_s_inner(j,i), nzut 2271 temp_pf(k,j,i) = lad_s(knzb_s_inner(j,i),j,i) 2272 ENDDO 2273 ENDDO 2274 ENDDO 2275 2276 CASE ( 'usm_canopy_khf' ) 2277 ! canopy kinematic heat flux 2278 DO i = nxl, nxr 2279 DO j = nys, nyn 2280 DO k = nzb_s_inner(j,i), nzut 2281 temp_pf(k,j,i) = pc_heating_rate(knzb_s_inner(j,i),j,i) 2282 ENDDO 2283 ENDDO 2284 ENDDO 2285 2286 CASE DEFAULT 2287 found = .FALSE. 2288 2289 END SELECT 2290 2291 ! fill out array local_pf which is subsequently treated by data_output_3d 2292 CALL exchange_horiz( temp_pf, nbgp ) 2293 DO j = nysg,nyng 2294 DO i = nxlg,nxrg 2295 DO k = nzb_do, nzt_do 2296 local_pf(i,j,k) = temp_pf(k,j,i) 2297 ENDDO 2298 ENDDO 2299 ENDDO 2300 2301 END SUBROUTINE usm_data_output_3d 2302 2303 2304 !! 2305 ! 2306 ! Description: 2307 !  2308 !> Soubroutine defines appropriate grid for netcdf variables. 2309 !> It is called out from subroutine netcdf. 2310 !! 2311 SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z ) 2312 2313 IMPLICIT NONE 2314 2315 CHARACTER (len=*), INTENT(IN) :: variable !< 2316 LOGICAL, INTENT(OUT) :: found !< 2317 CHARACTER (len=*), INTENT(OUT) :: grid_x !< 2318 CHARACTER (len=*), INTENT(OUT) :: grid_y !< 2319 CHARACTER (len=*), INTENT(OUT) :: grid_z !< 2320 2321 CHARACTER (len=20) :: var 2322 2323 var = TRIM(variable) 2324 IF ( var(1:11)=='usm_radnet_' .OR. var(1:13) =='usm_rad_insw_' .OR. & 2325 var(1:13) =='usm_rad_inlw_' .OR. var(1:16) =='usm_rad_inswdir_' .OR. & 2326 var(1:16) =='usm_rad_inswdif_' .OR. var(1:16) =='usm_rad_inswref_' .OR. & 2327 var(1:16) =='usm_rad_inlwdif_' .OR. var(1:16) =='usm_rad_inlwref_' .OR. & 2328 var(1:14) =='usm_rad_outsw_' .OR. var(1:14) =='usm_rad_outlw_' .OR. & 2329 var(1:11) =='usm_rad_hf_' .OR. & 2330 var(1:9) == 'usm_wshf_' .OR. var(1:9)== 'usm_wghf_' .OR. & 2331 var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. & 2332 var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. & 2333 var(1:7) =='usm_dif' .OR. var(1:11) =='usm_surfcat' .OR. & 2334 var(1:11) =='usm_surfalb' .OR. var(1:12) =='usm_surfemis' .OR. & 2335 var(1:7) == 'usm_lad' .OR. var(1:14) == 'usm_canopy_khf' ) THEN 2336 2337 found = .TRUE. 2338 grid_x = 'x' 2339 grid_y = 'y' 2340 grid_z = 'zu' 2341 ELSE 2342 found = .FALSE. 2343 grid_x = 'none' 2344 grid_y = 'none' 2345 grid_z = 'none' 2346 ENDIF 2347 2348 END SUBROUTINE usm_define_netcdf_grid 2349 2350 2351 !! 2352 !> Finds first model boundary crossed by a ray 2353 !! 2354 PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross) 2355 IMPLICIT NONE 2356 REAL(wp), DIMENSION(3), INTENT(in) :: origin !< ray origin 2357 REAL(wp), DIMENSION(3), INTENT(in) :: uvect !< ray unit vector 2358 INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross !< found boundary crossing (d, z, y, x) 2359 REAL(wp), DIMENSION(3) :: crossdist !< crossing distance 2360 INTEGER(iwp), DIMENSION(3) :: bdyd !< boundary direction 2361 REAL(wp) :: bdydim !< 2362 REAL(wp) :: dist !< 2363 INTEGER(iwp) :: seldim !< found fist crossing index 2364 INTEGER(iwp) :: d !< 2365 2366 bdydim = nzut + .5_wp !< top boundary 2367 bdyd(1) = isky 2368 crossdist(1) = (bdydim  origin(1)) / uvect(1) 2369 2370 IF ( uvect(2) >= 0._wp ) THEN 2371 bdydim = ny + .5_wp !< north global boundary 2372 bdyd(2) = inorthb 2373 ELSE 2374 bdydim = .5_wp !< south global boundary 2375 bdyd(2) = isouthb 2376 ENDIF 2377 crossdist(2) = (bdydim  origin(2)) / uvect(2) 2378 2379 IF ( uvect(3) >= 0._wp ) THEN 2380 bdydim = nx + .5_wp !< east global boundary 2381 bdyd(3) = ieastb 2382 ELSE 2383 bdydim = .5_wp !< west global boundary 2384 bdyd(3) = iwestb 2385 ENDIF 2386 crossdist(3) = (bdydim  origin(3)) / uvect(3) 2387 2388 seldim = minloc(crossdist, 1) 2389 dist = crossdist(seldim) 2390 d = bdyd(seldim) 2391 2392 bdycross(1) = d 2393 bdycross(2:4) = NINT( origin(:) + uvect(:)*dist & 2394 + .5_wp * (/ kdir(d), jdir(d), idir(d) /) ) 2395 END SUBROUTINE 2396 2397 2398 !! 2399 !> Determines whether two faces are oriented towards each other 2400 !! 2401 PURE LOGICAL FUNCTION usm_facing(x, y, z, d, x2, y2, z2, d2) 2402 IMPLICIT NONE 2403 INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d2 2404 2405 usm_facing = .FALSE. 2406 IF ( d==iroof .AND. d2==iroof ) RETURN 2407 IF ( d==isky .AND. d2==isky ) RETURN 2408 IF ( (d==isouth .OR. d==inorthb) .AND. (d2==isouth.OR.d2==inorthb) ) RETURN 2409 IF ( (d==inorth .OR. d==isouthb) .AND. (d2==inorth.OR.d2==isouthb) ) RETURN 2410 IF ( (d==iwest .OR. d==ieastb) .AND. (d2==iwest.OR.d2==ieastb) ) RETURN 2411 IF ( (d==ieast .OR. d==iwestb) .AND. (d2==ieast.OR.d2==iwestb) ) RETURN 2412 2413 SELECT CASE (d) 2414 CASE (iroof) !< ground, roof 2415 IF ( z2 < z ) RETURN 2416 CASE (isky) !< sky 2417 IF ( z2 > z ) RETURN 2418 CASE (isouth, inorthb) !< south facing 2419 IF ( y2 > y ) RETURN 2420 CASE (inorth, isouthb) !< north facing 2421 IF ( y2 < y ) RETURN 2422 CASE (iwest, ieastb) !< west facing 2423 IF ( x2 > x ) RETURN 2424 CASE (ieast, iwestb) !< east facing 2425 IF ( x2 < x ) RETURN 2426 END SELECT 2427 2428 SELECT CASE (d2) 2429 CASE (iroof) !< ground, roof 2430 IF ( z < z2 ) RETURN 2431 CASE (isky) !< sky 2432 IF ( z > z2 ) RETURN 2433 CASE (isouth, inorthb) !< south facing 2434 IF ( y > y2 ) RETURN 2435 CASE (inorth, isouthb) !< north facing 2436 IF ( y < y2 ) RETURN 2437 CASE (iwest, ieastb) !< west facing 2438 IF ( x > x2 ) RETURN 2439 CASE (ieast, iwestb) !< east facing 2440 IF ( x < x2 ) RETURN 2441 CASE (1) 2442 CONTINUE 2443 END SELECT 2444 2445 usm_facing = .TRUE. 2446 2447 END FUNCTION usm_facing 2448 2449 2450 !! 2451 ! Description: 2452 !  2453 !> Initialization of the wall surface model 2454 !! 2455 SUBROUTINE usm_init_material_model 2456 2457 IMPLICIT NONE 2458 2459 INTEGER(iwp) :: k, l !< running indices 2460 2461 CALL location_message( ' initialization of wall surface model', .TRUE. ) 2462 2463 ! Calculate wall grid spacings. 2464 ! Temperature is defined at the center of the wall layers, 2465 ! whereas gradients/fluxes are defined at the edges (_stag) 2466 DO l = nzb_wall, nzt_wall 2467 zwn(l) = zwn_default(l) 2468 ENDDO 2469 2470 ! apply for all particular wall grids 2471 DO l = startenergy, endenergy 2472 zw(:,l) = zwn(:) * thickness_wall(l) 2473 dz_wall(nzb_wall,l) = zw(nzb_wall,l) 2474 DO k = nzb_wall+1, nzt_wall 2475 dz_wall(k,l) = zw(k,l)  zw(k1,l) 2476 ENDDO 2477 2478 dz_wall(nzt_wall+1,l) = dz_wall(nzt_wall,l) 2479 2480 DO k = nzb_wall, nzt_wall1 2481 dz_wall_stag(k,l) = 0.5 * (dz_wall(k+1,l) + dz_wall(k,l)) 2482 ENDDO 2483 dz_wall_stag(nzt_wall,l) = dz_wall(nzt_wall,l) 2484 ENDDO 2485 2486 ddz_wall = 1.0_wp / dz_wall 2487 ddz_wall_stag = 1.0_wp / dz_wall_stag 2488 2489 CALL location_message( ' wall structures filed out', .TRUE. ) 2490 2491 CALL location_message( ' initialization of wall surface model finished', .TRUE. ) 2492 2493 END SUBROUTINE usm_init_material_model 2494 2495 2496 !! 2497 ! Description: 2498 !  2499 !> Initialization of the urban surface model 2500 !! 2501 SUBROUTINE usm_init_urban_surface 2502 2503 IMPLICIT NONE 2504 2505 INTEGER(iwp) :: i, j, k, l !< running indices 2506 REAL(wp) :: c, d, tin, exn 2507 2508 CALL cpu_log( log_point_s(78), 'usm_init', 'start' ) 2509 ! surface forcing have to be disabled for LSF 2510 ! in case of enabled urban surface module 2511 IF ( large_scale_forcing ) THEN 2512 lsf_surf = .FALSE. 2513 ENDIF 2514 2515 ! init anthropogenic sources of heat 2516 CALL usm_allocate_urban_surface() 2517 2518 ! read the surface_types array somewhere 2519 CALL usm_read_urban_surface_types() 2520 2521 ! init material heat model 2522 CALL usm_init_material_model() 2523 2524 IF ( usm_anthropogenic_heat ) THEN 2525 ! init anthropogenic sources of heat (from transportation for now) 2526 CALL usm_read_anthropogenic_heat() 2527 ENDIF 2528 2529 IF ( read_svf_on_init ) THEN 2530 ! read svf and svfsurf data from file 2531 CALL location_message( ' Start reading SVF from file', .TRUE. ) 2532 CALL usm_read_svf_from_file() 2533 CALL location_message( ' Reading SVF from file has finished', .TRUE. ) 2534 ELSE 2535 ! calculate SFV and CSF 2536 CALL location_message( ' Start calculation of SVF', .TRUE. ) 2537 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'start' ) 2538 CALL usm_calc_svf() 2539 CALL cpu_log( log_point_s(79), 'usm_calc_svf', 'stop' ) 2540 CALL location_message( ' Calculation of SVF has finished', .TRUE. ) 2541 ENDIF 2542 2543 IF ( write_svf_on_init ) THEN 2544 ! write svf and svfsurf data to file 2545 CALL usm_write_svf_to_file() 2546 ENDIF 2547 2548 IF ( plant_canopy ) THEN 2549 ! gridpcbl was only necessary for initialization 2550 DEALLOCATE( gridpcbl ) 2551 IF ( .NOT. ALLOCATED(pc_heating_rate) ) THEN 2552 ! then pc_heating_rate is allocated in init_plant_canopy 2553 ! in case of cthf /= 0 => we need to allocate it for our use here 2554 ALLOCATE( pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2555 ENDIF 2556 ENDIF 2557 2558 ! Intitialization of the surface and wall/ground/roof temperature 2559 2560 ! Initialization for restart runs 2561 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN 2562 2563 ! restore data from restart file 2564 CALL usm_read_restart_data() 2565 ELSE 2566 2567 ! Calculate initial surface temperature 2568 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 2569 2570 DO l = startenergy, endenergy 2571 k = surfl(iz,l) 2572 j = surfl(iy,l) 2573 i = surfl(ix,l) 2574 2575 ! Initial surface temperature set from pt of adjacent gridbox 2576 t_surf(l) = pt(k,j,i) * exn 2577 ENDDO 2578 2579 ! initial values for t_wall 2580 ! outer value is set to surface temperature 2581 ! inner value is set to wall_inner_temperature 2582 ! and profile is logaritmic (linear in nz) 2583 DO l = startenergy, endenergy 2584 IF ( isroof_surf(l) ) THEN 2585 tin = roof_inner_temperature 2586 ELSE IF ( surf(id,l)==iroof ) THEN 2587 tin = soil_inner_temperature 2588 ELSE 2589 tin = wall_inner_temperature 2590 ENDIF 2591 DO k = nzb_wall, nzt_wall+1 2592 c = REAL(knzb_wall,wp)/REAL(nzt_wall+1nzb_wall,wp) 2593 t_wall(k,:) = (1.0_wpc)*t_surf(:) + c*tin 2594 ENDDO 2595 ENDDO 2596 ENDIF 2597 2598 ! 2599 ! Possibly DO userdefined actions (e.g. define heterogeneous wall surface) 2600 CALL user_init_urban_surface 2601 2602 ! initialize prognostic values for the first timestep 2603 t_surf_p = t_surf 2604 t_wall_p = t_wall 2605 2606 ! Adjust radiative fluxes for urban surface at model start 2607 CALL usm_radiation 2608 2609 CALL cpu_log( log_point_s(78), 'usm_init', 'stop' ) 2610 2611 2612 END SUBROUTINE usm_init_urban_surface 2613 2614 2615 !! 2616 ! Description: 2617 !  2618 ! 2619 !> Wall model as part of the urban surface model. The model predicts wall 2620 !> temperature. 2621 !! 2622 SUBROUTINE usm_material_heat_model 2623 2624 2625 IMPLICIT NONE 2626 2627 INTEGER(iwp) :: i,j,k,l,kw !< running indices 2628 2629 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency 2630 2631 2632 DO l = startenergy, endenergy 2633 ! calculate frequently used parameters 2634 k = surfl(iz,l) 2635 j = surfl(iy,l) 2636 i = surfl(ix,l) 2637 2638 ! 2639 ! prognostic equation for ground/wall/roof temperature t_wall 2640 wtend(:) = 0.0_wp 2641 wtend(nzb_wall) = (1.0_wp/rho_c_wall(nzb_wall,l)) * & 2642 ( lambda_h(nzb_wall,l) * ( t_wall(nzb_wall+1,l) & 2643  t_wall(nzb_wall,l) ) * ddz_wall(nzb_wall+1,l) & 2644 + wghf_eb(l) ) * ddz_wall_stag(nzb_wall,l) 2645 2646 DO kw = nzb_wall+1, nzt_wall 2647 wtend(kw) = (1.0_wp/rho_c_wall(kw,l)) & 2648 * ( lambda_h(kw,l) & 2649 * ( t_wall(kw+1,l)  t_wall(kw,l) ) & 2650 * ddz_wall(kw+1,l) & 2651  lambda_h(kw1,l) & 2652 * ( t_wall(kw,l)  t_wall(kw1,l) ) & 2653 * ddz_wall(kw,l) & 2654 ) * ddz_wall_stag(kw,l) 2655 ENDDO 2656 2657 t_wall_p(nzb_wall:nzt_wall,l) = t_wall(nzb_wall:nzt_wall,l) & 2658 + dt_3d * ( tsc(2) & 2659 * wtend(nzb_wall:nzt_wall) + tsc(3) & 2660 * tt_wall_m(nzb_wall:nzt_wall,l) ) 2661 2662 ! 2663 ! calculate t_wall tendencies for the next RungeKutta step 2664 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2665 IF ( intermediate_timestep_count == 1 ) THEN 2666 DO kw = nzb_wall, nzt_wall 2667 tt_wall_m(kw,l) = wtend(kw) 2668 ENDDO 2669 ELSEIF ( intermediate_timestep_count < & 2670 intermediate_timestep_count_max ) THEN 2671 DO kw = nzb_wall, nzt_wall 2672 tt_wall_m(kw,l) = 9.5625_wp * wtend(kw) + 5.3125_wp & 2673 * tt_wall_m(kw,l) 2674 ENDDO 2675 ENDIF 2676 ENDIF 2677 ENDDO 2678 2679 END SUBROUTINE usm_material_heat_model 2680 2681 2682 !! 2683 ! Description: 2684 !  2685 !> This subroutine calculates interaction of the solar radiation 2686 !> with urban surface and updates surface, roofs and walls heatfluxes. 2687 !> It also updates rad_sw_out and rad_lw_out. 2688 !! 2689 SUBROUTINE usm_radiation 2690 2691 IMPLICIT NONE 2692 2693 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep 2694 INTEGER(iwp) :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, ipcgb 2695 INTEGER(iwp), DIMENSION(4) :: bdycross 2696 REAL(wp), DIMENSION(3,3) :: mrot !< grid rotation matrix (xyz) 2697 REAL(wp), DIMENSION(3,0:9) :: vnorm !< face direction normal vectors (xyz) 2698 REAL(wp), DIMENSION(3) :: sunorig !< grid rotated solar direction unit vector (xyz) 2699 REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx) 2700 REAL(wp), DIMENSION(0:9) :: costheta !< direct irradiance factor of solar angle 2701 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temp tendency 2702 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 2703 REAL(wp) :: rx, ry, rz 2704 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 2705 INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy 2706 2707 2708 IF ( plant_canopy ) THEN 2709 pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 2710 / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T) 2711 ENDIF 2712 2713 sun_direction = .TRUE. 2714 CALL calc_zenith !< required also for diffusion radiation 2715 2716 ! prepare rotated normal vectors and irradiance factor 2717 vnorm(1,:) = idir(:) 2718 vnorm(2,:) = jdir(:) 2719 vnorm(3,:) = kdir(:) 2720 mrot(1, :) = (/ cos(alpha), sin(alpha), 0._wp /) 2721 mrot(2, :) = (/ sin(alpha), cos(alpha), 0._wp /) 2722 mrot(3, :) = (/ 0._wp, 0._wp, 1._wp /) 2723 sunorig = (/ sun_dir_lon, sun_dir_lat, zenith(0) /) 2724 sunorig = matmul(mrot, sunorig) 2725 DO d = 0, 9 2726 costheta(d) = dot_product(sunorig, vnorm(:,d)) 2727 ENDDO 2728 2729 IF ( zenith(0) > 0 ) THEN 2730 ! now we will "squash" the sunorig vector by grid box size in 2731 ! each dimension, so that this new direction vector will allow us 2732 ! to traverse the ray path within grid coordinates directly 2733 sunorig_grid = (/ sunorig(3)/dz, sunorig(2)/dy, sunorig(1)/dx /) 2734 ! sunorig_grid = sunorig_grid / norm2(sunorig_grid) 2735 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) 2736 2737 IF ( plant_canopy ) THEN 2738 ! precompute effective box depth with prototype Leaf Area Density 2739 pc_box_dimshift = maxloc(sunorig, 1)  1 2740 CALL usm_box_absorb(cshift((/dx,dy,dz/), pc_box_dimshift), & 2741 60, prototype_lad, & 2742 cshift(sunorig, pc_box_dimshift), & 2743 pc_box_area, pc_abs_frac) 2744 pc_box_area = pc_box_area * sunorig(pc_box_dimshift+1) / sunorig(3) 2745 pc_abs_eff = log(1._wp  pc_abs_frac) / prototype_lad 2746 ENDIF 2747 ENDIF 2748 2749 ! split diffusion and direct part of the solar downward radiation 2750 ! comming from radiation model and store it in 2D arrays 2751 ! rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff 2752 IF ( split_diffusion_radiation ) THEN 2753 CALL usm_calc_diffusion_radiation 2754 ELSE 2755 rad_sw_in_diff = 0.0_wp 2756 rad_sw_in_dir(:,:) = rad_sw_in(0,:,:) 2757 rad_lw_in_diff(:,:) = rad_lw_in(0,:,:) 2758 ENDIF 2759 2760 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2761 ! First pass: direct + diffuse irradiance 2762 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2763 surfinswdir = 0._wp 2764 surfinswdif = 0._wp 2765 surfinlwdif = 0._wp 2766 surfins = 0._wp 2767 surfinl = 0._wp 2768 surfoutsl = 0._wp 2769 surfoutll = 0._wp 2770 2771 ! Set up thermal radiation from surfaces 2772 ! emiss_surf is defined only for surfaces for which energy balance is calculated 2773 surfoutll(startenergy:endenergy) = emiss_surf(startenergy:endenergy) * sigma_sb & 2774 * t_surf(startenergy:endenergy)**4 2775 2776 #if defined( __parallel ) 2777 ! might be optimized and gather only values relevant for current processor 2778 CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, & 2779 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2780 #else 2781 surfoutl(:) = surfoutll(:) 2782 #endif 2783 2784 isurf1 = 1 !< previous processed surface 2785 DO isvf = 1, nsvfl 2786 isurf = svfsurf(1, isvf) 2787 k = surfl(iz, isurf) 2788 j = surfl(iy, isurf) 2789 i = surfl(ix, isurf) 2790 isurfsrc = svfsurf(2, isvf) 2791 IF ( zenith(0) > 0 .AND. isurf /= isurf1 ) THEN 2792 ! locate the virtual surface where the direct solar ray crosses domain boundary 2793 ! (once per target surface) 2794 d = surfl(id, isurf) 2795 rz = REAL(k, wp)  0.5_wp * kdir(d) 2796 ry = REAL(j, wp)  0.5_wp * jdir(d) 2797 rx = REAL(i, wp)  0.5_wp * idir(d) 2798 2799 CALL usm_find_boundary_face( (/ rz, ry, rx /), sunorig_grid, bdycross) 2800 2801 isurf1 = isurf 2802 ENDIF 2803 2804 IF ( surf(id, isurfsrc) >= isky ) THEN 2805 ! diffuse rad from boundary surfaces. Since it is a simply 2806 ! calculated value, it is not assigned to surfref(s/l), 2807 ! instead it is used directly here 2808 ! we consider the radiation from the radiation model falling on surface 2809 ! as the radiation falling on the top of urban layer into the place of the source surface 2810 ! we consider it as a very reasonable simplification which allow as avoid 2811 ! necessity of other global range arrays and some all to all mpi communication 2812 surfinswdif(isurf) = surfinswdif(isurf) + rad_sw_in_diff(j,i) * svf(1,isvf) * svf(2,isvf) 2813 !< canopy shading is applied only to shortwave 2814 surfinlwdif(isurf) = surfinlwdif(isurf) + rad_lw_in_diff(j,i) * svf(1,isvf) 2815 ELSE 2816 ! for surfacetosurface factors we calculate thermal radiation in 1st pass 2817 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 2818 ENDIF 2819 2820 IF ( zenith(0) > 0 .AND. all( surf(:, isurfsrc) == bdycross ) ) THEN 2821 ! found svf between model boundary and the face => face isn't shaded 2822 surfinswdir(isurf) = rad_sw_in_dir(j, i) & 2823 * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0) 2824 2825 ENDIF 2826 ENDDO 2827 2828 IF ( plant_canopy ) THEN 2829 2830 pcbinsw(:) = 0._wp 2831 pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore 2832 ! 2833 ! pcsf first pass 2834 isurf1 = 1 !< previous processed pcgb 2835 DO isvf = nsvfl+1, nsvfcsfl 2836 ipcgb = svfsurf(1, isvf) 2837 i = pcbl(ix,ipcgb) 2838 j = pcbl(iy,ipcgb) 2839 k = pcbl(iz,ipcgb) 2840 isurfsrc = svfsurf(2, isvf) 2841 2842 IF ( zenith(0) > 0 .AND. ipcgb /= isurf1 ) THEN 2843 ! locate the virtual surface where the direct solar ray crosses domain boundary 2844 ! (once per target PC gridbox) 2845 rz = REAL(k, wp) 2846 ry = REAL(j, wp) 2847 rx = REAL(i, wp) 2848 CALL usm_find_boundary_face( (/ rz, ry, rx /), & 2849 sunorig_grid, bdycross) 2850 2851 isurf1 = ipcgb 2852 ENDIF 2853 2854 IF ( surf(id, isurfsrc) >= isky ) THEN 2855 ! Diffuse rad from boundary surfaces. See comments for svf above. 2856 pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * rad_sw_in_diff(j,i) 2857 ! canopy shading is applied only to shortwave, therefore no absorbtion for lw 2858 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * rad_lw_in_diff(j,i) 2859 !ELSE 2860 ! Thermal radiation in 1st pass 2861 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc) 2862 ENDIF 2863 2864 IF ( zenith(0) > 0 .AND. all( surf(:, isurfsrc) == bdycross ) ) THEN 2865 ! found svf between model boundary and the pcgb => pcgb isn't shaded 2866 pc_abs_frac = 1._wp  exp(pc_abs_eff * lad_s(k,j,i)) 2867 pcbinsw(ipcgb) = pcbinsw(ipcgb) & 2868 + rad_sw_in_dir(j, i) * pc_box_area * svf(2,isvf) * pc_abs_frac 2869 ENDIF 2870 ENDDO 2871 ENDIF 2872 surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy) 2873 surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy) 2874 surfinsw(:) = surfins(:) 2875 surfinlw(:) = surfinl(:) 2876 surfoutsw(:) = 0.0_wp 2877 surfoutlw(:) = surfoutll(:) 2878 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 2879  surfoutsw(startenergy:endenergy)  surfoutlw(startenergy:endenergy) 2880 2881 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2882 ! Next passes  reflections 2883 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2884 DO refstep = 1, nrefsteps 2885 2886 surfoutsl(startenergy:endenergy) = albedo_surf(startenergy:endenergy) * surfins(startenergy:endenergy) 2887 ! for nontransparent surfaces, longwave albedo is 1  emissivity 2888 surfoutll(startenergy:endenergy) = (1._wp  emiss_surf(startenergy:endenergy)) * surfinl(startenergy:endenergy) 2889 2890 #if defined( __parallel ) 2891 CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, & 2892 surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2893 CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, & 2894 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2895 #else 2896 surfouts(:) = surfoutsl(:) 2897 surfoutl(:) = surfoutll(:) 2898 #endif 2899 2900 ! reset for next pass input 2901 surfins(:) = 0._wp 2902 surfinl(:) = 0._wp 2903 2904 ! reflected radiation 2905 DO isvf = 1, nsvfl 2906 isurf = svfsurf(1, isvf) 2907 isurfsrc = svfsurf(2, isvf) 2908 2909 ! TODO: to remove if, use start+end for isvf 2910 IF ( surf(id, isurfsrc) < isky ) THEN 2911 surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 2912 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 2913 ENDIF 2914 ENDDO 2915 2916 ! radiation absorbed by plant canopy 2917 DO isvf = nsvfl+1, nsvfcsfl 2918 ipcgb = svfsurf(1, isvf) 2919 isurfsrc = svfsurf(2, isvf) 2920 2921 IF ( surf(id, isurfsrc) < isky ) THEN 2922 pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 2923 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc) 2924 ENDIF 2925 ENDDO 2926 2927 surfinsw(:) = surfinsw(:) + surfins(:) 2928 surfinlw(:) = surfinlw(:) + surfinl(:) 2929 surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy) 2930 surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy) 2931 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 2932  surfoutsw(startenergy:endenergy)  surfoutlw(startenergy:endenergy) 2933 2934 ENDDO 2935 2936 ! push heat flux absorbed by plant canopy to respective 3D arrays 2937 IF ( plant_canopy ) THEN 2938 pc_heating_rate(:,:,:) = 0._wp 2939 DO ipcgb = 1, npcbl 2940 j = pcbl(iy, ipcgb) 2941 i = pcbl(ix, ipcgb) 2942 k = pcbl(iz, ipcgb) 2943 kk = k  nzb_s_inner(j,i) ! lad arrays are defined flat 2944 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 2945 * pchf_prep(k) * pt(k, j, i) ! = dT/dt 2946 ENDDO 2947 ENDIF 2948 2949 ! return surface radiation to horizontal surfaces 2950 ! to rad_sw_in, rad_lw_in and rad_net for outputs 2951 !!!!!!!!!! 2952 ! we need the original radiation on urban top layer 2953 ! for calculation of MRT so we can't do adjustment here for now 2954 !!!!!!!!!! 2955 !!!DO isurf = 1, nsurfl 2956 !!! i = surfl(ix,isurf) 2957 !!! j = surfl(iy,isurf) 2958 !!! k = surfl(iz,isurf) 2959 !!! d = surfl(id,isurf) 2960 !!! IF ( d==iroof ) THEN 2961 !!! rad_sw_in(:,j,i) = surfinsw(isurf) 2962 !!! rad_lw_in(:,j,i) = surfinlw(isurf) 2963 !!! rad_net(j,i) = rad_sw_in(k,j,i)  rad_sw_out(k,j,i) + rad_lw_in(k,j,i)  rad_lw_out(k,j,i) 2964 !!! ENDIF 2965 !!!ENDDO 2966 2967 END SUBROUTINE usm_radiation 2968 1576 2969 1577 2970 !! … … 1626 3019 REAL(wp), PARAMETER :: grow_factor = 1.5_wp 1627 3020 1628 !Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also1629 !the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.3021 ! Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also 3022 ! the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor. 1630 3023 maxboxes = SUM(ABS(NINT(targ)  NINT(src))) + 1 1631 !IF ( debug_prints ) THEN1632 ! WRITE(9, *) 'Raytracing from ', src, ' to ', targ, ' with max ', maxboxes, ' boxes, ncsfl = ', ncsfl, ', ncsfla = ', ncsfla1633 ! FLUSH(9)1634 !ENDIF1635 3024 IF ( plant_canopy .AND. ncsfl + maxboxes > ncsfla ) THEN 1636 !use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)1637 !k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &1638 !/ log(grow_factor)), kind=wp))1639 !or use this code to simply always keep some extra space after growing3025 ! use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) 3026 ! k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) & 3027 ! / log(grow_factor)), kind=wp)) 3028 ! or use this code to simply always keep some extra space after growing 1640 3029 k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor) 1641 3030 1642 CALL merge_and_grow_csf(k)3031 CALL usm_merge_and_grow_csf(k) 1643 3032 ENDIF 1644 3033 … … 1648 3037 delta(:) = targ(:)  src(:) 1649 3038 distance = SQRT(SUM(delta(:)**2)) 1650 IF ( distance == 0._wp)THEN3039 IF ( distance == 0._wp ) THEN 1651 3040 visible = .TRUE. 1652 3041 RETURN … … 1657 3046 lastdist = 0._wp 1658 3047 1659 !Since all face coordinates have values *.5 and we'd like to use1660 !integers, all these have .5 added3048 ! Since all face coordinates have values *.5 and we'd like to use 3049 ! integers, all these have .5 added 1661 3050 DO d = 1, 3 1662 IF ( uvect(d) == 0._wp ) THEN3051 IF ( uvect(d) == 0._wp ) THEN 1663 3052 dimnext(d) = 999999999 1664 3053 dimdelta(d) = 999999999 1665 3054 dimnextdist(d) = 1.0E20_wp 1666 ELSE IF ( uvect(d) > 0._wp ) THEN3055 ELSE IF ( uvect(d) > 0._wp ) THEN 1667 3056 dimnext(d) = CEILING(src(d) + .5_wp) 1668 3057 dimdelta(d) = 1 … … 1675 3064 ENDDO 1676 3065 1677 ! dimnextdist(:) = (dimnext(:)  .5_wp  src(:)) / uvect(:)1678 ! will assign Infinity to those dimensions that have1679 ! uvect==0, which is supported by minloc1680 !kanani: with ifort compiler option fpe0 we get error "division by zero"; there must be a compilerindependent way to handle the setting of dimnextdist to infinity in case of uvect==0. > This is accounted for in the above DO loop1681 1682 3066 DO 1683 !along what dimension will the next wall crossing be?3067 ! along what dimension will the next wall crossing be? 1684 3068 seldim = minloc(dimnextdist, 1) 1685 3069 nextdist = dimnextdist(seldim) … … 1687 3071 1688 3072 crlen = nextdist  lastdist 1689 IF ( crlen > .001_wp)THEN3073 IF ( crlen > .001_wp ) THEN 1690 3074 crmid = (lastdist + nextdist) * .5_wp 1691 3075 box = NINT(src(:) + uvect(:) * crmid) 1692 3076 1693 !calculate index of the grid with global indices (box(2),box(3))1694 !in the array nzterr and plantt and id of the coresponding processor3077 ! calculate index of the grid with global indices (box(2),box(3)) 3078 ! in the array nzterr and plantt and id of the coresponding processor 1695 3079 px = box(3)/nnx 1696 3080 py = box(2)/nny 1697 3081 ip = px*pdims(2)+py 1698 3082 ig = ip*nnx*nny + (box(3)px*nnx)*nny + box(2)py*nny 1699 IF ( box(1) <= nzterr(ig) ) THEN3083 IF ( box(1) <= nzterr(ig) ) THEN 1700 3084 visible = .FALSE. 1701 IF ( ncsb > 0 ) THEN1702 !rewind written plant canopy sink factors  they are invalid3085 IF ( ncsb > 0 ) THEN 3086 ! rewind written plant canopy sink factors  they are invalid 1703 3087 ncsfl = ncsfl  ncsb 1704 3088 ENDIF … … 1706 3090 ENDIF 1707 3091 1708 IF ( plant_canopy)THEN1709 IF ( box(1) <= plantt(ig) ) THEN3092 IF ( plant_canopy ) THEN 3093 IF ( box(1) <= plantt(ig) ) THEN 1710 3094 #if defined( __parallel ) 1711 3095 lad_disp = (box(3)px*nnx)*(nny*nzu) + (box(2)py*nny)*nzu + box(1)nzub 1712 3096 IF ( usm_lad_rma ) THEN 1713 !Read LAD using MPI RMA3097 ! Read LAD using MPI RMA 1714 3098 CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' ) 1715 3099 CALL MPI_Get(lad_s_target, 1, MPI_REAL, ip, lad_disp, 1, MPI_REAL, & 1716 3100 win_lad, ierr) 1717 IF ( ierr /= 0 ) THEN3101 IF ( ierr /= 0 ) THEN 1718 3102 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get' 1719 3103 CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 ) … … 1721 3105 1722 3106 CALL MPI_Win_flush_local(ip, win_lad, ierr) 1723 IF ( ierr /= 0 ) THEN3107 IF ( ierr /= 0 ) THEN 1724 3108 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local' 1725 3109 CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 ) … … 1735 3119 * crlen*realdist) 1736 3120 1737 IF ( create_csf ) THEN1738 !write svf values into the array3121 IF ( create_csf ) THEN 3122 ! write svf values into the array 1739 3123 ncsb = ncsb + 1 1740 3124 ncsfl = ncsfl + 1 … … 1744 3128 acsf(ncsfl)%itz = box(1) 1745 3129 acsf(ncsfl)%isurfs = isrc 1746 acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !we postpone multiplication by transparency3130 acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) ! we postpone multiplication by transparency 1747 3131 acsf(ncsfl)%rtransp = REAL(transparency, wp) 1748 3132 ENDIF !< create_csf … … 1763 3147 END SUBROUTINE usm_raytrace 1764 3148 1765 1766 1767 !! 1768 !> Calculates radiation absorbed by box with given size and LAD. 1769 !> 1770 !> Simulates resol**2 rays (by equally spacing a bounding horizontal square 1771 !> conatining all possible rays that would cross the box) and calculates 1772 !> average transparency per ray. Returns fraction of absorbed radiation flux 1773 !> and area for which this fraction is effective. 1774 !! 1775 PURE SUBROUTINE usm_box_absorb(boxsize, resol, dens, uvec, area, absorb) 1776 IMPLICIT NONE 1777 1778 REAL(wp), DIMENSION(3), INTENT(in) :: & 1779 boxsize, & !< z, y, x size of box in m 1780 uvec !< z, y, x unit vector of incoming flux 1781 INTEGER(iwp), INTENT(in) :: & 1782 resol !< No. of rays in x and y dimensions 1783 REAL(wp), INTENT(in) :: & 1784 dens !< box density (e.g. Leaf Area Density) 1785 REAL(wp), INTENT(out) :: & 1786 area, & !< horizontal area for flux absorbtion 1787 absorb !< fraction of absorbed flux 1788 REAL(wp) :: & 1789 xshift, yshift, & 1790 xmin, xmax, ymin, ymax, & 1791 xorig, yorig, & 1792 dx1, dy1, dz1, dx2, dy2, dz2, & 1793 crdist, & 1794 transp 1795 INTEGER(iwp) :: & 1796 i, j 1797 1798 xshift = uvec(3) / uvec(1) * boxsize(1) 1799 xmin = min(0._wp, xshift) 1800 xmax = boxsize(3) + max(0._wp, xshift) 1801 yshift = uvec(2) / uvec(1) * boxsize(1) 1802 ymin = min(0._wp, yshift) 1803 ymax = boxsize(2) + max(0._wp, yshift) 1804 1805 transp = 0._wp 1806 DO i = 1, resol 1807 xorig = xmin + (xmaxxmin) * (i.5_wp) / resol 1808 DO j = 1, resol 1809 yorig = ymin + (ymaxymin) * (j.5_wp) / resol 1810 1811 dz1 = 0._wp 1812 dz2 = boxsize(1)/uvec(1) 1813 1814 IF ( uvec(2) > 0._wp ) THEN 1815 dy1 = yorig / uvec(2) !< crossing with y=0 1816 dy2 = (boxsize(2)yorig) / uvec(2) !< crossing with y=boxsize(2) 1817 ELSE IF (uvec(2) < 0._wp ) THEN 1818 dy1 = (boxsize(2)yorig) / uvec(2) !< crossing with y=boxsize(2) 1819 dy2 = yorig / uvec(2) !< crossing with y=0 1820 ELSE !uvec(2)==0 1821 dy1 = huge(1._wp) 1822 dy2 = huge(1._wp) 1823 ENDIF 1824 1825 IF ( uvec(3) > 0._wp ) THEN 1826 dx1 = xorig / uvec(3) !< crossing with x=0 1827 dx2 = (boxsize(3)xorig) / uvec(3) !< crossing with x=boxsize(3) 1828 ELSE IF (uvec(3) < 0._wp ) THEN 1829 dx1 = (boxsize(3)xorig) / uvec(3) !< crossing with x=boxsize(3) 1830 dx2 = xorig / uvec(3) !< crossing with x=0 1831 ELSE !uvec(1)==0 1832 dx1 = huge(1._wp) 1833 dx2 = huge(1._wp) 1834 ENDIF 1835 1836 crdist = max(0._wp, (min(dz2, dy2, dx2)  max(dz1, dy1, dx1))) 1837 transp = transp + exp(ext_coef * dens * crdist) 1838 ENDDO 1839 ENDDO 1840 transp = transp / resol**2 1841 area = (boxsize(3)+xshift)*(boxsize(2)+yshift) 1842 absorb = 1._wp  transp 1843 1844 END SUBROUTINE usm_box_absorb 1845 1846 1847 !! 1848 !> Finds first model boundary crossed by a ray 1849 !! 1850 PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross) 1851 IMPLICIT NONE 1852 REAL(wp), DIMENSION(3), INTENT(in) :: origin !< ray origin 1853 REAL(wp), DIMENSION(3), INTENT(in) :: uvect !< ray unit vector 1854 INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross !< found boundary crossing (d, z, y, x) 1855 REAL(wp), DIMENSION(3) :: crossdist !< crossing distance 1856 INTEGER(iwp), DIMENSION(3) :: bdyd !< boundary direction 1857 REAL(wp) :: bdydim !< 1858 REAL(wp) :: dist !< 1859 INTEGER(iwp) :: seldim !< found fist crossing index 1860 INTEGER(iwp) :: d !< 1861 1862 bdydim = nzut + .5_wp !< top boundary 1863 bdyd(1) = isky 1864 crossdist(1) = (bdydim  origin(1)) / uvect(1) 1865 1866 IF ( uvect(2) >= 0._wp ) THEN 1867 bdydim = ny + .5_wp !< north global boundary 1868 bdyd(2) = inorthb 1869 ELSE 1870 bdydim = .5_wp !< south global boundary 1871 bdyd(2) = isouthb 1872 ENDIF 1873 crossdist(2) = (bdydim  origin(2)) / uvect(2) 1874 1875 IF ( uvect(3) >= 0._wp ) THEN 1876 bdydim = nx + .5_wp !< east global boundary 1877 bdyd(3) = ieastb 1878 ELSE 1879 bdydim = .5_wp !< west global boundary 1880 bdyd(3) = iwestb 1881 ENDIF 1882 crossdist(3) = (bdydim  origin(3)) / uvect(3) 1883 1884 seldim = minloc(crossdist, 1) 1885 dist = crossdist(seldim) 1886 d = bdyd(seldim) 1887 1888 bdycross(1) = d 1889 bdycross(2:4) = NINT( origin(:) + uvect(:)*dist & 1890 + .5_wp * (/ kdir(d), jdir(d), idir(d) /) ) 1891 END SUBROUTINE 1892 1893 1894 !! 1895 !> Determines whether two faces are oriented towards each other 1896 !! 1897 PURE LOGICAL FUNCTION usm_facing(x, y, z, d, x2, y2, z2, d2) 1898 IMPLICIT NONE 1899 INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d2 1900 1901 usm_facing = .FALSE. 1902 IF ( d==iroof .AND. d2==iroof ) RETURN 1903 IF ( d==isky .AND. d2==isky ) RETURN 1904 IF ( (d==isouth .OR. d==inorthb) .AND. (d2==isouth.OR.d2==inorthb) ) RETURN 1905 IF ( (d==inorth .OR. d==isouthb) .AND. (d2==inorth.OR.d2==isouthb) ) RETURN 1906 IF ( (d==iwest .OR. d==ieastb) .AND. (d2==iwest.OR.d2==ieastb) ) RETURN 1907 IF ( (d==ieast .OR. d==iwestb) .AND. (d2==ieast.OR.d2==iwestb) ) RETURN 1908 1909 SELECT CASE (d) 1910 CASE (iroof) !< ground, roof 1911 IF ( z2 < z ) RETURN 1912 CASE (isky) !< sky 1913 IF ( z2 > z ) RETURN 1914 CASE (isouth, inorthb) !< south facing 1915 IF ( y2 > y ) RETURN 1916 CASE (inorth, isouthb) !< north facing 1917 IF ( y2 < y ) RETURN 1918 CASE (iwest, ieastb) !< west facing 1919 IF ( x2 > x ) RETURN 1920 CASE (ieast, iwestb) !< east facing 1921 IF ( x2 < x ) RETURN 1922 END SELECT 1923 1924 SELECT CASE (d2) 1925 CASE (iroof) !< ground, roof 1926 IF ( z < z2 ) RETURN 1927 CASE (isky) !< sky 1928 IF ( z > z2 ) RETURN 1929 CASE (isouth, inorthb) !< south facing 1930 IF ( y > y2 ) RETURN 1931 CASE (inorth, isouthb) !< north facing 1932 IF ( y < y2 ) RETURN 1933 CASE (iwest, ieastb) !< west facing 1934 IF ( x > x2 ) RETURN 1935 CASE (ieast, iwestb) !< east facing 1936 IF ( x < x2 ) RETURN 1937 CASE (1) 1938 CONTINUE 1939 END SELECT 1940 1941 usm_facing = .TRUE. 1942 1943 END FUNCTION usm_facing 1944 1945 1946 !! 1947 ! Description: 1948 !  1949 !> This subroutine calculates interaction of the solar radiation 1950 !> with urban surface and updates surface, roofs and walls heatfluxes. 1951 !> It also updates rad_sw_out and rad_lw_out. 1952 !! 1953 SUBROUTINE usm_radiation 1954 1955 IMPLICIT NONE 1956 1957 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep 1958 INTEGER(iwp) :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, ipcgb 1959 INTEGER(iwp), DIMENSION(4) :: bdycross 1960 REAL(wp), DIMENSION(3,3) :: mrot !< grid rotation matrix (xyz) 1961 REAL(wp), DIMENSION(3,0:9) :: vnorm !< face direction normal vectors (xyz) 1962 REAL(wp), DIMENSION(3) :: sunorig !< grid rotated solar direction unit vector (xyz) 1963 REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx) 1964 REAL(wp), DIMENSION(0:9) :: costheta !< direct irradiance factor of solar angle 1965 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temp tendency 1966 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 1967 REAL(wp) :: rx, ry, rz 1968 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 1969 INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy 1970 1971 1972 IF (plant_canopy) THEN 1973 pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 1974 / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T) 1975 ENDIF 1976 1977 sun_direction = .TRUE. 1978 CALL calc_zenith !< required also for diffusion radiation 1979 1980 ! prepare rotated normal vectors and irradiance factor 1981 vnorm(1,:) = idir(:) 1982 vnorm(2,:) = jdir(:) 1983 vnorm(3,:) = kdir(:) 1984 mrot(1, :) = (/ cos(alpha), sin(alpha), 0._wp /) 1985 mrot(2, :) = (/ sin(alpha), cos(alpha), 0._wp /) 1986 mrot(3, :) = (/ 0._wp, 0._wp, 1._wp /) 1987 sunorig = (/ sun_dir_lon, sun_dir_lat, zenith(0) /) 1988 sunorig = matmul(mrot, sunorig) 1989 DO d = 0, 9 1990 costheta(d) = dot_product(sunorig, vnorm(:,d)) 1991 ENDDO 1992 1993 IF (zenith(0) > 0) THEN 1994 ! now we will "squash" the sunorig vector by grid box size in 1995 ! each dimension, so that this new direction vector will allow us 1996 ! to traverse the ray path within grid coordinates directly 1997 sunorig_grid = (/ sunorig(3)/dz, sunorig(2)/dy, sunorig(1)/dx /) 1998 ! sunorig_grid = sunorig_grid / norm2(sunorig_grid) 1999 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) 2000 2001 IF ( plant_canopy ) THEN 2002 ! precompute effective box depth with prototype Leaf Area Density 2003 pc_box_dimshift = maxloc(sunorig, 1)  1 2004 CALL usm_box_absorb(cshift((/dx,dy,dz/), pc_box_dimshift), & 2005 60, prototype_lad, & 2006 cshift(sunorig, pc_box_dimshift), & 2007 pc_box_area, pc_abs_frac) 2008 pc_box_area = pc_box_area * sunorig(pc_box_dimshift+1) / sunorig(3) 2009 pc_abs_eff = log(1._wp  pc_abs_frac) / prototype_lad 2010 ENDIF 2011 ENDIF 2012 2013 ! split diffusion and direct part of the solar downward radiation 2014 ! comming from radiation model and store it in 2D arrays 2015 ! rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff 2016 IF (split_diffusion_radiation) THEN 2017 CALL usm_calc_diffusion_radiation 2018 ELSE 2019 rad_sw_in_diff = 0.0_wp 2020 rad_sw_in_dir(:,:) = rad_sw_in(0,:,:) 2021 rad_lw_in_diff(:,:) = rad_lw_in(0,:,:) 2022 ENDIF 2023 2024 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2025 ! First pass: direct + diffuse irradiance 2026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2027 surfinswdir = 0._wp 2028 surfinswdif = 0._wp 2029 surfinlwdif = 0._wp 2030 surfins = 0._wp 2031 surfinl = 0._wp 2032 surfoutsl = 0._wp 2033 surfoutll = 0._wp 2034 2035 ! Set up thermal radiation from surfaces 2036 ! emiss_surf is defined only for surfaces for which energy balance is calculated 2037 surfoutll(startenergy:endenergy) = emiss_surf(startenergy:endenergy) * sigma_sb & 2038 * t_surf(startenergy:endenergy)**4 2039 2040 #if defined( __parallel ) 2041 ! might be optimized and gather only values relevant for current processor 2042 CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, & 2043 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2044 #else 2045 surfoutl(:) = surfoutll(:) 2046 #endif 2047 2048 isurf1 = 1 !< previous processed surface 2049 DO isvf = 1, nsvfl 2050 isurf = svfsurf(1, isvf) 2051 k = surfl(iz, isurf) 2052 j = surfl(iy, isurf) 2053 i = surfl(ix, isurf) 2054 isurfsrc = svfsurf(2, isvf) 2055 IF ( zenith(0) > 0 .AND. isurf /= isurf1 ) THEN 2056 ! locate the virtual surface where the direct solar ray crosses domain boundary 2057 ! (once per target surface) 2058 d = surfl(id, isurf) 2059 rz = REAL(k, wp)  0.5_wp * kdir(d) 2060 ry = REAL(j, wp)  0.5_wp * jdir(d) 2061 rx = REAL(i, wp)  0.5_wp * idir(d) 2062 2063 CALL usm_find_boundary_face( (/ rz, ry, rx /), sunorig_grid, bdycross) 2064 2065 isurf1 = isurf 2066 ENDIF 2067 2068 IF ( surf(id, isurfsrc) >= isky ) THEN 2069 ! diffuse rad from boundary surfaces. Since it is a simply 2070 ! calculated value, it is not assigned to surfref(s/l), 2071 ! instead it is used directly here 2072 ! we consider the radiation from the radiation model falling on surface 2073 ! as the radiation falling on the top of urban layer into the place of the source surface 2074 ! we consider it as a very reasonable simplification which allow as avoid 2075 ! necessity of other global range arrays and some all to all mpi communication 2076 surfinswdif(isurf) = surfinswdif(isurf) + rad_sw_in_diff(j,i) * svf(1,isvf) * svf(2,isvf) 2077 !< canopy shading is applied only to shortwave 2078 surfinlwdif(isurf) = surfinlwdif(isurf) + rad_lw_in_diff(j,i) * svf(1,isvf) 2079 ELSE 2080 ! for surfacetosurface factors we calculate thermal radiation in 1st pass 2081 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 2082 ENDIF 2083 2084 IF ( zenith(0) > 0 .AND. all( surf(:, isurfsrc) == bdycross ) ) THEN 2085 ! found svf between model boundary and the face => face isn't shaded 2086 surfinswdir(isurf) = rad_sw_in_dir(j, i) & 2087 * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0) 2088 2089 ENDIF 2090 ENDDO 2091 2092 IF ( plant_canopy ) THEN 2093 2094 pcbinsw(:) = 0._wp 2095 pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore 2096 ! 2097 ! pcsf first pass 2098 isurf1 = 1 !< previous processed pcgb 2099 DO isvf = nsvfl+1, nsvfcsfl 2100 ipcgb = svfsurf(1, isvf) 2101 i = pcbl(ix,ipcgb) 2102 j = pcbl(iy,ipcgb) 2103 k = pcbl(iz,ipcgb) 2104 isurfsrc = svfsurf(2, isvf) 2105 2106 IF ( zenith(0) > 0 .AND. ipcgb /= isurf1 ) THEN 2107 ! locate the virtual surface where the direct solar ray crosses domain boundary 2108 ! (once per target PC gridbox) 2109 rz = REAL(k, wp) 2110 ry = REAL(j, wp) 2111 rx = REAL(i, wp) 2112 CALL usm_find_boundary_face( (/ rz, ry, rx /), & 2113 sunorig_grid, bdycross) 2114 2115 isurf1 = ipcgb 2116 ENDIF 2117 2118 IF ( surf(id, isurfsrc) >= isky ) THEN 2119 ! Diffuse rad from boundary surfaces. See comments for svf above. 2120 pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * rad_sw_in_diff(j,i) 2121 ! canopy shading is applied only to shortwave, therefore no absorbtion for lw 2122 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * rad_lw_in_diff(j,i) 2123 !ELSE 2124 ! Thermal radiation in 1st pass 2125 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc) 2126 ENDIF 2127 2128 IF ( zenith(0) > 0 .AND. all( surf(:, isurfsrc) == bdycross ) ) THEN 2129 ! found svf between model boundary and the pcgb => pcgb isn't shaded 2130 pc_abs_frac = 1._wp  exp(pc_abs_eff * lad_s(k,j,i)) 2131 pcbinsw(ipcgb) = pcbinsw(ipcgb) & 2132 + rad_sw_in_dir(j, i) * pc_box_area * svf(2,isvf) * pc_abs_frac 2133 ENDIF 2134 ENDDO 2135 ENDIF 2136 surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy) 2137 surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy) 2138 surfinsw(:) = surfins(:) 2139 surfinlw(:) = surfinl(:) 2140 surfoutsw(:) = 0.0_wp 2141 surfoutlw(:) = surfoutll(:) 2142 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 2143  surfoutsw(startenergy:endenergy)  surfoutlw(startenergy:endenergy) 2144 2145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2146 ! Next passes  reflections 2147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2148 DO refstep = 1, nrefsteps 2149 2150 surfoutsl(startenergy:endenergy) = albedo_surf(startenergy:endenergy) * surfins(startenergy:endenergy) 2151 ! for nontransparent surfaces, longwave albedo is 1  emissivity 2152 surfoutll(startenergy:endenergy) = (1._wp  emiss_surf(startenergy:endenergy)) * surfinl(startenergy:endenergy) 2153 2154 #if defined( __parallel ) 2155 CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, & 2156 surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2157 CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, & 2158 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 2159 #else 2160 surfouts(:) = surfoutsl(:) 2161 surfoutl(:) = surfoutll(:) 2162 #endif 2163 2164 ! reset for next pass input 2165 surfins(:) = 0._wp 2166 surfinl(:) = 0._wp 2167 2168 ! reflected radiation 2169 DO isvf = 1, nsvfl 2170 isurf = svfsurf(1, isvf) 2171 isurfsrc = svfsurf(2, isvf) 2172 2173 ! TODO: to remove if, use start+end for isvf 2174 IF ( surf(id, isurfsrc) < isky ) THEN 2175 surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 2176 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 2177 ENDIF 2178 ENDDO 2179 2180 ! radiation absorbed by plant canopy 2181 DO isvf = nsvfl+1, nsvfcsfl 2182 ipcgb = svfsurf(1, isvf) 2183 isurfsrc = svfsurf(2, isvf) 2184 2185 IF ( surf(id, isurfsrc) < isky ) THEN 2186 pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 2187 ! pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc) 2188 ENDIF 2189 ENDDO 2190 2191 surfinsw(:) = surfinsw(:) + surfins(:) 2192 surfinlw(:) = surfinlw(:) + surfinl(:) 2193 surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy) 2194 surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy) 2195 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 2196  surfoutsw(startenergy:endenergy)  surfoutlw(startenergy:endenergy) 2197 2198 ENDDO 2199 2200 ! push heat flux absorbed by plant canopy to respective 3D arrays 2201 IF (plant_canopy) THEN 2202 canopy_heat_flux(:,:,:) = 0._wp 2203 DO ipcgb = 1, npcbl 2204 j = pcbl(iy, ipcgb) 2205 i = pcbl(ix, ipcgb) 2206 k = pcbl(iz, ipcgb) 2207 kk = k  nzb_s_inner(j,i) ! lad arrays are defined flat 2208 canopy_heat_flux(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 2209 * pchf_prep(k) * pt(k, j, i) ! = dT/dt 2210 ENDDO 2211 ENDIF 2212 2213 ! return surface radiation to horizontal surfaces 2214 ! to rad_sw_in, rad_lw_in and rad_net for outputs 2215 !!!!!!!!!! 2216 ! we need the original radiation on urban top layer 2217 ! for calculation of MRT so we can't do adjustment here for now 2218 !!!!!!!!!! 2219 !!!DO isurf = 1, nsurfl 2220 !!! i = surfl(ix,isurf) 2221 !!! j = surfl(iy,isurf) 2222 !!! k = surfl(iz,isurf) 2223 !!! d = surfl(id,isurf) 2224 !!! IF ( d==iroof ) THEN 2225 !!! rad_sw_in(:,j,i) = surfinsw(isurf) 2226 !!! rad_lw_in(:,j,i) = surfinlw(isurf) 2227 !!! rad_net(j,i) = rad_sw_in(k,j,i)  rad_sw_out(k,j,i) + rad_lw_in(k,j,i)  rad_lw_out(k,j,i) 2228 !!! ENDIF 2229 !!!ENDDO 2230 2231 END SUBROUTINE usm_radiation 2232 2233 2234 2235 !! 2236 ! Description: 2237 !  2238 !> This subroutine splits direct and diffusion dw radiation 2239 !> It sould not be called in case the radiation model already does it 2240 !> It follows <CITATION> 2241 !! 2242 SUBROUTINE usm_calc_diffusion_radiation 2243 2244 REAL(wp), PARAMETER :: sol_const = 1367.0_wp !< solar conbstant 2245 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 2246 INTEGER(iwp) :: i, j 2247 REAL(wp), PARAMETER :: year_seconds = 86400._wp * 365._wp 2248 REAL(wp) :: year_angle !< angle 2249 REAL(wp) :: etr !< extraterestrial radiation 2250 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 2251 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 2252 REAL(wp) :: clearnessIndex !< clearness index 2253 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 2254 2255 2256 ! 2257 ! Calculate current day and time based on the initial values and simulation time 2258 year_angle = ((day_init*86400) + time_utc_init+time_since_reference_point) & 2259 / year_seconds * 2.0_wp * pi 2260 2261 etr = sol_const * (1.00011_wp + & 2262 0.034221_wp * cos(year_angle) + & 2263 0.001280_wp * sin(year_angle) + & 2264 0.000719_wp * cos(2.0_wp * year_angle) + & 2265 0.000077_wp * sin(2.0_wp * year_angle)) 2266 2267 ! 2268 ! Under a very low angle, we keep extraterestrial radiation at 2269 ! the last small value, therefore the clearness index will be pushed 2270 ! towards 0 while keeping full continuity. 2271 ! 2272 IF ( zenith(0) <= lowest_solarUp ) THEN 2273 corrected_solarUp = lowest_solarUp 2274 ELSE 2275 corrected_solarUp = zenith(0) 2276 ENDIF 2277 2278 horizontalETR = etr * corrected_solarUp 2279 2280 DO i = nxlg, nxrg 2281 DO j = nysg, nyng 2282 clearnessIndex = rad_sw_in(0,j,i) / horizontalETR 2283 diff_frac = 1.0_wp / (1.0_wp + exp(5.0033_wp + 8.6025_wp * clearnessIndex)) 2284 rad_sw_in_diff(j,i) = rad_sw_in(0,j,i) * diff_frac 2285 rad_sw_in_dir(j,i) = rad_sw_in(0,j,i) * (1.0_wp  diff_frac) 2286 rad_lw_in_diff(j,i) = rad_lw_in(0,j,i) 2287 ENDDO 2288 ENDDO 2289 2290 END SUBROUTINE usm_calc_diffusion_radiation 2291 2292 2293 !! 2294 ! Description: 2295 !  2296 !> Solver for the energy balance at the ground/roof/wall surface. 2297 !> It follows basic ideas and structure of lsm_energy_balance 2298 !> with many simplifications and adjustments. 2299 !> TODO better description 2300 !! 2301 SUBROUTINE usm_surface_energy_balance 2302 2303 IMPLICIT NONE 2304 2305 INTEGER(iwp) :: i, j, k, l, d ! running indices 2306 2307 REAL(wp) :: pt1 ! temperature at first grid box adjacent to surface 2308 REAL(wp) :: u1,v1,w1 ! near wall u,v,w 2309 REAL(wp) :: stend ! surface tendency 2310 REAL(wp) :: coef_1 ! first coeficient for prognostic equation 2311 REAL(wp) :: coef_2 ! second coeficient for prognostic equation 2312 REAL(wp) :: rho_cp ! rho_wall_surface * cp 2313 REAL(wp) :: r_a ! aerodynamic resistance for horizontal and vertical surfaces 2314 REAL(wp) :: f_shf ! factor for shf_eb 2315 REAL(wp) :: lambda_surface ! current value of lambda_surface (heat conductivity between air and wall) 2316 REAL(wp) :: Ueff ! effective wind speed for calculation of heat transfer coefficients 2317 REAL(wp) :: httc ! heat transfer coefficient 2318 REAL(wp), DIMENSION(nzub:nzut) :: exn ! value of the Exner function in layers 2319 2320 REAL(wp), DIMENSION(0:4) :: dxdir ! surface normal direction gridbox length 2321 REAL(wp) :: dtime ! simulated time of day (in UTC) 2322 INTEGER(iwp) :: dhour ! simulated hour of day (in UTC) 2323 REAL(wp) :: acoef ! actual coefficient of diurnal profile of anthropogenic heat 2324 2325 dxdir = (/dz,dy,dy,dx,dx/) 2326 2327 exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp !< Exner function 2328 2329 ! 2330 DO l = startenergy, endenergy 2331 ! Calculate frequently used parameters 2332 d = surfl(id,l) 2333 k = surfl(iz,l) 2334 j = surfl(iy,l) 2335 i = surfl(ix,l) 2336 2337 ! TODO  how to calculate lambda_surface for horizontal surfaces 2338 ! (lambda_surface is set according to stratification in land surface model) 2339 IF ( ol(j,i) >= 0.0_wp ) THEN 2340 lambda_surface = lambda_surf(l) 2341 ELSE 2342 lambda_surface = lambda_surf(l) 2343 ENDIF 2344 2345 pt1 = pt(k,j,i) 2346 2347 ! calculate rho * cp coefficient at surface layer 2348 rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) ) 2349 2350 ! calculate aerodyamic resistance. 2351 IF ( d == iroof ) THEN 2352 ! calculation for horizontal surfaces follows LSM formulation 2353 ! pt, us, ts are not available for the prognostic time step, 2354 ! data from the last time step is used here. 2355 2356 r_a = (pt1  t_surf(l)/exn(k)) / (ts(j,i) * us(j,i) + 1.0E10_wp) 2357 2358 ! make sure that the resistance does not drop to zero 2359 IF ( ABS(r_a) < 1.0E10_wp ) r_a = 1.0E10_wp 2360 2361 !!!!!!!!!!!!!!!!!!!! 2362 ! the parameterization is developed originally for larger scales 2363 ! (compare with remark in TUF3D) 2364 ! our first experiences show that the parameterization underestimates 2365 ! r_a in meter resolution. 2366 ! temporary solution  multiplication by magic constant :(. 2367 r_a = r_a * ra_horiz_coef 2368 2369 ! factor for shf_eb 2370 f_shf = rho_cp / r_a 2371 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2372 WRITE(9,'(f8.1,2i3,a,4i3,100000g20.5)') simulated_time,intermediate_timestep_count,myid,'wall_f_shf ', & 2373 i,j,k,d,f_shf, r_a,ra_horiz_coef,pt1,t_surf(l)/exn(k), ts(j,i), us(j,i) 2374 FLUSH(9) 2375 ENDIF 2376 ELSE 2377 2378 ! calculation of r_a for vertical surfaces 2379 ! 2380 ! heat transfer coefficient for forced convection along vertical walls 2381 ! follows formulation in TUF3d model (Krayenhoff & Voogt, 2006) 2382 ! 2383 ! H = httc (Tsfc  Tair) 2384 ! httc = rw * (11.8 + 4.2 * Ueff)  4.0 2385 ! 2386 ! rw: wall patch roughness relative to 1.0 for concrete 2387 ! Ueff: effective wind speed 2388 !  4.0 is a reduction of Rowley et al (1930) formulation based on 2389 ! Cole and Sturrock (1977) 2390 ! 2391 ! Ucan: Canyon wind speed 2392 ! wstar: convective velocity 2393 ! Qs: surface heat flux 2394 ! zH: height of the convective layer 2395 ! wstar = (g/Tcan*Qs*zH)**(1./3.) 2396 2397 ! staggered grid needs to be taken into consideration 2398 IF ( d == inorth ) THEN 2399 u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp 2400 v1 = v(k,j+1,i) 2401 ELSE IF ( d == isouth ) THEN 2402 u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp 2403 v1 = v(k,j,i) 2404 ELSE IF ( d == ieast ) THEN 2405 u1 = u(k,j,i+1) 2406 v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp 2407 ELSE IF ( d == iwest ) THEN 2408 u1 = u(k,j,i) 2409 v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp 2410 ELSE 2411 STOP 2412 ENDIF 2413 w1 = (w(k,j,i)+w(k1,j,i))*0.5_wp 2414 2415 Ueff = SQRT(u1**2 + v1**2 + w1**2) 2416 2417 httc = roughness_wall(l) * (11.8 + 4.2 * Ueff)  4.0 2418 2419 f_shf = httc 2420 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2421 WRITE(9,'(f8.1,2i3,a,4i3,100000g20.5)') simulated_time,intermediate_timestep_count,myid,'wall_f_shf ', & 2422 i,j,k,d,f_shf, Ueff, roughness_wall(l), pt1, t_surf(l)/exn(k) 2423 FLUSH(9) 2424 ENDIF 2425 ENDIF 2426 2427 ! add LW up so that it can be removed in prognostic equation 2428 rad_net_l(l) = surfinsw(l)  surfoutsw(l) + surfinlw(l)  surfoutlw(l) 2429 2430 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2431 WRITE(9,'(f8.1,2i3,a,4i3,100000f20.5)') simulated_time, intermediate_timestep_count, myid, 'wallrad ', & 2432 i,j,k,d, rad_net_l(l), surfinsw(l), surfoutsw(l), surfinlw(l), surfoutlw(l), t_surf(l) 2433 FLUSH(9) 2434 ENDIF 2435 ! numerator of the prognostic equation 2436 coef_1 = rad_net_l(l) + & !!!! coef +1 corresponds to lwout included in calculation of radnet_l 2437 (3.0_wp+1.0_wp) * emiss_surf(l) * sigma_sb * t_surf(l) ** 4 + & 2438 f_shf * pt1 + & 2439 lambda_surface * t_wall(nzb_wall,l) 2440 2441 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2442 WRITE(9,'(f8.1,2i3,a,4i3,100000g20.5)') simulated_time,intermediate_timestep_count,myid,'wallcoef1 ', & 2443 i,j,k,d, coef_1,rad_net_l(l),emiss_surf(l),sigma_sb, t_surf(l), & 2444 f_shf, pt1, lambda_surface, t_wall(nzb_wall,l) 2445 FLUSH(9) 2446 ENDIF 2447 2448 ! denominator of the prognostic equation 2449 coef_2 = 4.0_wp * emiss_surf(l) * sigma_sb * t_surf(l) ** 3 & 2450 + lambda_surface + f_shf / exn(k) 2451 2452 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2453 WRITE(9,'(f8.1,2i3,a,4i3,100000g20.5)') simulated_time,intermediate_timestep_count,myid,'wallcoef2 ',& 2454 i,j,k,d, coef_2, emiss_surf(l), sigma_sb, t_surf(l), lambda_surface, f_shf 2455 FLUSH(9) 2456 ENDIF 2457 2458 ! implicit solution when the surface layer has no heat capacity, 2459 ! otherwise use RK3 scheme. 2460 t_surf_p(l) = ( coef_1 * dt_3d * tsc(2) + c_surface(l) * t_surf(l) ) / & 2461 ( c_surface(l) + coef_2 * dt_3d * tsc(2) ) 2462 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2463 WRITE(9,'(f8.1,2i3,a,4i3,100000g20.5)') simulated_time,intermediate_timestep_count,myid,'walltsurf ', & 2464 i,j,k,d, t_surf_p(l), coef_1, dt_3d, tsc(2), c_surface(l), t_surf(l), coef_2, dt_3d 2465 FLUSH(9) 2466 ENDIF 2467 2468 ! add RK3 term 2469 t_surf_p(l) = t_surf_p(l) + dt_3d * tsc(3) * tt_surface_m(l) 2470 2471 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2472 WRITE(9,'(f8.1,2i3,a,4i3,100000f20.5)') simulated_time,intermediate_timestep_count,myid,'t_surf_p_wall3 ', & 2473 i,j,k,d,t_surf_p(l), dt_3d, tsc(3), tt_surface_m(l) 2474 FLUSH(9) 2475 ENDIF 2476 2477 ! calculate true tendency 2478 stend = (t_surf_p(l)  t_surf(l)  dt_3d * tsc(3) * tt_surface_m(l)) / (dt_3d * tsc(2)) 2479 2480 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2481 WRITE(9,'(f8.1,2i3,a,4i3,100000f20.5)') simulated_time,intermediate_timestep_count,myid,'t_surf_p_tend ', & 2482 i,j,k,d,stend,t_surf_p(l),t_surf(l), dt_3d, tsc(3), tt_surface_m(l), tsc(2) 2483 FLUSH(9) 2484 ENDIF 2485 2486 ! calculate t_surf tendencies for the next RungeKutta step 2487 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2488 IF ( intermediate_timestep_count == 1 ) THEN 2489 tt_surface_m(l) = stend 2490 ELSEIF ( intermediate_timestep_count < & 2491 intermediate_timestep_count_max ) THEN 2492 tt_surface_m(l) = 9.5625_wp * stend + 5.3125_wp & 2493 * tt_surface_m(l) 2494 ENDIF 2495 ENDIF 2496 2497 ! in case of fast changes in the skin temperature, it is required to 2498 ! update the radiative fluxes in order to keep the solution stable 2499 IF ( ABS( t_surf_p(l)  t_surf(l) ) > 1.0_wp ) THEN 2500 force_radiation_call_l = .TRUE. 2501 ENDIF 2502 2503 ! for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf. 2504 ! there is no equivalent surface gridpoint for vertical surfaces. 2505 ! pt(k,j,i) is calculated for all directions in diffusion_s 2506 ! using surface and wall heat fluxes 2507 IF ( d == iroof ) THEN 2508 pt(nzb_s_inner(j,i),j,i) = t_surf_p(l) / exn(k) 2509 ENDIF 2510 2511 ! calculate fluxes 2512 ! rad_net_l is never used!! 2513 rad_net_l(l) = rad_net_l(l) + 3.0_wp * sigma_sb & 2514 * t_surf(l)**4  4.0_wp * sigma_sb & 2515 * t_surf(l)**3 * t_surf_p(l) 2516 wghf_eb(l) = lambda_surface * (t_surf_p(l)  t_wall(nzb_wall,l)) 2517 2518 ! ground/wall/roof surface heat flux 2519 wshf_eb(l) =  f_shf * ( pt1  t_surf_p(l) ) 2520 2521 ! store kinematic surface heat fluxes for utilization in other processes 2522 ! diffusion_s, surface_layer_fluxes,... 2523 IF ( d == iroof ) THEN 2524 ! shf is used in diffusion_s and also 2525 ! for calculation of surface layer fluxes 2526 ! update for horizontal surfaces 2527 shf(j,i) = wshf_eb(l) / rho_cp 2528 ELSE 2529 ! surface heat flux for vertical surfaces 2530 ! used in diffusion_s 2531 wshf(l) = wshf_eb(l) / rho_cp 2532 ENDIF 2533 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2534 WRITE(9,'(f8.1,2i3,a,4i3,100000f20.5)') simulated_time, & 2535 intermediate_timestep_count, myid, 'shf2 ', & 2536 i,j,k,d, pt1, t_surf_p(l), stend, rho_cp, wshf_eb(l)/rho_cp, wshf_eb(l), & 2537 wghf_eb(l), lambda_surface, t_wall(nzb_wall,l) 2538 FLUSH(9) 2539 ENDIF 2540 2541 ENDDO 2542 2543 2544 IF ( usm_anthropogenic_heat .AND. & 2545 intermediate_timestep_count == intermediate_timestep_count_max ) THEN 2546 ! application of the additional anthropogenic heat sources 2547 ! we considere the traffic for now so all heat is absorbed 2548 ! to the first layer, generalization would be worth 2549 2550 ! calculation of actual profile coefficient 2551 ! ??? check time_since_reference_point ??? 2552 dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp) 2553 dhour = INT(dtime/3600.0_wp) 2554 ! linear interpolation of coeficient 2555 acoef = (REAL(dhour+1,wp)dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wpREAL(dhour,wp))*aheatprof(dhour+1) 2556 DO i = nxl, nxr 2557 DO j = nys, nyn 2558 IF ( aheat(j,i) > 0.0_wp ) THEN 2559 ! TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d 2560 ! given to anthropogenic heat aheat*acoef (W*m2) 2561 ! k = nzb_s_inner(j,i)+1 2562 ! pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz) 2563 ! Instead of this, we can adjust shf in case AH only at surface 2564 shf(j,i) = shf(j,i) + aheat(j,i)*acoef * ddx * ddy / rho_cp 2565 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2566 WRITE(9,'(f8.1,i3,a,3i3,100000f20.5)') simulated_time, myid, 'ah1 ', & 2567 i,j,k,shf(j,i), aheat(j,i)*acoef*ddx*ddy/rho_cp, & 2568 aheat(j,i),acoef, ddx,ddy,rho_cp 2569 FLUSH(9) 2570 ENDIF 2571 ENDIF 2572 ENDDO 2573 ENDDO 2574 ENDIF 2575 2576 ! pt and shf are defined on nxlg:nxrg,nysg:nyng 2577 ! get the borders from neighbours 2578 CALL exchange_horiz( pt, nbgp ) 2579 CALL exchange_horiz_2d( shf ) 2580 2581 2582 ! calculation of force_radiation_call: 2583 ! Make logical OR for all processes. 2584 ! Force radiation call if at least one processor forces it. 2585 IF ( intermediate_timestep_count == intermediate_timestep_count_max1 ) & 2586 THEN 2587 #if defined( __parallel ) 2588 IF ( collective_wait ) CALL mpi_barrier( comm2d, ierr ) 2589 CALL mpi_allreduce( force_radiation_call_l, force_radiation_call, & 2590 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 2591 #else 2592 force_radiation_call = force_radiation_call_l 2593 #endif 2594 force_radiation_call_l = .FALSE. 2595 ENDIF 2596 2597 2598 END SUBROUTINE usm_surface_energy_balance 2599 2600 3149 2601 3150 !! 2602 3151 ! Description: 2603 3152 !  2604 3153 ! 2605 !> Wall model as part of the urban surface model. The model predicts wall 2606 !> temperature. 2607 !! 2608 SUBROUTINE usm_material_heat_model 2609 2610 2611 IMPLICIT NONE 2612 2613 INTEGER(iwp) :: i,j,k,l,kw !< running indices 2614 2615 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency 2616 2617 2618 DO l = startenergy, endenergy 2619 ! calculate frequently used parameters 2620 k = surfl(iz,l) 2621 j = surfl(iy,l) 2622 i = surfl(ix,l) 2623 2624 ! 2625 ! prognostic equation for ground/wall/roof temperature t_wall 2626 wtend(:) = 0.0_wp 2627 wtend(nzb_wall) = (1.0_wp/rho_c_wall(nzb_wall,l)) * & 2628 ( lambda_h(nzb_wall,l) * ( t_wall(nzb_wall+1,l) & 2629  t_wall(nzb_wall,l) ) * ddz_wall(nzb_wall+1,l) & 2630 + wghf_eb(l) ) * ddz_wall_stag(nzb_wall,l) 2631 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2632 WRITE(9,'(f8.1,2i3,a,4i3,100000g12.5)') simulated_time,intermediate_timestep_count, & 2633 myid,'wallmodel1b ', & 2634 i,j,k,nzb_wall,wtend(nzb_wall), & 2635 t_wall(nzb_wall+1,l), t_wall(nzb_wall,l), & 2636 wghf_eb(l), & 2637 dz_wall(nzb_wall,l), dz_wall_stag(nzb_wall,l), & 2638 lambda_h(nzb_wall,l) 2639 FLUSH(9) 3154 !> This subroutine is part of the urban surface model. 3155 !> It reads daily heat produced by anthropogenic sources 3156 !> and the diurnal cycle of the heat. 3157 !! 3158 SUBROUTINE usm_read_anthropogenic_heat 3159 3160 INTEGER(iwp) :: i,j,ii 3161 REAL(wp) :: heat 3162 3163 ! allocation of array of sources of anthropogenic heat and their diural profile 3164 ALLOCATE( aheat(nys:nyn,nxl:nxr) ) 3165 ALLOCATE( aheatprof(0:24) ) 3166 3167 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3168 ! read daily amount of heat and its daily cycle 3169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3170 aheat = 0.0_wp 3171 DO ii = 0, io_blocks1 3172 IF ( ii == io_group ) THEN 3173 3174 ! open anthropogenic heat file 3175 OPEN( 151, file='ANTHROPOGENIC_HEAT'//TRIM(coupling_char), action='read', & 3176 status='old', form='formatted', err=11 ) 3177 i = 0 3178 j = 0 3179 DO 3180 READ( 151, *, err=12, end=13 ) i, j, heat 3181 IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN 3182 ! write heat into the array 3183 aheat(j,i) = heat 3184 ENDIF 3185 CYCLE 3186 12 WRITE(message_string,'(a,2i4)') 'error in file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' after line ',i,j 3187 CALL message( 'usm_read_anthropogenic_heat', 'PA0515', 0, 1, 0, 6, 0 ) 3188 ENDDO 3189 13 CLOSE(151) 3190 CYCLE 3191 11 message_string = 'file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' does not exist' 3192 CALL message( 'usm_read_anthropogenic_heat', 'PA0516', 1, 2, 0, 6, 0 ) 2640 3193 ENDIF 2641 3194 2642 DO kw = nzb_wall+1, nzt_wall 2643 wtend(kw) = (1.0_wp/rho_c_wall(kw,l)) & 2644 * ( lambda_h(kw,l) & 2645 * ( t_wall(kw+1,l)  t_wall(kw,l) ) & 2646 * ddz_wall(kw+1,l) & 2647  lambda_h(kw1,l) & 2648 * ( t_wall(kw,l)  t_wall(kw1,l) ) & 2649 * ddz_wall(kw,l) & 2650 ) * ddz_wall_stag(kw,l) 2651 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2652 WRITE(9,'(f8.1,2i3,a,4i3,100000g12.5)') simulated_time,intermediate_timestep_count,& 2653 myid,'wallmodel1c ', & 2654 i,j,k,kw,wtend(kw), & 2655 t_wall(kw+1,l), t_wall(kw,l), t_wall(kw1,l), & 2656 dz_wall(kw,l), dz_wall_stag(kw,l), & 2657 lambda_h(kw,l) 2658 FLUSH(9) 2659 ENDIF 2660 ENDDO 2661 2662 t_wall_p(nzb_wall:nzt_wall,l) = t_wall(nzb_wall:nzt_wall,l) & 2663 + dt_3d * ( tsc(2) & 2664 * wtend(nzb_wall:nzt_wall) + tsc(3) & 2665 * tt_wall_m(nzb_wall:nzt_wall,l) ) 2666 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2667 WRITE(9,'(f8.1,2i3,a,3i3,100000g15.5)') simulated_time,intermediate_timestep_count,& 2668 myid,'wallmodel2 ', & 2669 i,j,k, & 2670 t_wall_p(nzb_wall,l),t_wall(nzb_wall,l), & 2671 wtend(nzb_wall), tt_wall_m(nzb_wall,l) 2672 FLUSH(9) 3195 #if defined( __parallel ) && ! defined ( __check ) 3196 CALL mpi_barrier( comm2d, ierr ) 3197 #endif 3198 ENDDO 3199 3200 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3201 ! read diurnal profiles of heat sources 3202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3203 aheatprof = 0.0_wp 3204 DO ii = 0, io_blocks1 3205 IF ( ii == io_group ) THEN 3206 3207 ! open anthropogenic heat profile file 3208 OPEN( 151, file='ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char), action='read', & 3209 status='old', form='formatted', err=21 ) 3210 i = 0 3211 DO 3212 READ( 151, *, err=22, end=23 ) i, heat 3213 IF ( i >= 0 .AND. i <= 24 ) THEN 3214 ! write heat into the array 3215 aheatprof(i) = heat 3216 ENDIF 3217 CYCLE 3218 22 WRITE(message_string,'(a,i4)') 'error in file ANTHROPOGENIC_HEAT_PROFILE'// & 3219 TRIM(coupling_char)//' after line ',i 3220 CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 ) 3221 ENDDO 3222 aheatprof(24) = aheatprof(0) 3223 23 CLOSE(151) 3224 CYCLE 3225 21 message_string = 'file ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char)//' does not exist' 3226 CALL message( 'usm_read_anthropogenic_heat', 'PA0518', 1, 2, 0, 6, 0 ) 2673 3227 ENDIF 2674 3228 2675 ! 2676 ! calculate t_wall tendencies for the next RungeKutta step 2677 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2678 IF ( intermediate_timestep_count == 1 ) THEN 2679 DO kw = nzb_wall, nzt_wall 2680 tt_wall_m(kw,l) = wtend(kw) 2681 ENDDO 2682 ELSEIF ( intermediate_timestep_count < & 2683 intermediate_timestep_count_max ) THEN 2684 DO kw = nzb_wall, nzt_wall 2685 tt_wall_m(kw,l) = 9.5625_wp * wtend(kw) + 5.3125_wp & 2686 * tt_wall_m(kw,l) 2687 ENDDO 2688 ENDIF 2689 ENDIF 3229 #if defined( __parallel ) && ! defined ( __check ) 3230 CALL mpi_barrier( comm2d, ierr ) 3231 #endif 2690 3232 ENDDO 2691 2692 END SUBROUTINE usm_material_heat_model 2693 2694 !! 3233 3234 END SUBROUTINE usm_read_anthropogenic_heat 3235 3236 3237 !! 3238 ! 2695 3239 ! Description: 2696 3240 !  3241 !> Soubroutine reads t_surf and t_wall data from restart files 3242 !kanani: Renamed this routine according to corresponging routines in PALM 3243 !kanani: Modified the routine to match read_var_list, from where usm_read_restart_data 3244 ! shall be called in the future. This part has not been tested yet. (see virtual_flight_mod) 3245 ! Also, I had some trouble with the allocation of t_surf, since this is a pointer. 3246 ! So, I added some directives here. 3247 !! 3248 SUBROUTINE usm_read_restart_data 3249 3250 3251 IMPLICIT NONE 3252 3253 CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string 3254 3255 INTEGER :: i !< running index 3256 3257 3258 DO i = 0, io_blocks1 3259 IF ( i == io_group ) THEN 3260 READ ( 13 ) variable_chr 3261 DO WHILE ( TRIM( variable_chr ) /= '*** end usm ***' ) 3262 3263 SELECT CASE ( TRIM( variable_chr ) ) 3264 3265 CASE ( 't_surf' ) 3266 #if defined( __nopointer ) 3267 IF ( .NOT. ALLOCATED( t_surf ) ) & 3268 ALLOCATE( t_surf(startenergy:endenergy) ) 3269 READ ( 13 ) t_surf 3270 #else 3271 IF ( .NOT. ALLOCATED( t_surf_1 ) ) & 3272 ALLOCATE( t_surf_1(startenergy:endenergy) ) 3273 READ ( 13 ) t_surf_1 3274 #endif 3275 3276 CASE ( 't_wall' ) 3277 #if defined( __nopointer ) 3278 IF ( .NOT. ALLOCATED( t_wall ) ) & 3279 ALLOCATE( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 3280 READ ( 13 ) t_wall 3281 #else 3282 IF ( .NOT. ALLOCATED( t_wall_1 ) ) & 3283 ALLOCATE( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 3284 READ ( 13 ) t_wall_1 3285 #endif 3286 3287 CASE DEFAULT 3288 WRITE ( message_string, * ) 'unknown variable named "', & 3289 TRIM( variable_chr ), '" found in', & 3290 '&data from prior run on PE ', myid 3291 CALL message( 'user_read_restart_data', 'UI0012', 1, 2, 0, 6, 0 ) 3292 3293 END SELECT 3294 3295 READ ( 13 ) variable_chr 3296 3297 ENDDO 3298 ENDIF 3299 #if defined( __parallel ) 3300 CALL MPI_BARRIER( comm2d, ierr ) 3301 #endif 3302 ENDDO 3303 3304 END SUBROUTINE usm_read_restart_data 3305 3306 3307 !! 2697 3308 ! 2698 !> This function applies the kinematic wall heat fluxes2699 !> for walls in four directions for all gridboxes in urban layer.2700 !> It is called out from subroutine prognostic_equations.2701 !> TODO Compare performance with cycle runnig l=startwall,endwall...2702 !!2703 SUBROUTINE usm_wall_heat_flux2704 2705 IMPLICIT NONE2706 2707 INTEGER(iwp) :: i,j,k,d,l !< running indices2708 2709 ! DO i = nxl, nxr2710 ! DO j = nys, nyn2711 ! DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)2712 ! DO d = 1,42713 ! l = gridsurf(d,k,j,i)2714 ! IF ( l /= 0 ) THEN2715 ! tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)2716 ! ENDIF2717 ! ENDDO2718 ! ENDDO2719 ! ENDDO2720 ! ENDDO2721 2722 !ketelsen: To spare the 4 Dimensional array gridsurf, the complete lloop is executed for every i and j value2723 !ketelsen: This is not very elegant, but also not time critical2724 DO l = startenergy, endenergy2725 j = surfl(iy,l)2726 i = surfl(ix,l)2727 k = surfl(iz,l)2728 d = surfl(id,l)2729 tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)2730 ENDDO2731 2732 END SUBROUTINE usm_wall_heat_flux2733 2734 2735 !!2736 3309 ! Description: 2737 3310 !  2738 ! 2739 !> This function applies the kinematic wall heat fluxes 2740 !> for walls in four directions around the gridbox i,j. 2741 !> It is called out from subroutine prognostic_equations. 2742 !! 2743 SUBROUTINE usm_wall_heat_flux_ij(i,j) 2744 3311 !> Soubroutine reads svf and svfsurf data from saved file 3312 !! 3313 SUBROUTINE usm_read_svf_from_file 3314 2745 3315 IMPLICIT NONE 2746 2747 INTEGER(iwp), INTENT(in) :: i,j !< indices of grid box 2748 INTEGER(iwp) :: ii,jj,k,d,l 2749 2750 ! DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i) 2751 ! DO d = 1,4 2752 ! l = gridsurf(d,k,j,i) 2753 ! IF ( l /= 0 ) THEN 2754 ! tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d) 2755 ! ENDIF 2756 ! ENDDO 2757 ! ENDDO 2758 2759 2760 !ketelsen: To spare the 4 Dimensional array gridsurf, the complete lloop is executed for every i and j value 2761 !ketelsen: This is not very elegant, but also not time critical 2762 2763 DO l = startenergy, endenergy 2764 jj = surfl(iy,l) 2765 ii = surfl(ix,l) 2766 IF( ii == i .AND. jj == j) then 2767 k = surfl(iz,l) 2768 IF(k >= nzb_s_inner(j,i)+1 .AND. k <= nzb_s_outer(j,i)) then 2769 d = surfl(id,l) 2770 if(d >= 1 .and. d <= 4) then 2771 tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d) 2772 end if 2773 ENDIF 3316 INTEGER :: fsvf = 89 3317 INTEGER :: i 3318 CHARACTER(usm_version_len) :: usm_version_field 3319 CHARACTER(svf_code_len) :: svf_code_field 3320 3321 DO i = 0, io_blocks1 3322 IF ( i == io_group ) THEN 3323 OPEN ( fsvf, file=TRIM(svf_file_name)//TRIM(coupling_char)//myid_char, & 3324 form='unformatted', status='old' ) 3325 3326 ! read and check version 3327 READ ( fsvf ) usm_version_field 3328 IF ( TRIM(usm_version_field) /= TRIM(usm_version) ) THEN 3329 WRITE( message_string, * ) 'Version of binary SVF file "', & 3330 TRIM(usm_version_field), '" does not match ', & 3331 'the version of model "', TRIM(usm_version), '"' 3332 CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 ) 3333 ENDIF 3334 3335 ! read nsvfcsfl, nsvfl 3336 READ ( fsvf ) nsvfcsfl, nsvfl 3337 IF ( nsvfcsfl <= 0 ) THEN 3338 WRITE( message_string, * ) 'Wrong number of SVF or CSF' 3339 CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 ) 3340 ELSE 3341 WRITE(message_string,*) ' Number of SVF and CSF to read', nsvfcsfl, nsvfl 3342 CALL location_message( message_string, .TRUE. ) 3343 ENDIF 3344 3345 ALLOCATE(svf(ndsvf,nsvfcsfl)) 3346 ALLOCATE(svfsurf(ndsvf,nsvfcsfl)) 3347 3348 READ(fsvf) svf 3349 READ(fsvf) svfsurf 3350 READ ( fsvf ) svf_code_field 3351 3352 IF ( TRIM(svf_code_field) /= TRIM(svf_code) ) THEN 3353 WRITE( message_string, * ) 'Wrong structure of binary svf file' 3354 CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 ) 3355 ENDIF 3356 3357 CLOSE (fsvf) 3358 2774 3359 ENDIF 3360 #if defined( __parallel ) 3361 CALL mpi_barrier( comm2d, ierr ) 3362 #endif 2775 3363 ENDDO 2776 3364 2777 END SUBROUTINE usm_wall_heat_flux_ij 2778 2779 3365 END SUBROUTINE usm_read_svf_from_file 3366 2780 3367 2781 3368 !! … … 2792 3379 REAL(wp), DIMENSION(n_surface_params) :: wtp 2793 3380 2794 INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: us _par2795 REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: us _val3381 INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: usm_par 3382 REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: usm_val 2796 3383 INTEGER(iwp) :: k, l, d, iw, jw, kw, it, ip, ii, ij 2797 3384 INTEGER(iwp) :: i, j … … 2806 3393 REAL(wp) :: wealbedo3, wethick3, snalbedo3, snthick3 2807 3394 2808 2809 !read categories of walls and their parameters2810 3395 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3396 ! read categories of walls and their parameters 3397 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2811 3398 DO ii = 0, io_blocks1 2812 3399 IF ( ii == io_group ) THEN 2813 3400 2814 !open urban surface file3401 ! open urban surface file 2815 3402 OPEN( 151, file='SURFACE_PARAMETERS'//coupling_char, action='read', & 2816 3403 status='old', form='formatted', err=15 ) 2817 !first test and get n_surface_types3404 ! first test and get n_surface_types 2818 3405 k = 0 2819 3406 l = 0 … … 2829 3416 ALLOCATE( surface_type_codes(n_surface_types) ) 2830 3417 ALLOCATE( surface_params(n_surface_params, n_surface_types) ) 2831 !real reading3418 ! real reading 2832 3419 rewind( 151 ) 2833 3420 k = 0 … … 2850 3437 ENDDO 2851 3438 2852 2853 !read types of surfaces2854 2855 us _par = 03439 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3440 ! read types of surfaces 3441 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3442 usm_par = 0 2856 3443 DO ii = 0, io_blocks1 2857 3444 IF ( ii == io_group ) THEN 2858 3445 2859 3446 ! 2860 !open csv urban surface file3447 ! open csv urban surface file 2861 3448 OPEN( 151, file='URBAN_SURFACE'//TRIM(coupling_char), action='read', & 2862 3449 status='old', form='formatted', err=23 ) … … 2865 3452 DO 2866 3453 l = l+1 2867 !i, j, height, nz, roof, dirwe, dirsn, category, soilcat,2868 !weheight1, wecat1, snheight1, sncat1, weheight2, wecat2, snheight2, sncat2,2869 !weheight3, wecat3, snheight3, sncat33454 ! i, j, height, nz, roof, dirwe, dirsn, category, soilcat, 3455 ! weheight1, wecat1, snheight1, sncat1, weheight2, wecat2, snheight2, sncat2, 3456 ! weheight3, wecat3, snheight3, sncat3 2870 3457 READ( 151, *, err=21, end=25 ) i, j, height, nz, roof, dirwe, dirsn, & 2871 3458 category, albedo, thick, & … … 2877 3464 snheight3, sncat3, snalbedo3, snthick3 2878 3465 2879 IF ( i >= nxlg .AND. i <= nxrg .AND. j >= nysg .AND. j <= nyng )THEN2880 !write integer variables into array2881 us _par(:,j,i) = (/1, nz, roof, dirwe, dirsn, category,&3466 IF ( i >= nxlg .AND. i <= nxrg .AND. j >= nysg .AND. j <= nyng ) THEN 3467 ! write integer variables into array 3468 usm_par(:,j,i) = (/1, nz, roof, dirwe, dirsn, category, & 2882 3469 weheight1, wecat1, weheight2, wecat2, weheight3, wecat3, & 2883 3470 snheight1, sncat1, snheight2, sncat2, snheight3, sncat3 /) 2884 !write real values into array2885 us _val(:,j,i) = (/ albedo, thick,&3471 ! write real values into array 3472 usm_val(:,j,i) = (/ albedo, thick, & 2886 3473 wealbedo1, wethick1, wealbedo2, wethick2, & 2887 3474 wealbedo3, wethick3, snalbedo1, snthick1, & … … 2905 3492 2906 3493 ! 2907 !check completeness and formal correctness of the data3494 ! check completeness and formal correctness of the data 2908 3495 DO i = nxlg, nxrg 2909 3496 DO j = nysg, nyng 2910 IF ( us_par(0,j,i) /= 0 .AND. ( & !< incomplete data,supply default values later 2911 us_par(1,j,i) < nzb .OR. & 2912 us_par(1,j,i) > nzt .OR. & !< incorrect height (nz < nzb .OR. nz > nzt) 2913 us_par(2,j,i) < 0 .OR. & 2914 us_par(2,j,i) > 1 .OR. & !< incorrect roof sign 2915 us_par(3,j,i) < nzbnzt .OR. & 2916 us_par(3,j,i) > nztnzb .OR. & !< incorrect westeast wall direction sign 2917 us_par(4,j,i) < nzbnzt .OR. & 2918 us_par(4,j,i) > nztnzb .OR. & !< incorrect southnorth wall direction sign 2919 us_par(6,j,i) < nzb .OR. & 2920 us_par(6,j,i) > nzt .OR. & !< incorrect pedestrian level height for westeast wall 2921 us_par(8,j,i) > nzt .OR. & 2922 us_par(10,j,i) > nzt .OR. & !< incorrect wall or roof level height for westeast wall 2923 us_par(12,j,i) < nzb .OR. & 2924 us_par(12,j,i) > nzt .OR. & !< incorrect pedestrian level height for southnorth wall 2925 us_par(14,j,i) > nzt .OR. & 2926 us_par(16,j,i) > nzt & !< incorrect wall or roof level height for southnorth wall 2927 ) ) THEN 2928 ! incorrect input data 2929 IF ( debug_prints .AND. time_do3d < dt_3d ) THEN 2930 WRITE(9,*) 'Incorrect US input data at i,j=', i,j 2931 WRITE(9,*) 'us_par = ', us_par(:,j,i) 2932 ENDIF 3497 IF ( usm_par(0,j,i) /= 0 .AND. ( & !< incomplete data,supply default values later 3498 usm_par(1,j,i) < nzb .OR. & 3499 usm_par(1,j,i) > nzt .OR. & !< incorrect height (nz < nzb .OR. nz > nzt) 3500 usm_par(2,j,i) < 0 .OR. & 3501 usm_par(2,j,i) > 1 .OR. & !< incorrect roof sign 3502 usm_par(3,j,i) < nzbnzt .OR. & 3503 usm_par(3,j,i) > nztnzb .OR. & !< incorrect westeast wall direction sign 3504 usm_par(4,j,i) < nzbnzt .OR. & 3505 usm_par(4,j,i) > nztnzb .OR. & !< incorrect southnorth wall direction sign 3506 usm_par(6,j,i) < nzb .OR. & 3507 usm_par(6,j,i) > nzt .OR. & !< incorrect pedestrian level height for westeast wall 3508 usm_par(8,j,i) > nzt .OR. & 3509 usm_par(10,j,i) > nzt .OR. & !< incorrect wall or roof level height for westeast wall 3510 usm_par(12,j,i) < nzb .OR. & 3511 usm_par(12,j,i) > nzt .OR. & !< incorrect pedestrian level height for southnorth wall 3512 usm_par(14,j,i) > nzt .OR. & 3513 usm_par(16,j,i) > nzt & !< incorrect wall or roof level height for southnorth wall 3514 ) ) THEN 3515 ! incorrect input data 2933 3516 WRITE (message_string, "(A,2I5)") 'missing or incorrect data in file URBAN_SURFACE'// & 2934 3517 TRIM(coupling_char)//' for i,j=', i,j … … 2936 3519 ENDIF 2937 3520 2938 !IF ( us_alb(j,i) < 0.0_wp .OR. &2939 ! us_alb(j,i) > 1.0_wp ) & ! incorrect albedo of the land or roof2940 !THEN2941 ! incorrect albedo2942 !WRITE (message_string, "(A,2I5)") 'missing or incorrect albedo in file URBAN_SURFACE'//TRIM(coupling_char)//' for i,j=', i,j2943 !CALL message( 'usm_read_urban_surface', 'PA0504', 1, 2, 0, 6, 0 )2944 !ENDIF2945 2946 2947 3521 ENDDO 2948 3522 ENDDO 2949 3523 2950 !assign the surface types to local surface array3524 ! assign the surface types to local surface array 2951 3525 DO l = startenergy, endenergy 2952 3526 … … 2955 3529 j = surfl(iy,l) 2956 3530 i = surfl(ix,l) 2957 IF ( d == iroof ) THEN2958 !horizontal surface  land or roof3531 IF ( d == iroof ) THEN 3532 ! horizontal surface  land or roof 2959 3533 iw = i 2960 3534 jw = j 2961 IF ( us _par(5,jw,iw) == 0 )THEN2962 IF ( zu(kw) >= roof_height_limit ) THEN3535 IF ( usm_par(5,jw,iw) == 0 ) THEN 3536 IF ( zu(kw) >= roof_height_limit ) THEN 2963 3537 isroof_surf(l) = .TRUE. 2964 3538 surface_types(l) = roof_category !< default category for root surface … … 2970 3544 thickness_wall(l) = 1.0_wp 2971 3545 ELSE 2972 IF ( us _par(2,jw,iw)==0 )THEN3546 IF ( usm_par(2,jw,iw)==0 ) THEN 2973 3547 isroof_surf(l) = .FALSE. 2974 3548 thickness_wall(l) = 1.0_wp 2975 3549 ELSE 2976 3550 isroof_surf(l) = .TRUE. 2977 thickness_wall(l) = us _val(2,jw,iw)3551 thickness_wall(l) = usm_val(2,jw,iw) 2978 3552 ENDIF 2979 surface_types(l) = us _par(5,jw,iw)2980 albedo_surf(l) = us _val(1,jw,iw)3553 surface_types(l) = usm_par(5,jw,iw) 3554 albedo_surf(l) = usm_val(1,jw,iw) 2981 3555 ENDIF 2982 3556 ELSE … … 3004 3578 END SELECT 3005 3579 3006 IF ( kw <= us _par(ii,jw,iw) )THEN3007 !pedestrant zone3580 IF ( kw <= usm_par(ii,jw,iw) ) THEN 3581 ! pedestrant zone 3008 3582 isroof_surf(l) = .FALSE. 3009 IF ( us _par(ii+1,jw,iw) == 0 )THEN3583 IF ( usm_par(ii+1,jw,iw) == 0 ) THEN 3010 3584 surface_types(l) = pedestrant_category !< default category for wall surface in pedestrant zone 3011 3585 albedo_surf(l) = 1.0_wp 3012 3586 thickness_wall(l) = 1.0_wp 3013 3587 ELSE 3014 surface_types(l) = us _par(ii+1,jw,iw)3015 albedo_surf(l) = us _val(ij,jw,iw)3016 thickness_wall(l) = us _val(ij+1,jw,iw)3588 surface_types(l) = usm_par(ii+1,jw,iw) 3589 albedo_surf(l) = usm_val(ij,jw,iw) 3590 thickness_wall(l) = usm_val(ij+1,jw,iw) 3017 3591 ENDIF 3018 ELSE IF ( kw <= us _par(ii+2,jw,iw) )THEN3019 !wall zone3592 ELSE IF ( kw <= usm_par(ii+2,jw,iw) ) THEN 3593 ! wall zone 3020 3594 isroof_surf(l) = .FALSE. 3021 IF ( us _par(ii+3,jw,iw) == 0 )THEN3595 IF ( usm_par(ii+3,jw,iw) == 0 ) THEN 3022 3596 surface_types(l) = wall_category !< default category for wall surface 3023 3597 albedo_surf(l) = 1.0_wp 3024 3598 thickness_wall(l) = 1.0_wp 3025 3599 ELSE 3026 surface_types(l) = us _par(ii+3,jw,iw)3027 albedo_surf(l) = us _val(ij+2,jw,iw)3028 thickness_wall(l) = us _val(ij+3,jw,iw)3600 surface_types(l) = usm_par(ii+3,jw,iw) 3601 albedo_surf(l) = usm_val(ij+2,jw,iw) 3602 thickness_wall(l) = usm_val(ij+3,jw,iw) 3029 3603 ENDIF 3030 ELSE IF ( kw <= us _par(ii+4,jw,iw) )THEN3031 !roof zone3604 ELSE IF ( kw <= usm_par(ii+4,jw,iw) ) THEN 3605 ! roof zone 3032 3606 isroof_surf(l) = .TRUE. 3033 IF ( us _par(ii+5,jw,iw) == 0 )THEN3607 IF ( usm_par(ii+5,jw,iw) == 0 ) THEN 3034 3608 surface_types(l) = roof_category !< default category for roof surface 3035 3609 albedo_surf(l) = 1.0_wp 3036 3610 thickness_wall(l) = 1.0_wp 3037 3611 ELSE 3038 surface_types(l) = us _par(ii+5,jw,iw)3039 albedo_surf(l) = us _val(ij+4,jw,iw)3040 thickness_wall(l) = us _val(ij+5,jw,iw)3612 surface_types(l) = usm_par(ii+5,jw,iw) 3613 albedo_surf(l) = usm_val(ij+4,jw,iw) 3614 thickness_wall(l) = usm_val(ij+5,jw,iw) 3041 3615 ENDIF 3042 3616 ELSE 3043 ! something wrong 3044 !WRITE (message_string, "(A,3I5)") 'incorrect data in file URBAN_SURFACE'//TRIM(coupling_char)//' for i,j,k=', iw,jw,kw 3045 !WRITE(6,'(100i4)'), myid, l, d, ii, iw,jw,kw,us_par(:,jw,iw) 3046 !WRITE(6,'(i4,a)'), myid, message_string 3047 !FLUSH(6) 3617 ! something wrong 3048 3618 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 ) 3049 3619 ENDIF 3050 3620 ENDIF 3051 3621 3052 !find the type position3622 ! find the type position 3053 3623 it = surface_types(l) 3054 3624 ip = 99999 3055 3625 DO k = 1, n_surface_types 3056 IF ( surface_type_codes(k) == it ) THEN3626 IF ( surface_type_codes(k) == it ) THEN 3057 3627 ip = k 3058 3628 EXIT 3059 3629 ENDIF 3060 3630 ENDDO 3061 IF ( ip == 99999 ) THEN 3062 !PRINT*, myid, l, d, iw,jw,kw,it,n_surface_types,surface_type_codes 3063 !FLUSH(6) 3064 ! wall category not found 3631 IF ( ip == 99999 ) THEN 3632 ! wall category not found 3065 3633 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, ' not found for i,j,k=', iw,jw,kw 3066 3634 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 ) 3067 3635 ENDIF 3068 3636 3069 !Fill out the parameters of the wall3070 !wall surface:3637 ! Fill out the parameters of the wall 3638 ! wall surface: 3071 3639 3072 !albedo3073 IF ( albedo_surf(l) < 0.0_wp ) THEN3640 ! albedo 3641 IF ( albedo_surf(l) < 0.0_wp ) THEN 3074 3642 albedo_surf(l) = surface_params(ialbedo, ip) 3075 3643 ENDIF 3076 3644 3077 !emissivity of the wall3645 ! emissivity of the wall 3078 3646 emiss_surf(l) = surface_params(iemiss, ip) 3079 3647 3080 !heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )3648 ! heat conductivity Î»S between air and wall ( W mâ2 Kâ1 ) 3081 3649 lambda_surf(l) = surface_params(ilambdas, ip) 3082 !PRINT*, myid, 'lambda_surf=', l, lambda_surf(l)3083 3650 3084 !roughness relative to concrete3651 ! roughness relative to concrete 3085 3652 roughness_wall(l) = surface_params(irough, ip) 3086 3653 3087 !Surface skin layer heat capacity (J mâ2 Kâ1 )3654 ! Surface skin layer heat capacity (J mâ2 Kâ1 ) 3088 3655 c_surface(l) = surface_params(icsurf, ip) 3089 3656 3090 !wall material parameters:3657 ! wall material parameters: 3091 3658 3092 !thickness of the wall (m)3093 !missing values are replaced by default value for category3094 IF ( thickness_wall(l) <= 0.001_wp ) THEN3659 ! thickness of the wall (m) 3660 ! missing values are replaced by default value for category 3661 IF ( thickness_wall(l) <= 0.001_wp ) THEN 3095 3662 thickness_wall(l) = surface_params(ithick, ip) 3096 3663 ENDIF 3097 3664 3098 !volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )3665 ! volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 ) 3099 3666 rho_c_wall(:,l) = surface_params(irhoC, ip) 3100 3667 3101 !thermal conductivity Î»H of the wall (W mâ1 Kâ1 )3668 ! thermal conductivity Î»H of the wall (W mâ1 Kâ1 ) 3102 3669 lambda_h(:,l) = surface_params(ilambdah, ip) 3103 3670 3104 !IF ( debug_prints .AND. time_do3d < dt_3d ) THEN3105 ! WRITE(9,*) myid, 'readsurface3', l, d, albedo_surf(l), emiss_surf(l), lambda_surf(l), &3106 ! roughness_wall(l), c_surface(l), thickness_wall(l), rho_c_wall(0,l), lambda_h(0,l)3107 ! FLUSH(9)3108 !ENDIF3109 3671 ENDDO 3110 3672 … … 3113 3675 END SUBROUTINE usm_read_urban_surface_types 3114 3676 3115 3677 3116 3678 !! 3117 3679 ! Description: 3118 3680 !  3119 ! 3120 !> This subroutine is part of the urban surface model. 3121 !> It reads daily heat produced by anthropogenic sources 3122 !> and the diurnal cycle of the heat. 3123 !! 3124 SUBROUTINE usm_read_anthropogenic_heat 3125 3126 INTEGER(iwp) :: i,j,ii 3127 REAL(wp) :: heat 3128 3129 ! allocation of array of sources of anthropogenic heat and their diural profile 3130 ALLOCATE( aheat(nys:nyn,nxl:nxr) ) 3131 ALLOCATE( aheatprof(0:24) ) 3132 3133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3134 ! read daily amount of heat and its daily cycle 3135 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3136 aheat = 0.0_wp 3137 DO ii = 0, io_blocks1 3138 IF ( ii == io_group ) THEN 3139 3140 ! open anthropogenic heat file 3141 OPEN( 151, file='ANTHROPOGENIC_HEAT'//TRIM(coupling_char), action='read', & 3142 status='old', form='formatted', err=11 ) 3143 i = 0 3144 j = 0 3145 DO 3146 READ( 151, *, err=12, end=13 ) i, j, heat 3147 IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN 3148 ! write heat into the array 3149 aheat(j,i) = heat 3150 ENDIF 3151 CYCLE 3152 12 WRITE(message_string,'(a,2i4)') 'error in file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' after line ',i,j 3153 CALL message( 'usm_read_anthropogenic_heat', 'PA0515', 0, 1, 0, 6, 0 ) 3154 ENDDO 3155 13 CLOSE(151) 3156 CYCLE 3157 11 message_string = 'file ANTHROPOGENIC_HEAT'//TRIM(coupling_char)//' does not exist' 3158 CALL message( 'usm_read_anthropogenic_heat', 'PA0516', 1, 2, 0, 6, 0 ) 3681 !> Solver for the energy balance at the ground/roof/wall surface. 3682 !> It follows basic ideas and structure of lsm_energy_balance 3683 !> with many simplifications and adjustments. 3684 !> TODO better description 3685 !! 3686 SUBROUTINE usm_surface_energy_balance 3687 3688 IMPLICIT NONE 3689 3690 INTEGER(iwp) :: i, j, k, l, d !< running indices 3691 3692 REAL(wp) :: pt1 !< temperature at first grid box adjacent to surface 3693 REAL(wp) :: u1,v1,w1 !< near wall u,v,w 3694 REAL(wp) :: stend !< surface tendency 3695 REAL(wp) :: coef_1 !< first coeficient for prognostic equation 3696 REAL(wp) :: coef_2 !< second coeficient for prognostic equation 3697 REAL(wp) :: rho_cp !< rho_wall_surface * cp 3698 REAL(wp) :: r_a !< aerodynamic resistance for horizontal and vertical surfaces 3699 REAL(wp) :: f_shf !< factor for shf_eb 3700 REAL(wp) :: lambda_surface !< current value of lambda_surface (heat conductivity between air and wall) 3701 REAL(wp) :: Ueff !< effective wind speed for calculation of heat transfer coefficients 3702 REAL(wp) :: httc !< heat transfer coefficient 3703 REAL(wp), DIMENSION(nzub:nzut) :: exn !< value of the Exner function in layers 3704 3705 REAL(wp), DIMENSION(0:4) :: dxdir !< surface normal direction gridbox length 3706 REAL(wp) :: dtime !< simulated time of day (in UTC) 3707 INTEGER(iwp) :: dhour !< simulated hour of day (in UTC) 3708 REAL(wp) :: acoef !< actual coefficient of diurnal profile of anthropogenic heat 3709 3710 dxdir = (/dz,dy,dy,dx,dx/) 3711 3712 exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp !< Exner function 3713 3714 ! 3715 DO l = startenergy, endenergy 3716 ! Calculate frequently used parameters 3717 d = surfl(id,l) 3718 k = surfl(iz,l) 3719 j = surfl(iy,l) 3720 i = surfl(ix,l) 3721 3722 ! TODO  how to calculate lambda_surface for horizontal surfaces 3723 ! (lambda_surface is set according to stratification in land surface model) 3724 IF ( ol(j,i) >= 0.0_wp ) THEN 3725 lambda_surface = lambda_surf(l) 3726 ELSE 3727 lambda_surface = lambda_surf(l) 3159 3728 ENDIF 3160 3729 3161 #if defined( __parallel ) && ! defined ( __check ) 3162 CALL mpi_barrier( comm2d, ierr ) 3163 #endif 3164 ENDDO 3165 3166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3167 ! read diurnal profiles of heat sources 3168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3169 aheatprof = 0.0_wp 3170 DO ii = 0, io_blocks1 3171 IF ( ii == io_group ) THEN 3172 3173 ! open anthropogenic heat profile file 3174 OPEN( 151, file='ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char), action='read', & 3175 status='old', form='formatted', err=21 ) 3176 i = 0 3177 DO 3178 READ( 151, *, err=22, end=23 ) i, heat 3179 IF ( i >= 0 .AND. i <= 24 ) THEN 3180 ! write heat into the array 3181 aheatprof(i) = heat 3182 ENDIF 3183 CYCLE 3184 22 WRITE(message_string,'(a,i4)') 'error in file ANTHROPOGENIC_HEAT_PROFILE'// & 3185 TRIM(coupling_char)//' after line ',i 3186 CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 ) 3187 ENDDO 3188 aheatprof(24) = aheatprof(0) 3189 23 CLOSE(151) 3190 CYCLE 3191 21 message_string = 'file ANTHROPOGENIC_HEAT_PROFILE'//TRIM(coupling_char)//' does not exist' 3192 CALL message( 'usm_read_anthropogenic_heat', 'PA0518', 1, 2, 0, 6, 0 ) 3730 pt1 = pt(k,j,i) 3731 3732 ! calculate rho * cp coefficient at surface layer 3733 rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) ) 3734 3735 ! calculate aerodyamic resistance. 3736 IF ( d == iroof ) THEN 3737 ! calculation for horizontal surfaces follows LSM formulation 3738 ! pt, us, ts are not available for the prognostic time step, 3739 ! data from the last time step is used here. 3740 3741 r_a = (pt1  t_surf(l)/exn(k)) / (ts(j,i) * us(j,i) + 1.0E10_wp) 3742 3743 ! make sure that the resistance does not drop to zero 3744 IF ( ABS(r_a) < 1.0E10_wp ) r_a = 1.0E10_wp 3745 3746 ! the parameterization is developed originally for larger scales 3747 ! (compare with remark in TUF3D) 3748 ! our first experiences show that the parameterization underestimates 3749 ! r_a in meter resolution. 3750 ! temporary solution  multiplication by magic constant :(. 3751 r_a = r_a * ra_horiz_coef 3752 3753 ! factor for shf_eb 3754 f_shf = rho_cp / r_a 3755 ELSE 3756 ! calculation of r_a for vertical surfaces 3757 ! 3758 ! heat transfer coefficient for forced convection along vertical walls 3759 ! follows formulation in TUF3d model (Krayenhoff & Voogt, 2006) 3760 ! 3761 ! H = httc (Tsfc  Tair) 3762 ! httc = rw * (11.8 + 4.2 * Ueff)  4.0 3763 ! 3764 ! rw: wall patch roughness relative to 1.0 for concrete 3765 ! Ueff: effective wind speed 3766 !  4.0 is a reduction of Rowley et al (1930) formulation based on 3767 ! Cole and Sturrock (1977) 3768 ! 3769 ! Ucan: Canyon wind speed 3770 ! wstar: convective velocity 3771 ! Qs: surface heat flux 3772 ! zH: height of the convective layer 3773 ! wstar = (g/Tcan*Qs*zH)**(1./3.) 3774 3775 ! staggered grid needs to be taken into consideration 3776 IF ( d == inorth ) THEN