Changeset 2232 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/urban_surface_mod.f90 (modified) (87 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/urban_surface_mod.f90
r2214 r2232 21 21 ! Current revisions: 22 22 ! ------------------ 23 ! 23 ! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux; 24 ! insteat, heat fluxes are directly applied in diffusion_s. 24 25 ! 25 26 ! Former revisions: … … 103 104 !> reflection step using MPI_Alltoall instead of current MPI_Allgather. 104 105 !> 106 !> 3. Temporarily large values of surface heat flux can be observed, up to 107 !> 1.2 Km/s, which seem to be not realistic. 108 !> 105 109 !> @todo Check optimizations for RMA operations 106 110 !> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi) … … 111 115 112 116 USE arrays_3d, & 113 ONLY: zu, pt, pt_1, pt_2, p, ol, shf, ts, us,u, v, w, hyp, tend117 ONLY: zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend 114 118 115 119 USE cloud_parameters, & … … 137 141 USE indices, & 138 142 ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, & 139 nysg, nzb _s_inner, nzb_s_outer, nzb, nzt, nbgp143 nysg, nzb, nzt, nbgp, wall_flags_0 140 144 141 145 USE, INTRINSIC :: iso_c_binding … … 157 161 USE statistics, & 158 162 ONLY: hom, statistic_regions 163 164 USE surface_mod 159 165 160 166 … … 265 271 REAL(wp) :: rtransp !< 266 272 END TYPE 273 ! 274 !-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls. 275 TYPE t_surf_vertical 276 REAL(wp), DIMENSION(:), ALLOCATABLE :: t 277 END TYPE t_surf_vertical 278 ! 279 !-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls. 280 TYPE t_wall_vertical 281 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t 282 END TYPE t_wall_vertical 267 283 268 284 !-- arrays for calculation of svf and csf … … 320 336 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection 321 337 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection 322 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface 338 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface 323 339 324 340 !-- block variables needed for calculation of the plant canopy model inside the urban surface model … … 367 383 !-- surface and material model variables for walls, ground, roofs 368 384 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 369 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_types !< array of types of wall parameters370 385 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m) 371 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall !< 1/dz_wall372 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddz_wall_stag !< 1/dz_wall_stag373 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall !< wall grid spacing (center-center)374 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dz_wall_stag !< wall grid spacing (edge-edge)375 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw !< wall layer depths (m)376 386 377 387 #if defined( __nopointer ) 378 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf !< wall surface temperature (K) 379 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_p !< progn. wall surface temperature (K) 388 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h !< wall surface temperature (K) at horizontal walls 389 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_p !< progn. wall surface temperature (K) at horizontal walls 390 391 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v 392 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_p 380 393 #else 381 REAL(wp), DIMENSION(:), POINTER :: t_surf 382 REAL(wp), DIMENSION(:), POINTER :: t_surf_p 383 384 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_1 385 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_2 394 REAL(wp), DIMENSION(:), POINTER :: t_surf_h 395 REAL(wp), DIMENSION(:), POINTER :: t_surf_h_p 396 397 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_1 398 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_2 399 400 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v 401 TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v_p 402 403 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_1 404 TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_2 386 405 #endif 387 406 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K) … … 394 413 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 395 414 !-- parameters of the land, roof and wall surfaces 396 LOGICAL, DIMENSION(:), ALLOCATABLE :: isroof_surf !< is the surface the part of a roof397 415 REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface 398 416 !-- parameters of the wall surfaces 399 REAL(wp), DIMENSION(:), ALLOCATABLE :: c_surface !< heat capacity of the wall surface skin ( J mâ2 Kâ1 )400 417 REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface 401 REAL(wp), DIMENSION(:), ALLOCATABLE :: lambda_surf !< heat conductivity λS between air and surface ( W mâ2 Kâ1 )402 403 !-- parameters of the walls material404 REAL(wp), DIMENSION(:), ALLOCATABLE :: thickness_wall !< thickness of the wall, roof and soil layers405 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rho_c_wall !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)406 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lambda_h !< heat conductivity λT of the material ( W m-1 K-1 )407 REAL(wp), DIMENSION(:), ALLOCATABLE :: roughness_wall !< roughness relative to concrete408 409 !-- output wall heat flux arrays410 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf !< kinematic wall heat flux of sensible heat (needed for diffusion_s!<)411 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb !< wall heat flux of sensible heat in wall normal direction412 REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb413 REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb !< wall ground heat flux414 REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb415 418 416 419 #if defined( __nopointer ) 417 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall !< Wall temperature (K) 418 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_av !< Average of t_wall 419 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_p !< Prog. wall temperature (K) 420 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h !< Wall temperature (K) 421 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h_av !< Average of t_wall 422 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_wall_h_p !< Prog. wall temperature (K) 423 424 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v !< Wall temperature (K) 425 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av !< Average of t_wall 426 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_p !< Prog. wall temperature (K) 420 427 #else 421 REAL(wp), DIMENSION(:,:), POINTER :: t_wall, t_wall_p 422 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_av, t_wall_1, t_wall_2 428 REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p 429 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av, t_wall_h_1, t_wall_h_2 430 431 TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p 432 TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av, t_wall_v_1, t_wall_v_2 423 433 #endif 424 434 … … 488 498 MODULE PROCEDURE usm_swap_timelevel 489 499 END INTERFACE usm_swap_timelevel 490 491 INTERFACE usm_wall_heat_flux 492 MODULE PROCEDURE usm_wall_heat_flux 493 MODULE PROCEDURE usm_wall_heat_flux_ij 494 END INTERFACE usm_wall_heat_flux 495 500 496 501 INTERFACE usm_write_restart_data 497 502 MODULE PROCEDURE usm_write_restart_data … … 508 513 usm_energy_balance_land, usm_energy_balance_wall, nrefsteps, & 509 514 usm_init_urban_surface, usm_radiation, usm_read_restart_data, & 510 usm_wall_heat_flux, &511 515 usm_surface_energy_balance, usm_material_heat_model, & 512 516 usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, & … … 531 535 IMPLICIT NONE 532 536 533 INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids 534 INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb 535 INTEGER(iwp) :: procid 537 INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids, m 538 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 539 INTEGER(iwp) :: k_topo2 !< vertical index indicating topography top for given (j,i) 540 INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb 541 INTEGER(iwp) :: procid 536 542 537 543 … … 543 549 CALL location_message( '', .TRUE. ) 544 550 CALL location_message( ' allocation of needed arrays', .TRUE. ) 545 !-- find nzub, nzut, nzu 546 nzubl = minval(nzb_s_inner(nys:nyn,nxl:nxr)) 547 nzutl = maxval(nzb_s_inner(nys:nyn,nxl:nxr)) 551 ! 552 !-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be 553 !-- removed later). The following contruct finds the lowest / largest index 554 !-- for any upward-facing wall (see bit 12). 555 nzubl = MINVAL( & 556 MAXLOC( & 557 MERGE( 1, 0, & 558 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 559 ), DIM = 1 & 560 ) - 1 & 561 ) 562 nzutl = MAXVAL( & 563 MAXLOC( & 564 MERGE( 1, 0, & 565 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 566 ), DIM = 1 & 567 ) - 1 & 568 ) 548 569 nzubl = max(nzubl,nzb) 570 549 571 550 572 IF ( plant_canopy ) THEN … … 559 581 DO i = nxl, nxr 560 582 DO j = nys, nyn 583 ! 584 !-- Find topography top index 585 k_topo = MAXLOC( MERGE( 1, 0, & 586 BTEST( wall_flags_0(:,j,i), 12 ) & 587 ), DIM = 1 & 588 ) - 1 561 589 DO k = nzt+1, 0, -1 562 590 IF ( lad_s(k,j,i) /= 0.0_wp ) THEN 563 591 !-- we are at the top of the pcs 564 pct(j,i) = k + nzb_s_inner(j,i)592 pct(j,i) = k + k_topo 565 593 pch(j,i) = k 566 594 npcbl = npcbl + pch(j,i) … … 601 629 CALL location_message( ' calculation of indices for surfaces', .TRUE. ) 602 630 nsurfl = 0 603 !-- calculate land surface and roof 604 startland = nsurfl+1 605 nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1) 606 endland = nsurfl 607 nlands = endland-startland+1 608 609 !-- calculation of the walls 631 ! 632 !-- Number of land- and roof surfaces. Note, since horizontal surface elements 633 !-- are already counted in surface_mod, in case be simply reused here. 634 startland = 1 635 nsurfl = surf_usm_h%ns 636 endland = nsurfl 637 nlands = endland-startland+1 638 639 ! 640 !-- Number of vertical surfaces. As vertical surfaces are already 641 !-- counted in surface mod, it can be reused here. 610 642 startwall = nsurfl+1 611 DO i = nxl, nxr 612 DO j = nys, nyn 613 !-- test for walls 614 !-- (we don't use array flags because it isn't calculated in case of masking_method=.T.) 615 DO ids = 1, 4 !-- four wall directions 616 jr = min(max(j-jdir(ids),0),ny) 617 ir = min(max(i-idir(ids),0),nx) 618 nsurfl = nsurfl + max(0, nzb_s_inner(jr,ir)-nzb_s_inner(j,i)) 619 ENDDO 620 ENDDO 621 ENDDO 643 nsurfl = nsurfl + surf_usm_v(0)%ns + surf_usm_v(1)%ns + & 644 surf_usm_v(2)%ns + surf_usm_v(3)%ns 622 645 endwall = nsurfl 623 646 nwalls = endwall-startwall+1 624 625 !-- range of energy balance surfaces 647 648 649 !-- range of energy balance surfaces ! will be treated separately by surf_usm_h and surf_usm_v 626 650 nenergy = 0 627 651 IF ( usm_energy_balance_land ) THEN … … 641 665 !-- block of virtual surfaces 642 666 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 643 !-- calculate sky surfaces 667 !-- calculate sky surfaces ! not used so far! 644 668 startsky = nsurfl+1 645 669 nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1) … … 656 680 ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) ) 657 681 !-- calulation of the free borders of the domain 658 DO ids = 6,9 659 IF ( isborder(ids) ) THEN 660 !-- free border of the domain in direction ids 661 DO i = ijdb(1,ids), ijdb(2,ids) 662 DO j = ijdb(3,ids), ijdb(4,ids) 663 k = nzut - max(nzb_s_inner(j,i), nzb_s_inner(j-jdir(ids),i-idir(ids))) 664 nsurfl = nsurfl + k 665 ENDDO 666 ENDDO 667 ENDIF 682 DO ids = 6,9 683 IF ( isborder(ids) ) THEN 684 !-- free border of the domain in direction ids 685 DO i = ijdb(1,ids), ijdb(2,ids) 686 DO j = ijdb(3,ids), ijdb(4,ids) 687 688 k_topo = MAXLOC( MERGE( 1, 0, & 689 BTEST( wall_flags_0(:,j,i), 12 ) & 690 ), DIM = 1 & 691 ) - 1 692 k_topo2 = MAXLOC( MERGE( 1, 0, & 693 BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 ) & 694 ), DIM = 1 & 695 ) - 1 696 697 k = nzut - MAX( k_topo, k_topo2 ) 698 nsurfl = nsurfl + k 699 ENDDO 700 ENDDO 701 ENDIF 668 702 ENDDO 669 703 … … 676 710 DO i = nxl, nxr 677 711 DO j = nys, nyn 678 DO k = nzb_s_inner(j,i)+1, pct(j,i) 712 ! 713 !-- Find topography top index 714 k_topo = MAXLOC( MERGE( 1, 0, & 715 BTEST( wall_flags_0(:,j,i), 12 ) & 716 ), DIM = 1 & 717 ) - 1 718 DO k = k_topo + 1, pct(j,i) 679 719 ipcgb = ipcgb + 1 680 720 gridpcbl(k,j,i) = ipcgb … … 689 729 690 730 !-- fill surfl 691 ALLOCATE(surfl( 4,nsurfl))731 ALLOCATE(surfl(5,nsurfl)) 692 732 isurf = 0 693 733 … … 695 735 DO i = nxl, nxr 696 736 DO j = nys, nyn 697 isurf = isurf + 1 698 k = nzb_s_inner(j,i)+1 699 surfl(:,isurf) = (/iroof,k,j,i/) 737 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 738 k = surf_usm_h%k(m) 739 740 isurf = isurf + 1 741 surfl(:,isurf) = (/iroof,k,j,i,m/) 742 ENDDO 700 743 ENDDO 701 744 ENDDO … … 704 747 DO i = nxl, nxr 705 748 DO j = nys, nyn 706 DO ids = 1, 4 !> four wall directions 707 jr = min(max(j-jdir(ids),0),ny) 708 ir = min(max(i-idir(ids),0),nx) 709 DO k = nzb_s_inner(j,i)+1, nzb_s_inner(jr,ir) 710 isurf = isurf + 1 711 surfl(:,isurf) = (/ids,k,j,i/) 712 ENDDO 713 ENDDO 749 l = 0 750 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 751 k = surf_usm_v(l)%k(m) 752 753 isurf = isurf + 1 754 surfl(:,isurf) = (/2,k,j,i,m/) 755 ENDDO 756 l = 1 757 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 758 k = surf_usm_v(l)%k(m) 759 760 isurf = isurf + 1 761 surfl(:,isurf) = (/1,k,j,i,m/) 762 ENDDO 763 l = 2 764 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 765 k = surf_usm_v(l)%k(m) 766 767 isurf = isurf + 1 768 surfl(:,isurf) = (/4,k,j,i,m/) 769 ENDDO 770 l = 3 771 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 772 k = surf_usm_v(l)%k(m) 773 774 isurf = isurf + 1 775 surfl(:,isurf) = (/3,k,j,i,m/) 776 ENDDO 714 777 ENDDO 715 778 ENDDO … … 720 783 isurf = isurf + 1 721 784 k = nzut 722 surfl(:,isurf) = (/isky,k,j,i /)785 surfl(:,isurf) = (/isky,k,j,i,-1/) 723 786 ENDDO 724 787 ENDDO … … 730 793 DO i = ijdb(1,ids), ijdb(2,ids) 731 794 DO j = ijdb(3,ids), ijdb(4,ids) 732 DO k = max(nzb_s_inner(j,i),nzb_s_inner(j-jdir(ids),i-idir(ids)))+1, nzut 795 k_topo = MAXLOC( MERGE( 1, 0, & 796 BTEST( wall_flags_0(:,j,i), 12 ) & 797 ), DIM = 1 & 798 ) - 1 799 k_topo2 = MAXLOC( MERGE( 1, 0, & 800 BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 ) & 801 ), DIM = 1 & 802 ) - 1 803 DO k = MAX(k_topo,k_topo2)+1, nzut 733 804 isurf = isurf + 1 734 surfl(:,isurf) = (/ids,k,j,i /)805 surfl(:,isurf) = (/ids,k,j,i,-1/) 735 806 ENDDO 736 807 ENDDO … … 755 826 surfstart(numprocs) = k 756 827 nsurf = k 757 ALLOCATE(surf( 4,nsurf))828 ALLOCATE(surf(5,nsurf)) 758 829 759 830 #if defined( __parallel ) 760 CALL MPI_AllGatherv(surfl, nsurfl* 4, MPI_INTEGER, surf, nsurfs*4, surfstart*4, MPI_INTEGER, comm2d, ierr)831 CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, surfstart*5, MPI_INTEGER, comm2d, ierr) 761 832 #else 762 833 surf = surfl … … 787 858 ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual 788 859 ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual 789 ALLOCATE( surfhf(startenergy:endenergy) ) 790 ALLOCATE( rad_net_l(startenergy:endenergy) ) 860 861 862 863 ! 864 !-- Allocate radiation arrays which are part of the new data type. 865 !-- For horizontal surfaces. 866 ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns) ) 867 ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) ) 868 ! 869 !-- New 870 ALLOCATE( surf_usm_h%rad_in_sw(1:surf_usm_h%ns) ) 871 ALLOCATE( surf_usm_h%rad_out_sw(1:surf_usm_h%ns) ) 872 ALLOCATE( surf_usm_h%rad_in_lw(1:surf_usm_h%ns) ) 873 ALLOCATE( surf_usm_h%rad_out_lw(1:surf_usm_h%ns) ) 874 ! 875 !-- For vertical surfaces 876 DO l = 0, 3 877 ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) ) 878 ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) ) 879 ALLOCATE( surf_usm_v(l)%rad_in_sw(1:surf_usm_v(l)%ns) ) 880 ALLOCATE( surf_usm_v(l)%rad_out_sw(1:surf_usm_v(l)%ns) ) 881 ALLOCATE( surf_usm_v(l)%rad_in_lw(1:surf_usm_v(l)%ns) ) 882 ALLOCATE( surf_usm_v(l)%rad_out_lw(1:surf_usm_v(l)%ns) ) 883 ENDDO 791 884 792 885 !-- Wall surface model … … 794 887 795 888 !-- allocate array of wall types and wall parameters 796 ALLOCATE ( surface_types(startenergy:endenergy) ) 889 ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) ) 890 DO l = 0, 3 891 ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) ) 892 ENDDO 797 893 798 894 !-- broadband albedo of the land, roof and wall surface … … 802 898 ALLOCATE ( albedo_surf(nsurfl) ) 803 899 albedo_surf = 1.0_wp 900 ALLOCATE ( surf_usm_h%albedo_surf(1:surf_usm_h%ns) ) 901 DO l = 0, 3 902 ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) ) 903 ENDDO 804 904 805 !-- wall and roof surface parameters 806 ALLOCATE ( isroof_surf(startenergy:endenergy) ) 905 !-- wall and roof surface parameters. First for horizontal surfaces 807 906 ALLOCATE ( emiss_surf(startenergy:endenergy) ) 808 ALLOCATE ( lambda_surf(startenergy:endenergy) ) 809 ALLOCATE ( c_surface(startenergy:endenergy) ) 810 ALLOCATE ( roughness_wall(startenergy:endenergy) ) 907 908 ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) ) 909 ALLOCATE ( surf_usm_h%emiss_surf(1:surf_usm_h%ns) ) 910 ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) ) 911 ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) ) 912 ALLOCATE ( surf_usm_h%roughness_wall(1:surf_usm_h%ns) ) 913 ! 914 !-- For vertical surfaces. 915 DO l = 0, 3 916 ALLOCATE ( surf_usm_v(l)%emiss_surf(1:surf_usm_v(l)%ns) ) 917 ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) ) 918 ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) ) 919 ALLOCATE ( surf_usm_v(l)%roughness_wall(1:surf_usm_v(l)%ns) ) 920 ENDDO 811 921 812 !-- allocate wall and roof material parameters 813 ALLOCATE ( thickness_wall(startenergy:endenergy) ) 814 ALLOCATE ( lambda_h(nzb_wall:nzt_wall,startenergy:endenergy) ) 815 ALLOCATE ( rho_c_wall(nzb_wall:nzt_wall,startenergy:endenergy) ) 816 817 !-- allocate wall and roof layers sizes 922 !-- allocate wall and roof material parameters. First for horizontal surfaces 923 ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) ) 924 ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 925 ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 926 ! 927 !-- For vertical surfaces. 928 DO l = 0, 3 929 ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) ) 930 ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 931 ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 932 ENDDO 933 934 !-- allocate wall and roof layers sizes. For horizontal surfaces. 818 935 ALLOCATE ( zwn(nzb_wall:nzt_wall) ) 819 ALLOCATE ( dz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) ) 820 ALLOCATE ( ddz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) ) 821 ALLOCATE ( dz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) ) 822 ALLOCATE ( ddz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) ) 823 ALLOCATE ( zw(nzb_wall:nzt_wall, startenergy:endenergy) ) 824 825 !-- allocate wall and roof temperature arrays 936 ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 937 ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 938 ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 939 ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 940 ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 941 ! 942 !-- For vertical surfaces. 943 DO l = 0, 3 944 ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 945 ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 946 ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 947 ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 948 ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 949 ENDDO 950 951 !-- allocate wall and roof temperature arrays, for horizontal walls 826 952 #if defined( __nopointer ) 827 ALLOCATE ( t_surf(startenergy:endenergy) ) 828 ALLOCATE ( t_surf_p(startenergy:endenergy) ) 829 ALLOCATE ( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 830 ALLOCATE ( t_wall_p(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 953 ALLOCATE ( t_surf_h(1:surf_usm_h%ns) ) 954 ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) ) 955 ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 956 ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 957 958 ALLOCATE ( t_surf_h(1:surf_usm_h%ns) ) 959 ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) ) 960 ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 961 ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 831 962 #else 832 ALLOCATE ( t_surf_ 1(startenergy:endenergy) )833 ALLOCATE ( t_surf_ 2(startenergy:endenergy) )834 ALLOCATE ( t_wall_ 1(nzb_wall:nzt_wall+1,startenergy:endenergy) )835 ALLOCATE ( t_wall_ 2(nzb_wall:nzt_wall+1,startenergy:endenergy) )963 ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) ) 964 ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) ) 965 ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 966 ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 836 967 837 968 !-- initial assignment of the pointers 838 t_wall => t_wall_1; t_wall_p => t_wall_2839 t_surf => t_surf_1; t_surf_p => t_surf_2969 t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2 970 t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2 840 971 #endif 841 972 842 !-- allocate intermediate timestep arrays 843 ALLOCATE ( tt_surface_m(startenergy:endenergy) ) 844 ALLOCATE ( tt_wall_m(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 845 846 !-- allocate wall heat flux output array 847 ALLOCATE ( wshf(startwall:endwall) ) 848 ALLOCATE ( wshf_eb(startenergy:endenergy) ) 849 ALLOCATE ( wghf_eb(startenergy:endenergy) ) 850 851 !-- set inital values for prognostic quantities 852 tt_surface_m = 0.0_wp 853 tt_wall_m = 0.0_wp 854 855 wshf = 0.0_wp 856 wshf_eb = 0.0_wp 857 wghf_eb = 0.0_wp 973 !-- allocate wall and roof temperature arrays, for vertical walls 974 #if defined( __nopointer ) 975 DO l = 0, 3 976 ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) ) 977 ALLOCATE ( t_surf_v(l)%t_p(1:surf_usm_v(l)%ns) ) 978 ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 979 ALLOCATE ( t_wall_v(l)%t_p(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 980 ENDDO 981 #else 982 DO l = 0, 3 983 ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) ) 984 ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) ) 985 ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 986 ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 987 ENDDO 988 ! 989 !-- initial assignment of the pointers 990 t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2 991 t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2 992 #endif 993 ! 994 !-- Allocate intermediate timestep arrays. For horizontal surfaces. 995 ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns) ) 996 ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 997 ! 998 !-- Set inital values for prognostic quantities 999 IF ( ALLOCATED( surf_usm_h%tt_surface_m ) ) surf_usm_h%tt_surface_m = 0.0_wp 1000 IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp 1001 ! 1002 !-- Now, for vertical surfaces 1003 DO l = 0, 3 1004 ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns) ) 1005 ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 1006 IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) ) surf_usm_v(l)%tt_surface_m = 0.0_wp 1007 IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp 1008 ENDDO 1009 1010 !-- allocate wall heat flux output array and set initial values. For horizontal surfaces 1011 ! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed 1012 ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) ) 1013 ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) ) 1014 IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp 1015 IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp 1016 IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp 1017 ! 1018 !-- Now, for vertical surfaces 1019 DO l = 0, 3 1020 ! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed 1021 ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) ) 1022 ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) ) 1023 IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp 1024 IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp 1025 IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp 1026 ENDDO 858 1027 859 1028 END SUBROUTINE usm_allocate_urban_surface … … 874 1043 CHARACTER (len=*), INTENT(IN) :: variable 875 1044 876 INTEGER(iwp) :: i, j, k, l, ids, iwl,istat1045 INTEGER(iwp) :: i, j, k, l, m, ids, iwl,istat 877 1046 CHARACTER (len=varnamelength) :: var, surfid 878 1047 INTEGER(iwp), PARAMETER :: nd = 5 … … 910 1079 CASE ( 'usm_rad_net' ) 911 1080 !-- array of complete radiation balance 912 IF ( .NOT. ALLOCATED( rad_net_av) ) THEN913 ALLOCATE( rad_net_av(startenergy:endenergy) )914 rad_net_av = 0.0_wp1081 IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN 1082 ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) ) 1083 surf_usm_h%rad_net_av = 0.0_wp 915 1084 ENDIF 1085 DO l = 0, 3 1086 IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN 1087 ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) ) 1088 surf_usm_v(l)%rad_net_av = 0.0_wp 1089 ENDIF 1090 ENDDO 916 1091 917 1092 CASE ( 'usm_rad_insw' ) 918 1093 !-- array of sw radiation falling to surface after i-th reflection 919 IF ( .NOT. ALLOCATED(surf insw_av) ) THEN920 ALLOCATE( surf insw_av(startenergy:endenergy) )921 surf insw_av = 0.0_wp1094 IF ( .NOT. ALLOCATED(surf_usm_h%surfinsw_av) ) THEN 1095 ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) ) 1096 surf_usm_h%surfinsw_av = 0.0_wp 922 1097 ENDIF 1098 DO l = 0, 3 1099 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinsw_av) ) THEN 1100 ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) ) 1101 surf_usm_v(l)%surfinsw_av = 0.0_wp 1102 ENDIF 1103 ENDDO 923 1104 924 1105 CASE ( 'usm_rad_inlw' ) 925 1106 !-- array of lw radiation falling to surface after i-th reflection 926 IF ( .NOT. ALLOCATED(surf inlw_av) ) THEN927 ALLOCATE( surf inlw_av(startenergy:endenergy) )928 surf inlw_av = 0.0_wp1107 IF ( .NOT. ALLOCATED(surf_usm_h%surfinlw_av) ) THEN 1108 ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) ) 1109 surf_usm_h%surfinlw_av = 0.0_wp 929 1110 ENDIF 1111 DO l = 0, 3 1112 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinlw_av) ) THEN 1113 ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) ) 1114 surf_usm_v(l)%surfinlw_av = 0.0_wp 1115 ENDIF 1116 ENDDO 930 1117 931 1118 CASE ( 'usm_rad_inswdir' ) … … 952 1139 CASE ( 'usm_rad_inlwdif' ) 953 1140 !-- array of sw radiation falling to surface after i-th reflection 954 1141 IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN 955 1142 ALLOCATE( surfinlwdif_av(startenergy:endenergy) ) 956 1143 surfinlwdif_av = 0.0_wp … … 977 1164 surfoutlw_av = 0.0_wp 978 1165 ENDIF 979 980 1166 CASE ( 'usm_rad_ressw' ) 981 1167 !-- array of residua of sw radiation absorbed in surface after last reflection … … 984 1170 surfins_av = 0.0_wp 985 1171 ENDIF 986 1172 987 1173 CASE ( 'usm_rad_reslw' ) 988 1174 !-- array of residua of lw radiation absorbed in surface after last reflection … … 994 1180 CASE ( 'usm_rad_hf' ) 995 1181 !-- array of heat flux from radiation for surfaces after i-th reflection 996 IF ( .NOT. ALLOCATED(surf hf_av) ) THEN997 ALLOCATE( surf hf_av(startenergy:endenergy) )998 surf hf_av = 0.0_wp1182 IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN 1183 ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) ) 1184 surf_usm_h%surfhf_av = 0.0_wp 999 1185 ENDIF 1186 DO l = 0, 3 1187 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN 1188 ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) ) 1189 surf_usm_v(l)%surfhf_av = 0.0_wp 1190 ENDIF 1191 ENDDO 1000 1192 1001 1193 CASE ( 'usm_wshf' ) 1002 1194 !-- array of sensible heat flux from surfaces 1003 1195 !-- land surfaces 1004 IF ( .NOT. ALLOCATED( wshf_eb_av) ) THEN1005 ALLOCATE( wshf_eb_av(startenergy:endenergy) )1006 wshf_eb_av = 0.0_wp1196 IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN 1197 ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) ) 1198 surf_usm_h%wshf_eb_av = 0.0_wp 1007 1199 ENDIF 1200 DO l = 0, 3 1201 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN 1202 ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) ) 1203 surf_usm_v(l)%wshf_eb_av = 0.0_wp 1204 ENDIF 1205 ENDDO 1008 1206 1009 1207 CASE ( 'usm_wghf' ) 1010 1208 !-- array of heat flux from ground (wall, roof, land) 1011 IF ( .NOT. ALLOCATED( wghf_eb_av) ) THEN1012 ALLOCATE( wghf_eb_av(startenergy:endenergy) )1013 wghf_eb_av = 0.0_wp1209 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN 1210 ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) ) 1211 surf_usm_h%wghf_eb_av = 0.0_wp 1014 1212 ENDIF 1213 DO l = 0, 3 1214 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN 1215 ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) ) 1216 surf_usm_v(l)%wghf_eb_av = 0.0_wp 1217 ENDIF 1218 ENDDO 1015 1219 1016 1220 CASE ( 'usm_t_surf' ) 1017 1221 !-- surface temperature for surfaces 1018 IF ( .NOT. ALLOCATED( t_surf_av) ) THEN1019 ALLOCATE( t_surf_av(startenergy:endenergy) )1020 t_surf_av = 0.0_wp1222 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN 1223 ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) ) 1224 surf_usm_h%t_surf_av = 0.0_wp 1021 1225 ENDIF 1226 DO l = 0, 3 1227 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_av) ) THEN 1228 ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) ) 1229 surf_usm_v(l)%t_surf_av = 0.0_wp 1230 ENDIF 1231 ENDDO 1022 1232 1023 1233 CASE ( 'usm_t_wall' ) 1024 1234 !-- wall temperature for iwl layer of walls and land 1025 IF ( .NOT. ALLOCATED( t_wall_av) ) THEN1026 ALLOCATE( t_wall_av(nzb_wall:nzt_wall,startenergy:endenergy) )1027 t_wall_av = 0.0_wp1235 IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN 1236 ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 1237 surf_usm_h%t_wall_av = 0.0_wp 1028 1238 ENDIF 1239 DO l = 0, 3 1240 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN 1241 ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) ) 1242 surf_usm_v(l)%t_wall_av = 0.0_wp 1243 ENDIF 1244 ENDDO 1029 1245 1030 1246 CASE DEFAULT … … 1039 1255 CASE ( 'usm_rad_net' ) 1040 1256 !-- array of complete radiation balance 1041 DO l = startenergy, endenergy 1042 IF ( surfl(id,l) == ids ) THEN 1043 rad_net_av(l) = rad_net_av(l) + rad_net_l(l) 1044 ENDIF 1257 DO m = 1, surf_usm_h%ns 1258 surf_usm_h%rad_net_av(m) = & 1259 surf_usm_h%rad_net_av(m) + & 1260 surf_usm_h%rad_net_l(m) 1261 ENDDO 1262 DO l = 0, 3 1263 DO m = 1, surf_usm_v(l)%ns 1264 surf_usm_v(l)%rad_net_av(m) = & 1265 surf_usm_v(l)%rad_net_av(m) + & 1266 surf_usm_v(l)%rad_net_l(m) 1267 ENDDO 1045 1268 ENDDO 1046 1269 … … 1085 1308 ENDIF 1086 1309 ENDDO 1310 1087 1311 1088 1312 CASE ( 'usm_rad_inlwdif' ) … … 1090 1314 DO l = startenergy, endenergy 1091 1315 IF ( surfl(id,l) == ids ) THEN 1316 surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - & 1317 surfinswdir(l) - surfinswdif(l) 1318 ENDIF 1319 ENDDO 1320 ! 1321 CASE ( 'usm_rad_inlwref' ) 1322 !-- array of lw radiation falling to surface from reflections 1323 DO l = startenergy, endenergy 1324 IF ( surfl(id,l) == ids ) THEN 1092 1325 surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l) 1093 1326 ENDIF 1094 1327 ENDDO 1095 1328 1096 CASE ( 'usm_rad_ inlwref' )1097 !-- array of lw radiation falling to surface from reflections1329 CASE ( 'usm_rad_outsw' ) 1330 !-- array of sw radiation emitted from surface after i-th reflection 1098 1331 DO l = startenergy, endenergy 1099 1332 IF ( surfl(id,l) == ids ) THEN … … 1103 1336 ENDDO 1104 1337 1105 CASE ( 'usm_rad_outsw' )1106 !-- array of sw radiation emitted from surface after i-th reflection1107 DO l = startenergy, endenergy1108 IF ( surfl(id,l) == ids ) THEN1109 surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)1110 ENDIF1111 ENDDO1112 1113 1338 CASE ( 'usm_rad_outlw' ) 1114 1339 !-- array of lw radiation emitted from surface after i-th reflection 1115 1340 DO l = startenergy, endenergy 1116 1341 IF ( surfl(id,l) == ids ) THEN 1117 surfout lw_av(l) = surfoutlw_av(l) + surfoutlw(l)1342 surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l) 1118 1343 ENDIF 1119 1344 ENDDO … … 1123 1348 DO l = startenergy, endenergy 1124 1349 IF ( surfl(id,l) == ids ) THEN 1125 surf ins_av(l) = surfins_av(l) + surfins(l)1350 surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l) 1126 1351 ENDIF 1127 1352 ENDDO … … 1131 1356 DO l = startenergy, endenergy 1132 1357 IF ( surfl(id,l) == ids ) THEN 1133 surfin l_av(l) = surfinl_av(l) + surfinl(l)1358 surfins_av(l) = surfins_av(l) + surfins(l) 1134 1359 ENDIF 1135 1360 ENDDO … … 1137 1362 CASE ( 'usm_rad_hf' ) 1138 1363 !-- array of heat flux from radiation for surfaces after i-th reflection 1139 DO l = startenergy, endenergy 1140 IF ( surfl(id,l) == ids ) THEN 1141 surfhf_av(l) = surfhf_av(l) + surfhf(l) 1142 ENDIF 1364 DO m = 1, surf_usm_h%ns 1365 surf_usm_h%surfhf_av(m) = & 1366 surf_usm_h%surfhf_av(m) + & 1367 surf_usm_h%surfhf(m) 1368 ENDDO 1369 DO l = 0, 3 1370 DO m = 1, surf_usm_v(l)%ns 1371 surf_usm_v(l)%surfhf_av(m) = & 1372 surf_usm_v(l)%surfhf_av(m) + & 1373 surf_usm_v(l)%surfhf(m) 1374 ENDDO 1143 1375 ENDDO 1144 1376 1145 1377 CASE ( 'usm_wshf' ) 1146 1378 !-- array of sensible heat flux from surfaces (land, roof, wall) 1147 DO l = startenergy, endenergy 1148 IF ( surfl(id,l) == ids ) THEN 1149 wshf_eb_av(l) = wshf_eb_av(l) + wshf_eb(l) 1150 ENDIF 1379 DO m = 1, surf_usm_h%ns 1380 surf_usm_h%wshf_eb_av(m) = & 1381 surf_usm_h%wshf_eb_av(m) + & 1382 surf_usm_h%wshf_eb(m) 1383 ENDDO 1384 DO l = 0, 3 1385 DO m = 1, surf_usm_v(l)%ns 1386 surf_usm_v(l)%wshf_eb_av(m) = & 1387 surf_usm_v(l)%wshf_eb_av(m) + & 1388 surf_usm_v(l)%wshf_eb(m) 1389 ENDDO 1151 1390 ENDDO 1152 1391 1153 1392 CASE ( 'usm_wghf' ) 1154 1393 !-- array of heat flux from ground (wall, roof, land) 1155 DO l = startenergy, endenergy 1156 IF ( surfl(id,l) == ids ) THEN 1157 wghf_eb_av(l) = wghf_eb_av(l) + wghf_eb(l) 1158 ENDIF 1394 DO m = 1, surf_usm_h%ns 1395 surf_usm_h%wghf_eb_av(m) = & 1396 surf_usm_h%wghf_eb_av(m) + & 1397 surf_usm_h%wghf_eb(m) 1398 ENDDO 1399 DO l = 0, 3 1400 DO m = 1, surf_usm_v(l)%ns 1401 surf_usm_v(l)%wghf_eb_av(m) = & 1402 surf_usm_v(l)%wghf_eb_av(m) + & 1403 surf_usm_v(l)%wghf_eb(m) 1404 ENDDO 1159 1405 ENDDO 1160 1406 1161 1407 CASE ( 'usm_t_surf' ) 1162 1408 !-- surface temperature for surfaces 1163 DO l = startenergy, endenergy 1164 IF ( surfl(id,l) == ids ) THEN 1165 t_surf_av(l) = t_surf_av(l) + t_surf(l) 1166 ENDIF 1409 DO m = 1, surf_usm_h%ns 1410 surf_usm_h%t_surf_av(m) = & 1411 surf_usm_h%t_surf_av(m) + & 1412 t_surf_h(m) 1413 ENDDO 1414 DO l = 0, 3 1415 DO m = 1, surf_usm_v(l)%ns 1416 surf_usm_v(l)%t_surf_av(m) = & 1417 surf_usm_v(l)%t_surf_av(m) + & 1418 t_surf_v(l)%t(m) 1419 ENDDO 1167 1420 ENDDO 1168 1421 1169 1422 CASE ( 'usm_t_wall' ) 1170 1423 !-- wall temperature for iwl layer of walls and land 1171 DO l = startenergy, endenergy 1172 IF ( surfl(id,l) == ids ) THEN 1173 t_wall_av(iwl, l) = t_wall_av(iwl,l) + t_wall(iwl, l) 1174 ENDIF 1424 DO m = 1, surf_usm_h%ns 1425 surf_usm_h%t_wall_av(iwl,m) = & 1426 surf_usm_h%t_wall_av(iwl,m) + & 1427 t_wall_h(iwl,m) 1428 ENDDO 1429 DO l = 0, 3 1430 DO m = 1, surf_usm_v(l)%ns 1431 surf_usm_v(l)%t_wall_av(iwl,m) = & 1432 surf_usm_v(l)%t_wall_av(iwl,m) + & 1433 t_wall_v(l)%t(iwl,m) 1434 ENDDO 1175 1435 ENDDO 1176 1436 … … 1186 1446 CASE ( 'usm_rad_net' ) 1187 1447 !-- array of complete radiation balance 1188 DO l = startenergy, endenergy 1189 IF ( surfl(id,l) == ids ) THEN 1190 rad_net_av(l) = rad_net_av(l) / REAL( average_count_3d, kind=wp ) 1191 ENDIF 1448 DO m = 1, surf_usm_h%ns 1449 surf_usm_h%rad_net_av(m) = & 1450 surf_usm_h%rad_net_av(m) / & 1451 REAL( average_count_3d, kind=wp ) 1452 ENDDO 1453 DO l = 0, 3 1454 DO m = 1, surf_usm_v(l)%ns 1455 surf_usm_v(l)%rad_net_av(m) = & 1456 surf_usm_v(l)%rad_net_av(m) / & 1457 REAL( average_count_3d, kind=wp ) 1458 ENDDO 1192 1459 ENDDO 1193 1460 … … 1199 1466 ENDIF 1200 1467 ENDDO 1201 1468 1202 1469 CASE ( 'usm_rad_inlw' ) 1203 1470 !-- array of lw radiation falling to surface after i-th reflection … … 1207 1474 ENDIF 1208 1475 ENDDO 1209 1476 1210 1477 CASE ( 'usm_rad_inswdir' ) 1211 1478 !-- array of direct sw radiation falling to surface from sun … … 1215 1482 ENDIF 1216 1483 ENDDO 1217 1484 1218 1485 CASE ( 'usm_rad_inswdif' ) 1219 1486 !-- array of difusion sw radiation falling to surface from sky and borders of the domain … … 1223 1490 ENDIF 1224 1491 ENDDO 1225 1492 1226 1493 CASE ( 'usm_rad_inswref' ) 1227 1494 !-- array of sw radiation falling to surface from reflections … … 1231 1498 ENDIF 1232 1499 ENDDO 1233 1500 1234 1501 CASE ( 'usm_rad_inlwdif' ) 1235 1502 !-- array of sw radiation falling to surface after i-th reflection … … 1239 1506 ENDIF 1240 1507 ENDDO 1241 1508 1242 1509 CASE ( 'usm_rad_inlwref' ) 1243 1510 !-- array of lw radiation falling to surface from reflections … … 1247 1514 ENDIF 1248 1515 ENDDO 1249 1516 1250 1517 CASE ( 'usm_rad_outsw' ) 1251 1518 !-- array of sw radiation emitted from surface after i-th reflection … … 1255 1522 ENDIF 1256 1523 ENDDO 1257 1524 1258 1525 CASE ( 'usm_rad_outlw' ) 1259 1526 !-- array of lw radiation emitted from surface after i-th reflection … … 1263 1530 ENDIF 1264 1531 ENDDO 1265 1532 1266 1533 CASE ( 'usm_rad_ressw' ) 1267 1534 !-- array of residua of sw radiation absorbed in surface after last reflection … … 1282 1549 CASE ( 'usm_rad_hf' ) 1283 1550 !-- array of heat flux from radiation for surfaces after i-th reflection 1284 DO l = startenergy, endenergy 1285 IF ( surfl(id,l) == ids ) THEN 1286 surfhf_av(l) = surfhf_av(l) / REAL( average_count_3d, kind=wp ) 1287 ENDIF 1288 ENDDO 1289 1551 DO m = 1, surf_usm_h%ns 1552 surf_usm_h%surfhf_av(m) = & 1553 surf_usm_h%surfhf_av(m) / & 1554 REAL( average_count_3d, kind=wp ) 1555 ENDDO 1556 DO l = 0, 3 1557 DO m = 1, surf_usm_v(l)%ns 1558 surf_usm_v(l)%surfhf_av(m) = & 1559 surf_usm_v(l)%surfhf_av(m) / & 1560 REAL( average_count_3d, kind=wp ) 1561 ENDDO 1562 ENDDO 1563 1290 1564 CASE ( 'usm_wshf' ) 1291 1565 !-- array of sensible heat flux from surfaces (land, roof, wall) 1292 DO l = startenergy, endenergy 1293 IF ( surfl(id,l) == ids ) THEN 1294 wshf_eb_av(l) = wshf_eb_av(l) / REAL( average_count_3d, kind=wp ) 1295 ENDIF 1296 ENDDO 1297 1566 DO m = 1, surf_usm_h%ns 1567 surf_usm_h%wshf_eb_av(m) = & 1568 surf_usm_h%wshf_eb_av(m) / & 1569 REAL( average_count_3d, kind=wp ) 1570 ENDDO 1571 DO l = 0, 3 1572 DO m = 1, surf_usm_v(l)%ns 1573 surf_usm_v(l)%wshf_eb_av(m) = & 1574 surf_usm_v(l)%wshf_eb_av(m) / & 1575 REAL( average_count_3d, kind=wp ) 1576 ENDDO 1577 ENDDO 1578 1298 1579 CASE ( 'usm_wghf' ) 1299 1580 !-- array of heat flux from ground (wall, roof, land) 1300 DO l = startenergy, endenergy 1301 IF ( surfl(id,l) == ids ) THEN 1302 wghf_eb_av(l) = wghf_eb_av(l) / REAL( average_count_3d, kind=wp ) 1303 ENDIF 1304 ENDDO 1305 1581 DO m = 1, surf_usm_h%ns 1582 surf_usm_h%wghf_eb_av(m) = & 1583 surf_usm_h%wghf_eb_av(m) / & 1584 REAL( average_count_3d, kind=wp ) 1585 ENDDO 1586 DO l = 0, 3 1587 DO m = 1, surf_usm_v(l)%ns 1588 surf_usm_v(l)%wghf_eb_av(m) = & 1589 surf_usm_v(l)%wghf_eb_av(m) / & 1590 REAL( average_count_3d, kind=wp ) 1591 ENDDO 1592 ENDDO 1593 1306 1594 CASE ( 'usm_t_surf' ) 1307 1595 !-- surface temperature for surfaces 1308 DO l = startenergy, endenergy 1309 IF ( surfl(id,l) == ids ) THEN 1310 t_surf_av(l) = t_surf_av(l) / REAL( average_count_3d, kind=wp ) 1311 ENDIF 1312 ENDDO 1313 1596 DO m = 1, surf_usm_h%ns 1597 surf_usm_h%t_surf_av(m) = & 1598 surf_usm_h%t_surf_av(m) / & 1599 REAL( average_count_3d, kind=wp ) 1600 ENDDO 1601 DO l = 0, 3 1602 DO m = 1, surf_usm_v(l)%ns 1603 surf_usm_v(l)%t_surf_av(m) = & 1604 surf_usm_v(l)%t_surf_av(m) / & 1605 REAL( average_count_3d, kind=wp ) 1606 ENDDO 1607 ENDDO 1608 1314 1609 CASE ( 'usm_t_wall' ) 1315 1610 !-- wall temperature for iwl layer of walls and land 1316 DO l = startenergy, endenergy 1317 IF ( surfl(id,l) == ids ) THEN 1318 t_wall_av(iwl, l) = t_wall_av(iwl,l) / REAL( average_count_3d, kind=wp ) 1319 ENDIF 1611 DO m = 1, surf_usm_h%ns 1612 surf_usm_h%t_wall_av(iwl,m) = & 1613 surf_usm_h%t_wall_av(iwl,m) / & 1614 REAL( average_count_3d, kind=wp ) 1615 ENDDO 1616 DO l = 0, 3 1617 DO m = 1, surf_usm_v(l)%ns 1618 surf_usm_v(l)%t_wall_av(iwl,m) = & 1619 surf_usm_v(l)%t_wall_av(iwl,m) / & 1620 REAL( average_count_3d, kind=wp ) 1621 ENDDO 1320 1622 ENDDO 1321 1623 … … 1493 1795 TYPE(c_ptr) :: lad_s_rma_p !< allocated c pointer 1494 1796 INTEGER(kind=MPI_ADDRESS_KIND) :: size_lad_rma 1495 1797 ! 1496 1798 !-- calculation of the SVF 1497 1799 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 1498 1499 1800 ! 1500 1801 !-- precalculate face areas for different face directions using normal vector 1501 1802 DO d = 0, 9 … … 1524 1825 #if defined( __parallel ) 1525 1826 ALLOCATE( nzterrl(nys:nyn,nxl:nxr) ) 1526 nzterrl = nzb_s_inner(nys:nyn,nxl:nxr) 1827 nzterrl = MAXLOC( & 1828 MERGE( 1, 0, & 1829 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 1830 ), DIM = 1 & 1831 ) - 1 ! = nzb_s_inner(nys:nyn,nxl:nxr) 1527 1832 CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, & 1528 1833 nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr ) 1529 1834 DEALLOCATE(nzterrl) 1530 1835 #else 1531 nzterr = RESHAPE( nzb_s_inner(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) ) 1836 nzterr = RESHAPE( MAXLOC( & 1837 MERGE( 1, 0, & 1838 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 ) & 1839 ), DIM = 1 & 1840 ) - 1, & 1841 (/(nx+1)*(ny+1)/) & 1842 ) 1532 1843 #endif 1533 1844 IF ( plant_canopy ) THEN … … 1579 1890 DO i = nxl, nxr 1580 1891 DO j = nys, nyn 1581 k = nzb_s_inner(j, i) 1892 k = MAXLOC( & 1893 MERGE( 1, 0, & 1894 BTEST( wall_flags_0(:,j,i), 12 ) & 1895 ), DIM = 1 & 1896 ) - 1 1897 1582 1898 usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i) 1583 1899 ENDDO … … 1972 2288 CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 ) 1973 2289 1974 1975 2290 END SUBROUTINE usm_calc_svf 1976 2291 … … 2089 2404 INTEGER(iwp), DIMENSION(0:nd-1) :: dirend 2090 2405 INTEGER(iwp) :: ids,isurf,isvf,isurfs,isurflt 2091 INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat 2406 INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m 2407 INTEGER(iwp) :: k_topo !< topography top index 2092 2408 2093 2409 dirstart = (/ startland, startwall, startwall, startwall, startwall /) … … 2139 2455 CASE ( 'usm_surfz' ) 2140 2456 !-- array of lw radiation falling to local surface after i-th reflection 2141 DO isurf = dirstart(ids), dirend(ids) 2142 IF ( surfl(id,isurf) == ids ) THEN 2143 IF ( surfl(id,isurf) == iroof ) THEN 2144 temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) = & 2145 max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)), & 2146 REAL(surfl(iz,isurf),wp)) 2147 ELSE 2148 temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) = & 2149 max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)), & 2150 REAL(surfl(iz,isurf),wp)+1.0_wp) 2151 ENDIF 2152 ENDIF 2457 DO m = 1, surf_usm_h%ns 2458 i = surf_usm_h%i(m) 2459 j = surf_usm_h%j(m) 2460 k = surf_usm_h%k(m) 2461 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) ) 2462 ENDDO 2463 DO l = 0, 3 2464 DO m = 1, surf_usm_v(l)%ns 2465 i = surf_usm_v(l)%i(m) 2466 j = surf_usm_v(l)%j(m) 2467 k = surf_usm_v(l)%k(m) 2468 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp ) 2469 ENDDO 2153 2470 ENDDO 2154 2471 2155 2472 CASE ( 'usm_surfcat' ) 2156 2473 !-- surface category 2157 DO isurf = dirstart(ids), dirend(ids) 2158 IF ( surfl(id,isurf) == ids ) THEN 2159 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surface_types(isurf) 2160 ENDIF 2474 DO m = 1, surf_usm_h%ns 2475 i = surf_usm_h%i(m) 2476 j = surf_usm_h%j(m) 2477 k = surf_usm_h%k(m) 2478 temp_pf(k,j,i) = surf_usm_h%surface_types(m) 2479 ENDDO 2480 DO l = 0, 3 2481 DO m = 1, surf_usm_v(l)%ns 2482 i = surf_usm_v(l)%i(m) 2483 j = surf_usm_v(l)%j(m) 2484 k = surf_usm_v(l)%k(m) 2485 temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m) 2486 ENDDO 2161 2487 ENDDO 2162 2488 2163 2489 CASE ( 'usm_surfalb' ) 2164 2490 !-- surface albedo 2165 DO isurf = dirstart(ids), dirend(ids) 2166 IF ( surfl(id,isurf) == ids ) THEN 2167 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = albedo_surf(isurf) 2168 ENDIF 2491 DO m = 1, surf_usm_h%ns 2492 i = surf_usm_h%i(m) 2493 j = surf_usm_h%j(m) 2494 k = surf_usm_h%k(m) 2495 temp_pf(k,j,i) = surf_usm_h%albedo_surf(m) 2496 ENDDO 2497 DO l = 0, 3 2498 DO m = 1, surf_usm_v(l)%ns 2499 i = surf_usm_v(l)%i(m) 2500 j = surf_usm_v(l)%j(m) 2501 k = surf_usm_v(l)%k(m) 2502 temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m) 2503 ENDDO 2169 2504 ENDDO 2170 2505 2171 2506 CASE ( 'usm_surfemis' ) 2172 2507 !-- surface albedo 2173 DO isurf = dirstart(ids), dirend(ids) 2174 IF ( surfl(id,isurf) == ids ) THEN 2175 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = emiss_surf(isurf) 2176 ENDIF 2508 DO m = 1, surf_usm_h%ns 2509 i = surf_usm_h%i(m) 2510 j = surf_usm_h%j(m) 2511 k = surf_usm_h%k(m) 2512 temp_pf(k,j,i) = surf_usm_h%emiss_surf(m) 2177 2513 ENDDO 2178 2514 DO l = 0, 3 2515 DO m = 1, surf_usm_v(l)%ns 2516 i = surf_usm_v(l)%i(m) 2517 j = surf_usm_v(l)%j(m) 2518 k = surf_usm_v(l)%k(m) 2519 temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m) 2520 ENDDO 2521 ENDDO 2522 ! 2523 !-- Not adjusted so far 2179 2524 CASE ( 'usm_svf', 'usm_dif' ) 2180 2525 !-- shape view factors or iradiance factors to selected surface … … 2197 2542 CASE ( 'usm_rad_net' ) 2198 2543 !-- array of complete radiation balance 2199 DO isurf = dirstart(ids), dirend(ids) 2200 IF ( surfl(id,isurf) == ids ) THEN 2201 IF ( av == 0 ) THEN 2202 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_l(isurf) 2203 ELSE 2204 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_av(isurf) 2205 ENDIF 2206 ENDIF 2207 ENDDO 2544 IF ( av == 0 ) THEN 2545 DO m = 1, surf_usm_h%ns 2546 i = surf_usm_h%i(m) 2547 j = surf_usm_h%j(m) 2548 k = surf_usm_h%k(m) 2549 temp_pf(k,j,i) = surf_usm_h%rad_net_l(m) 2550 ENDDO 2551 DO l = 0, 3 2552 DO m = 1, surf_usm_v(l)%ns 2553 i = surf_usm_v(l)%i(m) 2554 j = surf_usm_v(l)%j(m) 2555 k = surf_usm_v(l)%k(m) 2556 temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m) 2557 ENDDO 2558 ENDDO 2559 ELSE 2560 DO m = 1, surf_usm_h%ns 2561 i = surf_usm_h%i(m) 2562 j = surf_usm_h%j(m) 2563 k = surf_usm_h%k(m) 2564 temp_pf(k,j,i) = surf_usm_h%rad_net_av(m) 2565 ENDDO 2566 DO l = 0, 3 2567 DO m = 1, surf_usm_v(l)%ns 2568 i = surf_usm_v(l)%i(m) 2569 j = surf_usm_v(l)%j(m) 2570 k = surf_usm_v(l)%k(m) 2571 temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m) 2572 ENDDO 2573 ENDDO 2574 ENDIF 2208 2575 2209 2576 CASE ( 'usm_rad_insw' ) … … 2268 2635 ENDDO 2269 2636 2270 CASE ( 'usm_rad_inlwdif' )2271 !-- array of sw radiation falling to surface after i-th reflection2272 DO isurf = dirstart(ids), dirend(ids)2273 IF ( surfl(id,isurf) == ids ) THEN2274 IF ( av == 0 ) THEN2275 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)2276 ELSE2277 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)2278 ENDIF2279 ENDIF2280 ENDDO2281 2282 2637 CASE ( 'usm_rad_inlwref' ) 2283 2638 !-- array of lw radiation falling to surface from reflections … … 2339 2694 ENDIF 2340 2695 ENDDO 2341 2696 2342 2697 CASE ( 'usm_rad_hf' ) 2343 2698 !-- array of heat flux from radiation for surfaces after all reflections 2344 DO isurf = dirstart(ids), dirend(ids) 2345 IF ( surfl(id,isurf) == ids ) THEN 2346 IF ( av == 0 ) THEN 2347 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf(isurf) 2348 ELSE 2349 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf_av(isurf) 2350 ENDIF 2351 ENDIF 2352 ENDDO 2353 2699 IF ( av == 0 ) THEN 2700 DO m = 1, surf_usm_h%ns 2701 i = surf_usm_h%i(m) 2702 j = surf_usm_h%j(m) 2703 k = surf_usm_h%k(m) 2704 temp_pf(k,j,i) = surf_usm_h%surfhf(m) 2705 ENDDO 2706 DO l = 0, 3 2707 DO m = 1, surf_usm_v(l)%ns 2708 i = surf_usm_v(l)%i(m) 2709 j = surf_usm_v(l)%j(m) 2710 k = surf_usm_v(l)%k(m) 2711 temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m) 2712 ENDDO 2713 ENDDO 2714 ELSE 2715 DO m = 1, surf_usm_h%ns 2716 i = surf_usm_h%i(m) 2717 j = surf_usm_h%j(m) 2718 k = surf_usm_h%k(m) 2719 temp_pf(k,j,i) = surf_usm_h%surfhf_av(m) 2720 ENDDO 2721 DO l = 0, 3 2722 DO m = 1, surf_usm_v(l)%ns 2723 i = surf_usm_v(l)%i(m) 2724 j = surf_usm_v(l)%j(m) 2725 k = surf_usm_v(l)%k(m) 2726 temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m) 2727 ENDDO 2728 ENDDO 2729 ENDIF 2730 2354 2731 CASE ( 'usm_wshf' ) 2355 2732 !-- array of sensible heat flux from surfaces 2356 !-- horizontal surfaces 2357 DO isurf = dirstart(ids), dirend(ids) 2358 IF ( surfl(id,isurf) == ids ) THEN 2359 IF ( av == 0 ) THEN 2360 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb(isurf) 2361 ELSE 2362 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb_av(isurf) 2363 ENDIF 2364 ENDIF 2365 ENDDO 2733 IF ( av == 0 ) THEN 2734 DO m = 1, surf_usm_h%ns 2735 i = surf_usm_h%i(m) 2736 j = surf_usm_h%j(m) 2737 k = surf_usm_h%k(m) 2738 temp_pf(k,j,i) = surf_usm_h%wshf_eb(m) 2739 ENDDO 2740 DO l = 0, 3 2741 DO m = 1, surf_usm_v(l)%ns 2742 i = surf_usm_v(l)%i(m) 2743 j = surf_usm_v(l)%j(m) 2744 k = surf_usm_v(l)%k(m) 2745 temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m) 2746 ENDDO 2747 ENDDO 2748 ELSE 2749 DO m = 1, surf_usm_h%ns 2750 i = surf_usm_h%i(m) 2751 j = surf_usm_h%j(m) 2752 k = surf_usm_h%k(m) 2753 temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m) 2754 ENDDO 2755 DO l = 0, 3 2756 DO m = 1, surf_usm_v(l)%ns 2757 i = surf_usm_v(l)%i(m) 2758 j = surf_usm_v(l)%j(m) 2759 k = surf_usm_v(l)%k(m) 2760 temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m) 2761 ENDDO 2762 ENDDO 2763 ENDIF 2764 2366 2765 2367 2766 CASE ( 'usm_wghf' ) 2368 2767 !-- array of heat flux from ground (land, wall, roof) 2369 DO isurf = dirstart(ids), dirend(ids) 2370 IF ( surfl(id,isurf) == ids ) THEN 2371 IF ( av == 0 ) THEN 2372 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb(isurf) 2373 ELSE 2374 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb_av(isurf) 2375 ENDIF 2376 ENDIF 2377 ENDDO 2768 IF ( av == 0 ) THEN 2769 DO m = 1, surf_usm_h%ns 2770 i = surf_usm_h%i(m) 2771 j = surf_usm_h%j(m) 2772 k = surf_usm_h%k(m) 2773 temp_pf(k,j,i) = surf_usm_h%wghf_eb(m) 2774 ENDDO 2775 DO l = 0, 3 2776 DO m = 1, surf_usm_v(l)%ns 2777 i = surf_usm_v(l)%i(m) 2778 j = surf_usm_v(l)%j(m) 2779 k = surf_usm_v(l)%k(m) 2780 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m) 2781 ENDDO 2782 ENDDO 2783 ELSE 2784 DO m = 1, surf_usm_h%ns 2785 i = surf_usm_h%i(m) 2786 j = surf_usm_h%j(m) 2787 k = surf_usm_h%k(m) 2788 temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m) 2789 ENDDO 2790 DO l = 0, 3 2791 DO m = 1, surf_usm_v(l)%ns 2792 i = surf_usm_v(l)%i(m) 2793 j = surf_usm_v(l)%j(m) 2794 k = surf_usm_v(l)%k(m) 2795 temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m) 2796 ENDDO 2797 ENDDO 2798 ENDIF 2378 2799 2379 2800 CASE ( 'usm_t_surf' ) 2380 2801 !-- surface temperature for surfaces 2381 DO isurf = max(startenergy,dirstart(ids)), min(endenergy,dirend(ids)) 2382 IF ( surfl(id,isurf) == ids ) THEN 2383 IF ( av == 0 ) THEN 2384 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf(isurf) 2385 ELSE 2386 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf_av(isurf) 2387 ENDIF 2388 ENDIF 2389 ENDDO 2802 IF ( av == 0 ) THEN 2803 DO m = 1, surf_usm_h%ns 2804 i = surf_usm_h%i(m) 2805 j = surf_usm_h%j(m) 2806 k = surf_usm_h%k(m) 2807 temp_pf(k,j,i) = t_surf_h(m) 2808 ENDDO 2809 DO l = 0, 3 2810 DO m = 1, surf_usm_v(l)%ns 2811 i = surf_usm_v(l)%i(m) 2812 j = surf_usm_v(l)%j(m) 2813 k = surf_usm_v(l)%k(m) 2814 temp_pf(k,j,i) = t_surf_v(l)%t(m) 2815 ENDDO 2816 ENDDO 2817 ELSE 2818 DO m = 1, surf_usm_h%ns 2819 i = surf_usm_h%i(m) 2820 j = surf_usm_h%j(m) 2821 k = surf_usm_h%k(m) 2822 temp_pf(k,j,i) = surf_usm_h%t_surf_av(m) 2823 ENDDO 2824 DO l = 0, 3 2825 DO m = 1, surf_usm_v(l)%ns 2826 i = surf_usm_v(l)%i(m) 2827 j = surf_usm_v(l)%j(m) 2828 k = surf_usm_v(l)%k(m) 2829 temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m) 2830 ENDDO 2831 ENDDO 2832 ENDIF 2390 2833 2391 2834 CASE ( 'usm_t_wall' ) 2392 2835 !-- wall temperature for iwl layer of walls and land 2393 DO isurf = dirstart(ids), dirend(ids) 2394 IF ( surfl(id,isurf) == ids ) THEN 2395 IF ( av == 0 ) THEN 2396 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall(iwl,isurf) 2397 ELSE 2398 temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall_av(iwl,isurf) 2399 ENDIF 2400 ENDIF 2401 ENDDO 2836 IF ( av == 0 ) THEN 2837 DO m = 1, surf_usm_h%ns 2838 i = surf_usm_h%i(m) 2839 j = surf_usm_h%j(m) 2840 k = surf_usm_h%k(m) 2841 temp_pf(k,j,i) = t_wall_h(iwl,m) 2842 ENDDO 2843 DO l = 0, 3 2844 DO m = 1, surf_usm_v(l)%ns 2845 i = surf_usm_v(l)%i(m) 2846 j = surf_usm_v(l)%j(m) 2847 k = surf_usm_v(l)%k(m) 2848 temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m) 2849 ENDDO 2850 ENDDO 2851 ELSE 2852 DO m = 1, surf_usm_h%ns 2853 i = surf_usm_h%i(m) 2854 j = surf_usm_h%j(m) 2855 k = surf_usm_h%k(m) 2856 temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m) 2857 ENDDO 2858 DO l = 0, 3 2859 DO m = 1, surf_usm_v(l)%ns 2860 i = surf_usm_v(l)%i(m) 2861 j = surf_usm_v(l)%j(m) 2862 k = surf_usm_v(l)%k(m) 2863 temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m) 2864 ENDDO 2865 ENDDO 2866 ENDIF 2402 2867 2403 2868 CASE DEFAULT … … 2408 2873 !-- fill out array local_pf which is subsequently treated by data_output_3d 2409 2874 CALL exchange_horiz( temp_pf, nbgp ) 2875 ! 2876 !-- To Do: why reversed loop order 2410 2877 DO j = nysg,nyng 2411 2878 DO i = nxlg,nxrg … … 2589 3056 IMPLICIT NONE 2590 3057 2591 INTEGER(iwp) :: k, l !< running indices3058 INTEGER(iwp) :: k, l, m !< running indices 2592 3059 2593 3060 CALL location_message( ' initialization of wall surface model', .TRUE. ) … … 2596 3063 !-- Temperature is defined at the center of the wall layers, 2597 3064 !-- whereas gradients/fluxes are defined at the edges (_stag) 2598 DO l= nzb_wall, nzt_wall2599 zwn( l) = zwn_default(l)3065 DO k = nzb_wall, nzt_wall 3066 zwn(k) = zwn_default(k) 2600 3067 ENDDO 2601 2602 !-- apply for all particular wall grids 2603 DO l = startenergy, endenergy 2604 zw(:,l) = zwn(:) * thickness_wall(l) 2605 dz_wall(nzb_wall,l) = zw(nzb_wall,l) 3068 ! 3069 !-- apply for all particular surface grids. First for horizontal surfaces 3070 DO m = 1, surf_usm_h%ns 3071 surf_usm_h%zw(:,m) = zwn(:) * & 3072 surf_usm_h%thickness_wall(m) 3073 surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m) 2606 3074 DO k = nzb_wall+1, nzt_wall 2607 dz_wall(k,l) = zw(k,l) - zw(k-1,l) 3075 surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - & 3076 surf_usm_h%zw(k-1,m) 2608 3077 ENDDO 2609 3078 2610 dz_wall(nzt_wall+1,l) = dz_wall(nzt_wall,l)3079 surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m) 2611 3080 2612 3081 DO k = nzb_wall, nzt_wall-1 2613 dz_wall_stag(k,l) = 0.5 * (dz_wall(k+1,l) + dz_wall(k,l)) 3082 surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( & 3083 surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) ) 2614 3084 ENDDO 2615 dz_wall_stag(nzt_wall,l) = dz_wall(nzt_wall,l)3085 surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m) 2616 3086 ENDDO 2617 2618 ddz_wall = 1.0_wp / dz_wall 2619 ddz_wall_stag = 1.0_wp / dz_wall_stag 3087 surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall 3088 surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag 3089 ! 3090 !-- For vertical surfaces 3091 DO l = 0, 3 3092 DO m = 1, surf_usm_v(l)%ns 3093 surf_usm_v(l)%zw(:,m) = zwn(:) * & 3094 surf_usm_v(l)%thickness_wall(m) 3095 surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m) 3096 DO k = nzb_wall+1, nzt_wall 3097 surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - & 3098 surf_usm_v(l)%zw(k-1,m) 3099 ENDDO 3100 3101 surf_usm_v(l)%dz_wall(nzt_wall+1,m) = & 3102 surf_usm_v(l)%dz_wall(nzt_wall,m) 3103 3104 DO k = nzb_wall, nzt_wall-1 3105 surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( & 3106 surf_usm_v(l)%dz_wall(k+1,m) + & 3107 surf_usm_v(l)%dz_wall(k,m) ) 3108 ENDDO 3109 surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = & 3110 surf_usm_v(l)%dz_wall(nzt_wall,m) 3111 ENDDO 3112 surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall 3113 surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag 3114 ENDDO 3115 2620 3116 2621 3117 CALL location_message( ' wall structures filed out', .TRUE. ) … … 2635 3131 IMPLICIT NONE 2636 3132 2637 INTEGER(iwp) :: i, j, k, l !< running indices3133 INTEGER(iwp) :: i, j, k, l, m !< running indices 2638 3134 REAL(wp) :: c, d, tin, exn 2639 3135 … … 2699 3195 ELSE 2700 3196 2701 !-- Calculate initial surface temperature 3197 !-- Calculate initial surface temperature from pt of adjacent gridbox 2702 3198 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 2703 3199 2704 DO l = startenergy, endenergy 2705 k = surfl(iz,l) 2706 j = surfl(iy,l) 2707 i = surfl(ix,l) 2708 2709 !-- Initial surface temperature set from pt of adjacent gridbox 2710 t_surf(l) = pt(k,j,i) * exn 3200 ! 3201 !-- At horizontal surfaces. Please note, t_surf_h is defined on a 3202 !-- different data type, but with the same dimension. 3203 DO m = 1, surf_usm_h%ns 3204 i = surf_usm_h%i(m) 3205 j = surf_usm_h%j(m) 3206 k = surf_usm_h%k(m) 3207 3208 t_surf_h(m) = pt(k,j,i) * exn 2711 3209 ENDDO 3210 ! 3211 !-- At vertical surfaces. 3212 DO l = 0, 3 3213 DO m = 1, surf_usm_v(l)%ns 3214 i = surf_usm_v(l)%i(m) 3215 j = surf_usm_v(l)%j(m) 3216 k = surf_usm_v(l)%k(m) 3217 3218 t_surf_v(l)%t(m) = pt(k,j,i) * exn 3219 ENDDO 3220 ENDDO 3221 2712 3222 2713 3223 !-- initial values for t_wall 2714 3224 !-- outer value is set to surface temperature 2715 3225 !-- inner value is set to wall_inner_temperature 2716 !-- and profile is logaritmic (linear in nz) 2717 DO l = startenergy, endenergy 2718 IF ( isroof_surf(l) ) THEN 2719 tin = roof_inner_temperature 2720 ELSE IF ( surf(id,l)==iroof ) THEN 2721 tin = soil_inner_temperature 2722 ELSE 2723 tin = wall_inner_temperature 2724 ENDIF 2725 DO k = nzb_wall, nzt_wall+1 2726 c = REAL(k-nzb_wall,wp)/REAL(nzt_wall+1-nzb_wall,wp) 2727 t_wall(k,:) = (1.0_wp-c)*t_surf(:) + c*tin 2728 ENDDO 3226 !-- and profile is logaritmic (linear in nz). 3227 !-- Horizontal surfaces 3228 DO m = 1, surf_usm_h%ns 3229 ! 3230 !-- Roof 3231 IF ( surf_usm_h%isroof_surf(m) ) THEN 3232 tin = roof_inner_temperature 3233 ! 3234 !-- Normal land surface 3235 ELSE 3236 tin = soil_inner_temperature 3237 ENDIF 3238 3239 DO k = nzb_wall, nzt_wall+1 3240 c = REAL( k - nzb_wall, wp ) / & 3241 REAL( nzt_wall + 1 - nzb_wall , wp ) 3242 3243 t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin 3244 ENDDO 2729 3245 ENDDO 3246 ! 3247 !-- Vertical surfaces 3248 DO l = 0, 3 3249 DO m = 1, surf_usm_v(l)%ns 3250 ! 3251 !-- Inner wall 3252 tin = wall_inner_temperature 3253 3254 DO k = nzb_wall, nzt_wall+1 3255 c = REAL( k - nzb_wall, wp ) / & 3256 REAL( nzt_wall + 1 - nzb_wall , wp ) 3257 3258 t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) + & 3259 c * tin 3260 ENDDO 3261 ENDDO 3262 ENDDO 3263 2730 3264 ENDIF 2731 3265 2732 3266 !-- 2733 !-- 3267 !-- Possibly DO user-defined actions (e.g. define heterogeneous wall surface) 2734 3268 CALL user_init_urban_surface 2735 3269 2736 3270 !-- initialize prognostic values for the first timestep 2737 t_surf_p = t_surf 2738 t_wall_p = t_wall 3271 t_surf_h_p = t_surf_h 3272 t_surf_v_p = t_surf_v 3273 3274 t_wall_h_p = t_wall_h 3275 t_wall_v_p = t_wall_v 2739 3276 2740 3277 !-- Adjust radiative fluxes for urban surface at model start … … 2759 3296 IMPLICIT NONE 2760 3297 2761 INTEGER(iwp) :: i,j,k,l,kw !< running indices3298 INTEGER(iwp) :: i,j,k,l,kw, m !< running indices 2762 3299 2763 3300 REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend !< tendency 2764 3301 2765 2766 DO l = startenergy, endenergy 2767 !-- calculate frequently used parameters 2768 k = surfl(iz,l) 2769 j = surfl(iy,l) 2770 i = surfl(ix,l) 2771 2772 ! 2773 !-- prognostic equation for ground/wall/roof temperature t_wall 2774 wtend(:) = 0.0_wp 2775 wtend(nzb_wall) = (1.0_wp/rho_c_wall(nzb_wall,l)) * & 2776 ( lambda_h(nzb_wall,l) * ( t_wall(nzb_wall+1,l) & 2777 - t_wall(nzb_wall,l) ) * ddz_wall(nzb_wall+1,l) & 2778 + wghf_eb(l) ) * ddz_wall_stag(nzb_wall,l) 3302 ! 3303 !-- For horizontal surfaces 3304 DO m = 1, surf_usm_h%ns 3305 ! 3306 !-- Obtain indices 3307 i = surf_usm_h%i(m) 3308 j = surf_usm_h%j(m) 3309 k = surf_usm_h%k(m) 3310 ! 3311 !-- prognostic equation for ground/roof temperature t_wall_h 3312 wtend(:) = 0.0_wp 3313 wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) * & 3314 ( surf_usm_h%lambda_h(nzb_wall,m) * & 3315 ( t_wall_h(nzb_wall+1,m) & 3316 - t_wall_h(nzb_wall,m) ) * & 3317 surf_usm_h%ddz_wall(nzb_wall+1,m) & 3318 + surf_usm_h%wghf_eb(m) ) * & 3319 surf_usm_h%ddz_wall_stag(nzb_wall,m) 2779 3320 2780 2781 wtend(kw) = (1.0_wp/rho_c_wall(kw,l))&2782 * ( lambda_h(kw,l)&2783 * ( t_wall (kw+1,l) - t_wall(kw,l) )&2784 * ddz_wall(kw+1,l)&2785 - lambda_h(kw-1,l)&2786 * ( t_wall (kw,l) - t_wall(kw-1,l) )&2787 * ddz_wall(kw,l)&2788 ) * ddz_wall_stag(kw,l)3321 DO kw = nzb_wall+1, nzt_wall 3322 wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m)) & 3323 * ( surf_usm_h%lambda_h(kw,m) & 3324 * ( t_wall_h(kw+1,m) - t_wall_h(kw,m) ) & 3325 * surf_usm_h%ddz_wall(kw+1,m) & 3326 - surf_usm_h%lambda_h(kw-1,m) & 3327 * ( t_wall_h(kw,m) - t_wall_h(kw-1,m) ) & 3328 * surf_usm_h%ddz_wall(kw,m) & 3329 ) * surf_usm_h%ddz_wall_stag(kw,m) 2789 3330 ENDDO 2790 3331 2791 t_wall_p(nzb_wall:nzt_wall,l) = t_wall(nzb_wall:nzt_wall,l)&2792 + dt_3d * ( tsc(2)&2793 * wtend(nzb_wall:nzt_wall) + tsc(3)&2794 * tt_wall_m(nzb_wall:nzt_wall,l) )3332 t_wall_h_p(nzb_wall:nzt_wall,m) = t_wall_h(nzb_wall:nzt_wall,m) & 3333 + dt_3d * ( tsc(2) & 3334 * wtend(nzb_wall:nzt_wall) + tsc(3) & 3335 * surf_usm_h%tt_wall_m(nzb_wall:nzt_wall,m) ) 2795 3336 2796 ! 2797 !-- calculate t_wall tendencies for the next Runge-Kutta step 2798 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2799 IF ( intermediate_timestep_count == 1 ) THEN 3337 ! 3338 !-- calculate t_wall tendencies for the next Runge-Kutta step 3339 IF ( timestep_scheme(1:5) == 'runge' ) THEN 3340 IF ( intermediate_timestep_count == 1 ) THEN 3341 DO kw = nzb_wall, nzt_wall 3342 surf_usm_h%tt_wall_m(kw,m) = wtend(kw) 3343 ENDDO 3344 ELSEIF ( intermediate_timestep_count < & 3345 intermediate_timestep_count_max ) THEN 2800 3346 DO kw = nzb_wall, nzt_wall 2801 tt_wall_m(kw,l) = wtend(kw) 3347 surf_usm_h%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) + & 3348 5.3125_wp * surf_usm_h%tt_wall_m(kw,m) 2802 3349 ENDDO 2803 ELSEIF ( intermediate_timestep_count < & 2804 intermediate_timestep_count_max ) THEN 2805 DO kw = nzb_wall, nzt_wall 2806 tt_wall_m(kw,l) = -9.5625_wp * wtend(kw) + 5.3125_wp & 2807 * tt_wall_m(kw,l) 2808 ENDDO 2809 ENDIF 2810 ENDIF 3350 ENDIF 3351 ENDIF 3352 ENDDO 3353 ! 3354 !-- For vertical surfaces 3355 DO l = 0, 3 3356 DO m = 1, surf_usm_v(l)%ns 3357 ! 3358 !-- Obtain indices 3359 i = surf_usm_v(l)%i(m) 3360 j = surf_usm_v(l)%j(m) 3361 k = surf_usm_v(l)%k(m) 3362 ! 3363 !-- prognostic equation for wall temperature t_wall_v 3364 wtend(:) = 0.0_wp 3365 wtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_wall(nzb_wall,m)) * & 3366 ( surf_usm_v(l)%lambda_h(nzb_wall,m) * & 3367 ( t_wall_v(l)%t(nzb_wall+1,m) & 3368 - t_wall_v(l)%t(nzb_wall,m) ) * & 3369 surf_usm_v(l)%ddz_wall(nzb_wall+1,m) & 3370 + surf_usm_v(l)%wghf_eb(m) ) * & 3371 surf_usm_v(l)%ddz_wall_stag(nzb_wall,m) 3372 3373 DO kw = nzb_wall+1, nzt_wall 3374 wtend(kw) = (1.0_wp / surf_usm_v(l)%rho_c_wall(kw,m)) & 3375 * ( surf_usm_v(l)%lambda_h(kw,m) & 3376 * ( t_wall_v(l)%t(kw+1,m) - t_wall_v(l)%t(kw,m) )& 3377 * surf_usm_v(l)%ddz_wall(kw+1,m) & 3378 - surf_usm_v(l)%lambda_h(kw-1,m) & 3379 * ( t_wall_v(l)%t(kw,m) - t_wall_v(l)%t(kw-1,m) )& 3380 * surf_usm_v(l)%ddz_wall(kw,m) & 3381 ) * surf_usm_v(l)%ddz_wall_stag(kw,m) 3382 ENDDO 3383 3384 t_wall_v_p(l)%t(nzb_wall:nzt_wall,m) = & 3385 t_wall_v(l)%t(nzb_wall:nzt_wall,m) & 3386 + dt_3d * ( tsc(2) & 3387 * wtend(nzb_wall:nzt_wall) + tsc(3) & 3388 * surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall,m) ) 3389 3390 ! 3391 !-- calculate t_wall tendencies for the next Runge-Kutta step 3392 IF ( timestep_scheme(1:5) == 'runge' ) THEN 3393 IF ( intermediate_timestep_count == 1 ) THEN 3394 DO kw = nzb_wall, nzt_wall 3395 surf_usm_v(l)%tt_wall_m(kw,m) = wtend(kw) 3396 ENDDO 3397 ELSEIF ( intermediate_timestep_count < & 3398 intermediate_timestep_count_max ) THEN 3399 DO kw = nzb_wall, nzt_wall 3400 surf_usm_v(l)%tt_wall_m(kw,m) = & 3401 - 9.5625_wp * wtend(kw) + & 3402 5.3125_wp * surf_usm_v(l)%tt_wall_m(kw,m) 3403 ENDDO 3404 ENDIF 3405 ENDIF 3406 ENDDO 2811 3407 ENDDO 2812 3408 … … 2857 3453 !-- Read user-defined namelist 2858 3454 READ ( 11, urban_surface_par ) 2859 2860 3455 ! 2861 3456 !-- Set flag that indicates that the land surface model is switched on 2862 3457 urban_surface = .TRUE. 2863 2864 3458 2865 3459 10 CONTINUE … … 2879 3473 IMPLICIT NONE 2880 3474 2881 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep 3475 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll 2882 3476 INTEGER(iwp) :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb 2883 3477 INTEGER(iwp), DIMENSION(4) :: bdycross … … 2892 3486 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 2893 3487 INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy 3488 INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /) 2894 3489 2895 3490 … … 2949 3544 !-- First pass: direct + diffuse irradiance 2950 3545 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2951 surfinswdir = 0._wp 2952 surfinswdif = 0._wp 2953 surfinlwdif = 0._wp 2954 surfins = 0._wp2955 surfinl = 0._wp2956 surfoutsl = 0._wp2957 surfoutll = 0._wp3546 surfinswdir = 0._wp !nsurfl 3547 surfinswdif = 0._wp !nsurfl 3548 surfinlwdif = 0._wp !nsurfl 3549 surfins = 0._wp !nsurfl 3550 surfinl = 0._wp !nsurfl 3551 surfoutsl(:) = 0.0_wp !start-end 3552 surfoutll(:) = 0.0_wp !start-end 2958 3553 2959 3554 !-- Set up thermal radiation from surfaces 2960 3555 !-- emiss_surf is defined only for surfaces for which energy balance is calculated 2961 surfoutll(startenergy:endenergy) = emiss_surf(startenergy:endenergy) * sigma_sb & 2962 * t_surf(startenergy:endenergy)**4 3556 !-- Workaround: reorder surface data type back on 1D array including all surfaces, 3557 !-- which implies to reorder horizontal and vertical surfaces 3558 ! 3559 !-- Horizontal walls 3560 mm = 1 3561 DO i = nxl, nxr 3562 DO j = nys, nyn 3563 3564 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 3565 surfoutll(mm) = surf_usm_h%emiss_surf(m) * sigma_sb & 3566 * t_surf_h(m)**4 3567 albedo_surf(mm) = surf_usm_h%albedo_surf(m) 3568 emiss_surf(mm) = surf_usm_h%emiss_surf(m) 3569 mm = mm + 1 3570 ENDDO 3571 ENDDO 3572 ENDDO 3573 ! 3574 !-- Vertical walls 3575 DO i = nxl, nxr 3576 DO j = nys, nyn 3577 DO ll = 0, 3 3578 l = reorder(ll) 3579 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 3580 surfoutll(mm) = surf_usm_v(l)%emiss_surf(m) * sigma_sb & 3581 * t_surf_v(l)%t(m)**4 3582 albedo_surf(mm) = surf_usm_v(l)%albedo_surf(m) 3583 emiss_surf(mm) = surf_usm_v(l)%emiss_surf(m) 3584 mm = mm + 1 3585 ENDDO 3586 ENDDO 3587 ENDDO 3588 ENDDO 2963 3589 2964 3590 #if defined( __parallel ) 2965 3591 !-- might be optimized and gather only values relevant for current processor 3592 2966 3593 CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, & 2967 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 3594 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global 2968 3595 #else 2969 surfoutl(:) = surfoutll(:) 3596 surfoutl(:) = surfoutll(:) !nsurf global 2970 3597 #endif 2971 3598 … … 3005 3632 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 3006 3633 ENDIF 3007 3008 IF ( zenith(0) > 0 .AND. all( surf( :,isurfsrc) == bdycross ) ) THEN3634 3635 IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN 3009 3636 !-- found svf between model boundary and the face => face isn't shaded 3010 surfinswdir(isurf) = rad_sw_in_dir(j, 3637 surfinswdir(isurf) = rad_sw_in_dir(j,i) & 3011 3638 * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0) 3012 3639 … … 3050 3677 ENDIF 3051 3678 3052 IF ( zenith(0) > 0 .AND. all( surf( :,isurfsrc) == bdycross ) ) THEN3679 IF ( zenith(0) > 0 .AND. all( surf(1:4,isurfsrc) == bdycross ) ) THEN 3053 3680 !-- found svf between model boundary and the pcgb => pcgb isn't shaded 3054 3681 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) … … 3058 3685 ENDDO 3059 3686 ENDIF 3687 3060 3688 surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy) 3061 3689 surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy) … … 3064 3692 surfoutsw(:) = 0.0_wp 3065 3693 surfoutlw(:) = surfoutll(:) 3066 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &3067 - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)3694 ! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 3695 ! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy) 3068 3696 3069 3697 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 3117 3745 surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy) 3118 3746 surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy) 3119 surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &3120 - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)3747 ! surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) & 3748 ! - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy) 3121 3749 3122 3750 ENDDO … … 3129 3757 i = pcbl(ix, ipcgb) 3130 3758 k = pcbl(iz, ipcgb) 3131 kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat 3759 kk = k - MAXLOC( & 3760 MERGE( 1, 0, & 3761 BTEST( wall_flags_0(:,j,i), 12 ) & 3762 ), DIM = 1 & 3763 ) - 1 ! kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat 3132 3764 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 3133 3765 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt 3134 3766 ENDDO 3135 3767 ENDIF 3768 ! 3769 !-- Transfer radiation arrays required for energy balance to the respective data types 3770 DO i = startenergy, endenergy 3771 m = surfl(5,i) 3772 ! 3773 !-- upward-facing 3774 IF ( surfl(1,i) == 0 ) THEN 3775 surf_usm_h%rad_in_sw(m) = surfinsw(i) 3776 surf_usm_h%rad_out_sw(m) = surfoutsw(i) 3777 surf_usm_h%rad_in_lw(m) = surfinlw(i) 3778 surf_usm_h%rad_out_lw(m) = surfoutlw(i) 3779 ! 3780 !-- southward-facding 3781 ELSEIF ( surfl(1,i) == 1 ) THEN 3782 surf_usm_v(1)%rad_in_sw(m) = surfinsw(i) 3783 surf_usm_v(1)%rad_out_sw(m) = surfoutsw(i) 3784 surf_usm_v(1)%rad_in_lw(m) = surfinlw(i) 3785 surf_usm_v(1)%rad_out_lw(m) = surfoutlw(i) 3786 ! 3787 !-- northward-facding 3788 ELSEIF ( surfl(1,i) == 2 ) THEN 3789 surf_usm_v(0)%rad_in_sw(m) = surfinsw(i) 3790 surf_usm_v(0)%rad_out_sw(m) = surfoutsw(i) 3791 surf_usm_v(0)%rad_in_lw(m) = surfinlw(i) 3792 surf_usm_v(0)%rad_out_lw(m) = surfoutlw(i) 3793 ! 3794 !-- westward-facding 3795 ELSEIF ( surfl(1,i) == 3 ) THEN 3796 surf_usm_v(3)%rad_in_sw(m) = surfinsw(i) 3797 surf_usm_v(3)%rad_out_sw(m) = surfoutsw(i) 3798 surf_usm_v(3)%rad_in_lw(m) = surfinlw(i) 3799 surf_usm_v(3)%rad_out_lw(m) = surfoutlw(i) 3800 ! 3801 !-- eastward-facing 3802 ELSEIF ( surfl(1,i) == 4 ) THEN 3803 surf_usm_v(2)%rad_in_sw(m) = surfinsw(i) 3804 surf_usm_v(2)%rad_out_sw(m) = surfoutsw(i) 3805 surf_usm_v(2)%rad_in_lw(m) = surfinlw(i) 3806 surf_usm_v(2)%rad_out_lw(m) = surfoutlw(i) 3807 ENDIF 3808 3809 ENDDO 3810 3811 3812 DO m = 1, surf_usm_h%ns 3813 surf_usm_h%surfhf(m) = surf_usm_h%rad_in_sw(m) + & 3814 surf_usm_h%rad_in_lw(m) - & 3815 surf_usm_h%rad_out_sw(m) - & 3816 surf_usm_h%rad_out_lw(m) 3817 ENDDO 3818 3819 DO l = 0, 3 3820 DO m = 1, surf_usm_v(l)%ns 3821 surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_in_sw(m) + & 3822 surf_usm_v(l)%rad_in_lw(m) - & 3823 surf_usm_v(l)%rad_out_sw(m) - & 3824 surf_usm_v(l)%rad_out_lw(m) 3825 ENDDO 3826 ENDDO 3136 3827 3137 3828 !-- return surface radiation to horizontal surfaces … … 3207 3898 REAL(wp), PARAMETER :: grow_factor = 1.5_wp !< factor of expansion of grow arrays 3208 3899 3209 3900 ! 3210 3901 !-- Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also 3211 3902 !-- the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor. … … 3351 4042 visible = .TRUE. 3352 4043 3353 3354 4044 END SUBROUTINE usm_raytrace 3355 4045 … … 3470 4160 SELECT CASE ( TRIM( variable_chr ) ) 3471 4161 3472 CASE ( 't_surf ' )4162 CASE ( 't_surf_h' ) 3473 4163 #if defined( __nopointer ) 3474 IF ( .NOT. ALLOCATED( t_surf ) )&3475 ALLOCATE( t_surf (startenergy:endenergy) )3476 READ ( 13 ) t_surf 4164 IF ( .NOT. ALLOCATED( t_surf_h ) ) & 4165 ALLOCATE( t_surf_h(1:surf_usm_h%ns) ) 4166 READ ( 13 ) t_surf_h 3477 4167 #else 3478 IF ( .NOT. ALLOCATED( t_surf_ 1 ) )&3479 ALLOCATE( t_surf_ 1(startenergy:endenergy) )3480 READ ( 13 ) t_surf_ 14168 IF ( .NOT. ALLOCATED( t_surf_h_1 ) ) & 4169 ALLOCATE( t_surf_h_1(1:surf_usm_h%ns) ) 4170 READ ( 13 ) t_surf_h_1 3481 4171 #endif 3482 3483 CASE ( 't_wall' ) 4172 CASE ( 't_surf_v(0)' ) 4173 #if defined( __nopointer ) 4174 IF ( .NOT. ALLOCATED( t_surf_v(0)%t ) ) & 4175 ALLOCATE( t_surf_v(0)%t(1:surf_usm_v(0)%ns) ) 4176 READ ( 13 ) t_surf_v(0)%t 4177 #else 4178 IF ( .NOT. ALLOCATED( t_surf_v_1(0)%t ) ) & 4179 ALLOCATE( t_surf_v_1(0)%t(1:surf_usm_v(0)%ns) ) 4180 READ ( 13 ) t_surf_v_1(0)%t 4181 #endif 4182 CASE ( 't_surf_v(1)' ) 4183 #if defined( __nopointer ) 4184 IF ( .NOT. ALLOCATED( t_surf_v(1)%t ) ) & 4185 ALLOCATE( t_surf_v(1)%t(1:surf_usm_v(1)%ns) ) 4186 READ ( 13 ) t_surf_v(1)%t 4187 #else 4188 IF ( .NOT. ALLOCATED( t_surf_v_1(1)%t ) ) & 4189 ALLOCATE( t_surf_v_1(1)%t(1:surf_usm_v(1)%ns) ) 4190 READ ( 13 ) t_surf_v_1(1)%t 4191 #endif 4192 CASE ( 't_surf_v(2)' ) 4193 #if defined( __nopointer ) 4194 IF ( .NOT. ALLOCATED( t_surf_v(2)%t ) ) & 4195 ALLOCATE( t_surf_v(2)%t(1:surf_usm_v(2)%ns) ) 4196 READ ( 13 ) t_surf_v(2)%t 4197 #else 4198 IF ( .NOT. ALLOCATED( t_surf_v_1(2)%t ) ) & 4199 ALLOCATE( t_surf_v_1(2)%t(1:surf_usm_v(2)%ns) ) 4200 READ ( 13 ) t_surf_v_1(2)%t 4201 #endif 4202 CASE ( 't_surf_v(3)' ) 4203 #if defined( __nopointer ) 4204 IF ( .NOT. ALLOCATED( t_surf_v(3)%t ) ) & 4205 ALLOCATE( t_surf_v(3)%t(1:surf_usm_v(3)%ns) ) 4206 READ ( 13 ) t_surf_v(3)%t 4207 #else 4208 IF ( .NOT. ALLOCATED( t_surf_v_1(3)%t ) ) & 4209 ALLOCATE( t_surf_v_1(3)%t(1:surf_usm_v(3)%ns) ) 4210 READ ( 13 ) t_surf_v_1(3)%t 4211 #endif 4212 CASE ( 't_wall_h' ) 3484 4213 #if defined( __nopointer ) 3485 IF ( .NOT. ALLOCATED( t_wall ) )&3486 ALLOCATE( t_wall (nzb_wall:nzt_wall+1,startenergy:endenergy) )3487 READ ( 13 ) t_wall 4214 IF ( .NOT. ALLOCATED( t_wall_h ) ) & 4215 ALLOCATE( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 4216 READ ( 13 ) t_wall_h 3488 4217 #else 3489 IF ( .NOT. ALLOCATED( t_wall_1 ) ) & 3490 ALLOCATE( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) ) 3491 READ ( 13 ) t_wall_1 4218 IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) & 4219 ALLOCATE( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 4220 READ ( 13 ) t_wall_h_1 4221 #endif 4222 CASE ( 't_wall_v(0)' ) 4223 #if defined( __nopointer ) 4224 IF ( .NOT. ALLOCATED( t_wall_v(0)%t ) ) & 4225 ALLOCATE( t_wall_v(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) ) 4226 READ ( 13 ) t_wall_v(0)%t 4227 #else 4228 IF ( .NOT. ALLOCATED( t_wall_v_1(0)%t ) ) & 4229 ALLOCATE( t_wall_v_1(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) ) 4230 READ ( 13 ) t_wall_v_1(0)%t 4231 #endif 4232 CASE ( 't_wall_v(1)' ) 4233 #if defined( __nopointer ) 4234 IF ( .NOT. ALLOCATED( t_wall_v(1)%t ) ) & 4235 ALLOCATE( t_wall_v(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) ) 4236 READ ( 13 ) t_wall_v(1)%t 4237 #else 4238 IF ( .NOT. ALLOCATED( t_wall_v_1(0)%t ) ) & 4239 ALLOCATE( t_wall_v_1(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) ) 4240 READ ( 13 ) t_wall_v_1(1)%t 4241 #endif 4242 CASE ( 't_wall_v(2)' ) 4243 #if defined( __nopointer ) 4244 IF ( .NOT. ALLOCATED( t_wall_v(2)%t ) ) & 4245 ALLOCATE( t_wall_v(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) ) 4246 READ ( 13 ) t_wall_v(2)%t 4247 #else 4248 IF ( .NOT. ALLOCATED( t_wall_v_1(2)%t ) ) & 4249 ALLOCATE( t_wall_v_1(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) ) 4250 READ ( 13 ) t_wall_v_1(2)%t 4251 #endif 4252 CASE ( 't_wall_v(3)' ) 4253 #if defined( __nopointer ) 4254 IF ( .NOT. ALLOCATED( t_wall_v(3)%t ) ) & 4255 ALLOCATE( t_wall_v(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) ) 4256 READ ( 13 ) t_wall_v(3)%t 4257 #else 4258 IF ( .NOT. ALLOCATED( t_wall_v_1(3)%t ) ) & 4259 ALLOCATE( t_wall_v_1(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) ) 4260 READ ( 13 ) t_wall_v_1(3)%t 3492 4261 #endif 3493 4262 … … 3593 4362 INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: usm_par 3594 4363 REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: usm_val 3595 INTEGER(iwp) :: k, l, d, iw, jw, kw, it, ip, ii, ij 4364 INTEGER(iwp) :: k, l, d, iw, jw, kw, it, ip, ii, ij, m 3596 4365 INTEGER(iwp) :: i, j 3597 4366 INTEGER(iwp) :: nz, roof, dirwe, dirsn … … 3733 4502 ENDDO 3734 4503 ENDDO 3735 3736 !-- assign the surface types to local surface array 3737 DO l = startenergy, endenergy 3738 3739 d = surfl(id,l) 3740 kw = surfl(iz,l) 3741 j = surfl(iy,l) 3742 i = surfl(ix,l) 3743 IF ( d == iroof ) THEN 3744 !-- horizontal surface - land or roof 3745 iw = i 3746 jw = j 3747 IF ( usm_par(5,jw,iw) == 0 ) THEN 3748 IF ( zu(kw) >= roof_height_limit ) THEN 3749 isroof_surf(l) = .TRUE. 3750 surface_types(l) = roof_category !< default category for root surface 3751 ELSE 3752 isroof_surf(l) = .FALSE. 3753 surface_types(l) = land_category !< default category for land surface 3754 ENDIF 3755 albedo_surf(l) = -1.0_wp 3756 thickness_wall(l) = -1.0_wp 3757 ELSE 3758 IF ( usm_par(2,jw,iw)==0 ) THEN 3759 isroof_surf(l) = .FALSE. 3760 thickness_wall(l) = -1.0_wp 3761 ELSE 3762 isroof_surf(l) = .TRUE. 3763 thickness_wall(l) = usm_val(2,jw,iw) 3764 ENDIF 3765 surface_types(l) = usm_par(5,jw,iw) 3766 albedo_surf(l) = usm_val(1,jw,iw) 3767 ENDIF 3768 ELSE 3769 SELECT CASE (d) 3770 CASE (iwest) 3771 iw = i 3772 jw = j 3773 ii = 6 3774 ij = 3 3775 CASE (ieast) 3776 iw = i-1 3777 jw = j 3778 ii = 6 3779 ij = 3 3780 CASE (isouth) 3781 iw = i 3782 jw = j 3783 ii = 12 3784 ij = 9 3785 CASE (inorth) 3786 iw = i 3787 jw = j-1 3788 ii = 12 3789 ij = 9 3790 END SELECT 3791 3792 IF ( kw <= usm_par(ii,jw,iw) ) THEN 3793 !-- pedestrant zone 3794 isroof_surf(l) = .FALSE. 3795 IF ( usm_par(ii+1,jw,iw) == 0 ) THEN 3796 surface_types(l) = pedestrant_category !< default category for wall surface in pedestrant zone 3797 albedo_surf(l) = -1.0_wp 3798 thickness_wall(l) = -1.0_wp 3799 ELSE 3800 surface_types(l) = usm_par(ii+1,jw,iw) 3801 albedo_surf(l) = usm_val(ij,jw,iw) 3802 thickness_wall(l) = usm_val(ij+1,jw,iw) 3803 ENDIF 3804 ELSE IF ( kw <= usm_par(ii+2,jw,iw) ) THEN 3805 !-- wall zone 3806 isroof_surf(l) = .FALSE. 3807 IF ( usm_par(ii+3,jw,iw) == 0 ) THEN 3808 surface_types(l) = wall_category !< default category for wall surface 3809 albedo_surf(l) = -1.0_wp 3810 thickness_wall(l) = -1.0_wp 3811 ELSE 3812 surface_types(l) = usm_par(ii+3,jw,iw) 3813 albedo_surf(l) = usm_val(ij+2,jw,iw) 3814 thickness_wall(l) = usm_val(ij+3,jw,iw) 3815 ENDIF 3816 ELSE IF ( kw <= usm_par(ii+4,jw,iw) ) THEN 3817 !-- roof zone 3818 isroof_surf(l) = .TRUE. 3819 IF ( usm_par(ii+5,jw,iw) == 0 ) THEN 3820 surface_types(l) = roof_category !< default category for roof surface 3821 albedo_surf(l) = -1.0_wp 3822 thickness_wall(l) = -1.0_wp 3823 ELSE 3824 surface_types(l) = usm_par(ii+5,jw,iw) 3825 albedo_surf(l) = usm_val(ij+4,jw,iw) 3826 thickness_wall(l) = usm_val(ij+5,jw,iw) 3827 ENDIF 3828 ELSE 3829 !-- something wrong 3830 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 ) 3831 ENDIF 3832 ENDIF 3833 3834 !-- find the type position 3835 it = surface_types(l) 3836 ip = -99999 3837 DO k = 1, n_surface_types 3838 IF ( surface_type_codes(k) == it ) THEN 4504 ! 4505 !-- Assign the surface types to the respective data type. 4506 !-- First, for horizontal upward-facing surfaces. 4507 DO m = 1, surf_usm_h%ns 4508 iw = surf_usm_h%i(m) 4509 jw = surf_usm_h%j(m) 4510 kw = surf_usm_h%k(m) 4511 4512 IF ( usm_par(5,jw,iw) == 0 ) THEN 4513 IF ( zu(kw) >= roof_height_limit ) THEN 4514 surf_usm_h%isroof_surf(m) = .TRUE. 4515 surf_usm_h%surface_types(m) = roof_category !< default category for root surface 4516 ELSE 4517 surf_usm_h%isroof_surf(m) = .FALSE. 4518 surf_usm_h%surface_types(m) = land_category !< default category for land surface 4519 ENDIF 4520 surf_usm_h%albedo_surf(m) = -1.0_wp 4521 surf_usm_h%thickness_wall(m) = -1.0_wp 4522 ELSE 4523 IF ( usm_par(2,jw,iw)==0 ) THEN 4524 surf_usm_h%isroof_surf(m) = .FALSE. 4525 surf_usm_h%thickness_wall(m) = -1.0_wp 4526 ELSE 4527 surf_usm_h%isroof_surf(m) = .TRUE. 4528 surf_usm_h%thickness_wall(m) = usm_val(2,jw,iw) 4529 ENDIF 4530 surf_usm_h%surface_types(m) = usm_par(5,jw,iw) 4531 surf_usm_h%albedo_surf(m) = usm_val(1,jw,iw) 4532 ENDIF 4533 ! 4534 !-- Find the type position 4535 it = surf_usm_h%surface_types(m) 4536 ip = -99999 4537 DO k = 1, n_surface_types 4538 IF ( surface_type_codes(k) == it ) THEN 4539 ip = k 4540 EXIT 4541 ENDIF 4542 ENDDO 4543 IF ( ip == -99999 ) THEN 4544 !-- wall category not found 4545 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, & 4546 ' not found for i,j,k=', iw,jw,kw 4547 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 ) 4548 ENDIF 4549 ! 4550 !-- Albedo 4551 IF ( surf_usm_h%albedo_surf(m) < 0.0_wp ) THEN 4552 surf_usm_h%albedo_surf(m) = surface_params(ialbedo,ip) 4553 ENDIF 4554 ! 4555 !-- emissivity of the wall 4556 surf_usm_h%emiss_surf(m) = surface_params(iemiss,ip) 4557 ! 4558 !-- heat conductivity λS between air and wall ( W mâ2 Kâ1 ) 4559 surf_usm_h%lambda_surf(m) = surface_params(ilambdas,ip) 4560 ! 4561 !-- roughness relative to concrete 4562 surf_usm_h%roughness_wall(m) = surface_params(irough,ip) 4563 ! 4564 !-- Surface skin layer heat capacity (J mâ2 Kâ1 ) 4565 surf_usm_h%c_surface(m) = surface_params(icsurf,ip) 4566 ! 4567 !-- wall material parameters: 4568 !-- thickness of the wall (m) 4569 !-- missing values are replaced by default value for category 4570 IF ( surf_usm_h%thickness_wall(m) <= 0.001_wp ) THEN 4571 surf_usm_h%thickness_wall(m) = surface_params(ithick,ip) 4572 ENDIF 4573 ! 4574 !-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 ) 4575 surf_usm_h%rho_c_wall(:,m) = surface_params(irhoC,ip) 4576 ! 4577 !-- thermal conductivity λH of the wall (W mâ1 Kâ1 ) 4578 surf_usm_h%lambda_h(:,m) = surface_params(ilambdah,ip) 4579 4580 ENDDO 4581 ! 4582 !-- For vertical surface elements ( 0 -- northward-facing, 1 -- southward-facing, 4583 !-- 2 -- eastward-facing, 3 -- westward-facing ) 4584 DO l = 0, 3 4585 DO m = 1, surf_usm_v(l)%ns 4586 i = surf_usm_v(l)%i(m) 4587 j = surf_usm_v(l)%j(m) 4588 kw = surf_usm_v(l)%k(m) 4589 4590 IF ( l == 3 ) THEN ! westward facing 4591 iw = i 4592 jw = j 4593 ii = 6 4594 ij = 3 4595 ELSEIF ( l == 2 ) THEN 4596 iw = i-1 4597 jw = j 4598 ii = 6 4599 ij = 3 4600 ELSEIF ( l == 1 ) THEN 4601 iw = i 4602 jw = j 4603 ii = 12 4604 ij = 9 4605 ELSEIF ( l == 0 ) THEN 4606 iw = i 4607 jw = j-1 4608 ii = 12 4609 ij = 9 4610 ENDIF 4611 4612 IF ( kw <= usm_par(ii,jw,iw) ) THEN 4613 !-- pedestrant zone 4614 IF ( usm_par(ii+1,jw,iw) == 0 ) THEN 4615 surf_usm_v(l)%surface_types(m) = pedestrant_category !< default category for wall surface in pedestrant zone 4616 surf_usm_v(l)%albedo_surf(m) = -1.0_wp 4617 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 4618 ELSE 4619 surf_usm_v(l)%surface_types(m) = usm_par(ii+1,jw,iw) 4620 surf_usm_v(l)%albedo_surf(m) = usm_val(ij,jw,iw) 4621 surf_usm_v(l)%thickness_wall(m) = usm_val(ij+1,jw,iw) 4622 ENDIF 4623 ELSE IF ( kw <= usm_par(ii+2,jw,iw) ) THEN 4624 !-- wall zone 4625 IF ( usm_par(ii+3,jw,iw) == 0 ) THEN 4626 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface 4627 surf_usm_v(l)%albedo_surf(m) = -1.0_wp 4628 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 4629 ELSE 4630 surf_usm_v(l)%surface_types(m) = usm_par(ii+3,jw,iw) 4631 surf_usm_v(l)%albedo_surf(m) = usm_val(ij+2,jw,iw) 4632 surf_usm_v(l)%thickness_wall(m) = usm_val(ij+3,jw,iw) 4633 ENDIF 4634 ELSE IF ( kw <= usm_par(ii+4,jw,iw) ) THEN 4635 !-- roof zone 4636 IF ( usm_par(ii+5,jw,iw) == 0 ) THEN 4637 surf_usm_v(l)%surface_types(m) = roof_category !< default category for roof surface 4638 surf_usm_v(l)%albedo_surf(m) = -1.0_wp 4639 surf_usm_v(l)%thickness_wall(m) = -1.0_wp 4640 ELSE 4641 surf_usm_v(l)%surface_types(m) = usm_par(ii+5,jw,iw) 4642 surf_usm_v(l)%albedo_surf(m) = usm_val(ij+4,jw,iw) 4643 surf_usm_v(l)%thickness_wall(m) = usm_val(ij+5,jw,iw) 4644 ENDIF 4645 ELSE 4646 !-- something wrong 4647 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 ) 4648 ENDIF 4649 4650 ! 4651 !-- Find the type position 4652 it = surf_usm_v(l)%surface_types(m) 4653 ip = -99999 4654 DO k = 1, n_surface_types 4655 IF ( surface_type_codes(k) == it ) THEN 3839 4656 ip = k 3840 4657 EXIT 3841 ENDIF 3842 ENDDO 3843 IF ( ip == -99999 ) THEN 3844 !-- wall category not found 3845 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, ' not found for i,j,k=', iw,jw,kw 3846 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 ) 3847 ENDIF 3848 3849 !-- Fill out the parameters of the wall 3850 !-- wall surface: 3851 3852 !-- albedo 3853 IF ( albedo_surf(l) < 0.0_wp ) THEN 3854 albedo_surf(l) = surface_params(ialbedo, ip) 3855 ENDIF 3856 3857 !-- emissivity of the wall 3858 emiss_surf(l) = surface_params(iemiss, ip) 3859 3860 !-- heat conductivity λS between air and wall ( W mâ2 Kâ1 ) 3861 lambda_surf(l) = surface_params(ilambdas, ip) 3862 3863 !-- roughness relative to concrete 3864 roughness_wall(l) = surface_params(irough, ip) 3865 3866 !-- Surface skin layer heat capacity (J mâ2 Kâ1 ) 3867 c_surface(l) = surface_params(icsurf, ip) 3868 3869 !-- wall material parameters: 3870 3871 !-- thickness of the wall (m) 3872 !-- missing values are replaced by default value for category 3873 IF ( thickness_wall(l) <= 0.001_wp ) THEN 3874 thickness_wall(l) = surface_params(ithick, ip) 3875 ENDIF 3876 3877 !-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 ) 3878 rho_c_wall(:,l) = surface_params(irhoC, ip) 3879 3880 !-- thermal conductivity λH of the wall (W mâ1 Kâ1 ) 3881 lambda_h(:,l) = surface_params(ilambdah, ip) 3882 4658 ENDIF 4659 ENDDO 4660 IF ( ip == -99999 ) THEN 4661 !-- wall category not found 4662 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, & 4663 ' not found for i,j,k=', iw,jw,kw 4664 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 ) 4665 ENDIF 4666 ! 4667 !-- Albedo 4668 IF ( surf_usm_v(l)%albedo_surf(m) < 0.0_wp ) THEN 4669 surf_usm_v(l)%albedo_surf(m) = surface_params(ialbedo,ip) 4670 ENDIF 4671 ! 4672 !-- emissivity of the wall 4673 surf_usm_v(l)%emiss_surf(m) = surface_params(iemiss,ip) 4674 ! 4675 !-- heat conductivity λS between air and wall ( W mâ2 Kâ1 ) 4676 surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip) 4677 ! 4678 !-- roughness relative to concrete 4679 surf_usm_v(l)%roughness_wall(m) = surface_params(irough,ip) 4680 ! 4681 !-- Surface skin layer heat capacity (J mâ2 Kâ1 ) 4682 surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip) 4683 ! 4684 !-- wall material parameters: 4685 !-- thickness of the wall (m) 4686 !-- missing values are replaced by default value for category 4687 IF ( surf_usm_v(l)%thickness_wall(m) <= 0.001_wp ) THEN 4688 surf_usm_v(l)%thickness_wall(m) = surface_params(ithick,ip) 4689 ENDIF 4690 ! 4691 !-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 ) 4692 surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip) 4693 ! 4694 !-- thermal conductivity λH of the wall (W mâ1 Kâ1 ) 4695 surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip) 4696 4697 ENDDO 3883 4698 ENDDO 3884 4699 … … 3900 4715 IMPLICIT NONE 3901 4716 3902 INTEGER(iwp) :: i, j, k, l, d 4717 INTEGER(iwp) :: i, j, k, l, d, m !< running indices 3903 4718 3904 4719 REAL(wp) :: pt1 !< temperature at first grid box adjacent to surface … … 3924 4739 3925 4740 exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp !< Exner function 4741 ! 4742 !-- First, treat horizontal surface elements 4743 4744 DO m = 1, surf_usm_h%ns 4745 ! 4746 !-- Get indices of respective grid point 4747 i = surf_usm_h%i(m) 4748 j = surf_usm_h%j(m) 4749 k = surf_usm_h%k(m) 4750 ! 4751 !-- TODO - how to calculate lambda_surface for horizontal surfaces 4752 !-- (lambda_surface is set according to stratification in land surface model) 4753 !-- MS: ??? 4754 IF ( surf_usm_h%ol(m) >= 0.0_wp ) THEN 4755 lambda_surface = surf_usm_h%lambda_surf(m) 4756 ELSE 4757 lambda_surface = surf_usm_h%lambda_surf(m) 4758 ENDIF 3926 4759 3927 !-- 3928 DO l = startenergy, endenergy 3929 !-- Calculate frequently used parameters 3930 d = surfl(id,l) 3931 k = surfl(iz,l) 3932 j = surfl(iy,l) 3933 i = surfl(ix,l) 3934 3935 !-- TODO - how to calculate lambda_surface for horizontal surfaces 3936 !-- (lambda_surface is set according to stratification in land surface model) 3937 IF ( ol(j,i) >= 0.0_wp ) THEN 3938 lambda_surface = lambda_surf(l) 3939 ELSE 3940 lambda_surface = lambda_surf(l) 3941 ENDIF 4760 pt1 = pt(k,j,i) 4761 ! 4762 !-- calculate rho * cp coefficient at surface layer 4763 rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) ) 4764 ! 4765 !-- Calculate aerodyamic resistance. 4766 !-- Calculation for horizontal surfaces follows LSM formulation 4767 !-- pt, us, ts are not available for the prognostic time step, 4768 !-- data from the last time step is used here. 4769 4770 r_a = ( pt1 - t_surf_h(m) / exn(k) ) / & 4771 ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp ) 4772 4773 !-- make sure that the resistance does not drop to zero 4774 IF ( ABS(r_a) < 1.0E-10_wp ) r_a = 1.0E-10_wp 4775 4776 !-- the parameterization is developed originally for larger scales 4777 !-- (compare with remark in TUF-3D) 4778 !-- our first experiences show that the parameterization underestimates 4779 !-- r_a in meter resolution. 4780 !-- temporary solution - multiplication by magic constant :-(. 4781 r_a = r_a * ra_horiz_coef 4782 4783 !-- factor for shf_eb 4784 f_shf = rho_cp / r_a 4785 4786 !-- add LW up so that it can be removed in prognostic equation 4787 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_in_sw(m) - & 4788 surf_usm_h%rad_out_sw(m) + & 4789 surf_usm_h%rad_in_lw(m) - & 4790 surf_usm_h%rad_out_lw(m) 4791 4792 !-- numerator of the prognostic equation 4793 coef_1 = surf_usm_h%rad_net_l(m) + & 4794 ( 3.0_wp + 1.0_wp ) * surf_usm_h%emiss_surf(m) * sigma_sb * & 4795 t_surf_h(m) ** 4 + & 4796 f_shf * pt1 + & 4797 lambda_surface * t_wall_h(nzb_wall,m) 4798 4799 !-- denominator of the prognostic equation 4800 coef_2 = 4.0_wp * surf_usm_h%emiss_surf(m) * sigma_sb * & 4801 t_surf_h(m) ** 3 & 4802 + lambda_surface + f_shf / exn(k) 4803 4804 !-- implicit solution when the surface layer has no heat capacity, 4805 !-- otherwise use RK3 scheme. 4806 t_surf_h_p(m) = ( coef_1 * dt_3d * tsc(2) + & 4807 surf_usm_h%c_surface(m) * t_surf_h(m) ) / & 4808 ( surf_usm_h%c_surface(m) + coef_2 * dt_3d * tsc(2) ) 4809 4810 !-- add RK3 term 4811 t_surf_h_p(m) = t_surf_h_p(m) + dt_3d * tsc(3) * & 4812 surf_usm_h%tt_surface_m(m) 3942 4813 3943 pt1 = pt(k,j,i) 3944 3945 !-- calculate rho * cp coefficient at surface layer 3946 rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) ) 3947 3948 !-- calculate aerodyamic resistance. 3949 IF ( d == iroof ) THEN 3950 !-- calculation for horizontal surfaces follows LSM formulation 3951 !-- pt, us, ts are not available for the prognostic time step, 3952 !-- data from the last time step is used here. 4814 !-- calculate true tendency 4815 stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) * & 4816 surf_usm_h%tt_surface_m(m)) / ( dt_3d * tsc(2) ) 4817 4818 !-- calculate t_surf tendencies for the next Runge-Kutta step 4819 IF ( timestep_scheme(1:5) == 'runge' ) THEN 4820 IF ( intermediate_timestep_count == 1 ) THEN 4821 surf_usm_h%tt_surface_m(m) = stend 4822 ELSEIF ( intermediate_timestep_count < & 4823 intermediate_timestep_count_max ) THEN 4824 surf_usm_h%tt_surface_m(m) = -9.5625_wp * stend + & 4825 5.3125_wp * surf_usm_h%tt_surface_m(m) 4826 ENDIF 4827 ENDIF 4828 4829 !-- in case of fast changes in the skin temperature, it is required to 4830 !-- update the radiative fluxes in order to keep the solution stable 4831 IF ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) THEN 4832 force_radiation_call_l = .TRUE. 4833 ENDIF 4834 ! 4835 !-- for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf. 4836 !-- there is no equivalent surface gridpoint for vertical surfaces. 4837 !-- pt(k,j,i) is calculated for all directions in diffusion_s 4838 !-- using surface and wall heat fluxes 4839 pt(k-1,j,i) = t_surf_h_p(m) / exn(k) ! not for vertical surfaces 4840 4841 !-- calculate fluxes 4842 !-- rad_net_l is never used! 4843 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) + & 4844 3.0_wp * sigma_sb * & 4845 t_surf_h(m)**4 - 4.0_wp * sigma_sb * & 4846 t_surf_h(m)**3 * t_surf_h_p(m) 4847 surf_usm_h%wghf_eb(m) = lambda_surface * & 4848 ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) ) 4849 ! 4850 !-- ground/wall/roof surface heat flux 4851 surf_usm_h%wshf_eb(m) = - f_shf * ( pt1 - t_surf_h_p(m) ) 4852 ! 4853 !-- store kinematic surface heat fluxes for utilization in other processes 4854 !-- diffusion_s, surface_layer_fluxes,... 4855 surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / rho_cp 4856 4857 ENDDO 4858 ! 4859 !-- Now, treat vertical surface elements 4860 DO l = 0, 3 4861 DO m = 1, surf_usm_v(l)%ns 4862 ! 4863 !-- Get indices of respective grid point 4864 i = surf_usm_v(l)%i(m) 4865 j = surf_usm_v(l)%j(m) 4866 k = surf_usm_v(l)%k(m) 4867 ! 4868 !-- TODO - how to calculate lambda_surface for horizontal (??? do you mean verical ???) surfaces 4869 !-- (lambda_surface is set according to stratification in land surface model). 4870 !-- Please note, for vertical surfaces no ol is defined, since 4871 !-- stratification is not considered in this case. 4872 lambda_surface = surf_usm_v(l)%lambda_surf(m) 4873 4874 pt1 = pt(k,j,i) 4875 ! 4876 !-- calculate rho * cp coefficient at surface layer 4877 rho_cp = cp * hyp(k) / ( r_d * pt1 * exn(k) ) 4878 4879 !-- Calculation of r_a for vertical surfaces 4880 !-- 4881 !-- heat transfer coefficient for forced convection along vertical walls 4882 !-- follows formulation in TUF3d model (Krayenhoff & Voogt, 2006) 4883 !-- 4884 !-- H = httc (Tsfc - Tair) 4885 !-- httc = rw * (11.8 + 4.2 * Ueff) - 4.0 4886 !-- 4887 !-- rw: wall patch roughness relative to 1.0 for concrete 4888 !-- Ueff: effective wind speed 4889 !-- - 4.0 is a reduction of Rowley et al (1930) formulation based on 4890 !-- Cole and Sturrock (1977) 4891 !-- 4892 !-- Ucan: Canyon wind speed 4893 !-- wstar: convective velocity 4894 !-- Qs: surface heat flux 4895 !-- zH: height of the convective layer 4896 !-- wstar = (g/Tcan*Qs*zH)**(1./3.) 3953 4897 3954 r_a = (pt1 - t_surf(l)/exn(k)) / (ts(j,i) * us(j,i) + 1.0E-10_wp) 4898 !-- Effective velocity components must always 4899 !-- be defined at scalar grid point. The wall normal component is 4900 !-- obtained by simple linear interpolation. ( An alternative would 4901 !-- be an logarithmic interpolation. ) 4902 u1 = ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp 4903 v1 = ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp 4904 w1 = ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp 3955 4905 3956 !-- make sure that the resistance does not drop to zero 3957 IF ( ABS(r_a) < 1.0E-10_wp ) r_a = 1.0E-10_wp 3958 3959 !-- the parameterization is developed originally for larger scales 3960 !-- (compare with remark in TUF-3D) 3961 !-- our first experiences show that the parameterization underestimates 3962 !-- r_a in meter resolution. 3963 !-- temporary solution - multiplication by magic constant :-(. 3964 r_a = r_a * ra_horiz_coef 3965 3966 !-- factor for shf_eb 3967 f_shf = rho_cp / r_a 3968 ELSE 3969 !-- calculation of r_a for vertical surfaces 3970 !-- 3971 !-- heat transfer coefficient for forced convection along vertical walls 3972 !-- follows formulation in TUF3d model (Krayenhoff & Voogt, 2006) 3973 !-- 3974 !-- H = httc (Tsfc - Tair) 3975 !-- httc = rw * (11.8 + 4.2 * Ueff) - 4.0 3976 !-- 3977 !-- rw: wall patch roughness relative to 1.0 for concrete 3978 !-- Ueff: effective wind speed 3979 !-- - 4.0 is a reduction of Rowley et al (1930) formulation based on 3980 !-- Cole and Sturrock (1977) 3981 !-- 3982 !-- Ucan: Canyon wind speed 3983 !-- wstar: convective velocity 3984 !-- Qs: surface heat flux 3985 !-- zH: height of the convective layer 3986 !-- wstar = (g/Tcan*Qs*zH)**(1./3.) 3987 3988 !-- staggered grid needs to be taken into consideration 3989 IF ( d == inorth ) THEN 3990 u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp 3991 v1 = v(k,j+1,i) 3992 ELSE IF ( d == isouth ) THEN 3993 u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp 3994 v1 = v(k,j,i) 3995 ELSE IF ( d == ieast ) THEN 3996 u1 = u(k,j,i+1) 3997 v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp 3998 ELSE IF ( d == iwest ) THEN 3999 u1 = u(k,j,i) 4000 v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp 4001 ELSE 4002 STOP 4003 ENDIF 4004 w1 = (w(k,j,i)+w(k-1,j,i))*0.5_wp 4005 4006 Ueff = SQRT(u1**2 + v1**2 + w1**2) 4007 httc = roughness_wall(l) * (11.8 + 4.2 * Ueff) - 4.0 4008 f_shf = httc 4009 ENDIF 4010 4011 !-- add LW up so that it can be removed in prognostic equation 4012 rad_net_l(l) = surfinsw(l) - surfoutsw(l) + surfinlw(l) - surfoutlw(l) 4013 4014 !-- numerator of the prognostic equation 4015 coef_1 = rad_net_l(l) + & ! coef +1 corresponds to -lwout included in calculation of radnet_l 4016 (3.0_wp+1.0_wp) * emiss_surf(l) * sigma_sb * t_surf(l) ** 4 + & 4017 f_shf * pt1 + & 4018 lambda_surface * t_wall(nzb_wall,l) 4019 4020 !-- denominator of the prognostic equation 4021 coef_2 = 4.0_wp * emiss_surf(l) * sigma_sb * t_surf(l) ** 3 & 4022 + lambda_surface + f_shf / exn(k) 4023 4024 !-- implicit solution when the surface layer has no heat capacity, 4025 !-- otherwise use RK3 scheme. 4026 t_surf_p(l) = ( coef_1 * dt_3d * tsc(2) + c_surface(l) * t_surf(l) ) / & 4027 ( c_surface(l) + coef_2 * dt_3d * tsc(2) ) 4028 4029 !-- add RK3 term 4030 t_surf_p(l) = t_surf_p(l) + dt_3d * tsc(3) * tt_surface_m(l) 4906 Ueff = SQRT( u1**2 + v1**2 + w1**2 ) 4907 httc = surf_usm_v(l)%roughness_wall(m) * & 4908 ( 11.8_wp + 4.2_wp * Ueff ) - 4.0_wp 4909 f_shf = httc 4910 4911 !-- add LW up so that it can be removed in prognostic equation 4912 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_in_sw(m) - & 4913 surf_usm_v(l)%rad_out_sw(m) + & 4914 surf_usm_v(l)%rad_in_lw(m) - & 4915 surf_usm_v(l)%rad_out_lw(m) 4916 4917 !-- numerator of the prognostic equation 4918 coef_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout included in calculation of radnet_l 4919 ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emiss_surf(m) * sigma_sb * & 4920 t_surf_v(l)%t(m) ** 4 + & 4921 f_shf * pt1 + & 4922 lambda_surface * t_wall_v(l)%t(nzb_wall,m) 4923 4924 !-- denominator of the prognostic equation 4925 coef_2 = 4.0_wp * surf_usm_v(l)%emiss_surf(m) * sigma_sb * & 4926 t_surf_v(l)%t(m) ** 3 & 4927 + lambda_surface + f_shf / exn(k) 4928 4929 !-- implicit solution when the surface layer has no heat capacity, 4930 !-- otherwise use RK3 scheme. 4931 t_surf_v_p(l)%t(m) = ( coef_1 * dt_3d * tsc(2) + & 4932 surf_usm_v(l)%c_surface(m) * t_surf_v(l)%t(m) ) / & 4933 ( surf_usm_v(l)%c_surface(m) + coef_2 * dt_3d * tsc(2) ) 4934 4935 !-- add RK3 term 4936 t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) * & 4937 surf_usm_v(l)%tt_surface_m(m) 4031 4938 4032 !-- calculate true tendency 4033 stend = (t_surf_p(l) - t_surf(l) - dt_3d * tsc(3) * tt_surface_m(l)) / (dt_3d * tsc(2)) 4034 4035 !-- calculate t_surf tendencies for the next Runge-Kutta step 4036 IF ( timestep_scheme(1:5) == 'runge' ) THEN 4037 IF ( intermediate_timestep_count == 1 ) THEN 4038 tt_surface_m(l) = stend 4039 ELSEIF ( intermediate_timestep_count < & 4040 intermediate_timestep_count_max ) THEN 4041 tt_surface_m(l) = -9.5625_wp * stend + 5.3125_wp & 4042 * tt_surface_m(l) 4043 ENDIF 4044 ENDIF 4045 4046 !-- in case of fast changes in the skin temperature, it is required to 4047 !-- update the radiative fluxes in order to keep the solution stable 4048 IF ( ABS( t_surf_p(l) - t_surf(l) ) > 1.0_wp ) THEN 4049 force_radiation_call_l = .TRUE. 4050 ENDIF 4051 4052 !-- for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf. 4053 !-- there is no equivalent surface gridpoint for vertical surfaces. 4054 !-- pt(k,j,i) is calculated for all directions in diffusion_s 4055 !-- using surface and wall heat fluxes 4056 IF ( d == iroof ) THEN 4057 pt(nzb_s_inner(j,i),j,i) = t_surf_p(l) / exn(k) 4058 ENDIF 4059 4060 !-- calculate fluxes 4061 !-- rad_net_l is never used! 4062 rad_net_l(l) = rad_net_l(l) + 3.0_wp * sigma_sb & 4063 * t_surf(l)**4 - 4.0_wp * sigma_sb & 4064 * t_surf(l)**3 * t_surf_p(l) 4065 wghf_eb(l) = lambda_surface * (t_surf_p(l) - t_wall(nzb_wall,l)) 4066 4067 !-- ground/wall/roof surface heat flux 4068 wshf_eb(l) = - f_shf * ( pt1 - t_surf_p(l) ) 4069 4070 !-- store kinematic surface heat fluxes for utilization in other processes 4071 !-- diffusion_s, surface_layer_fluxes,... 4072 IF ( d == iroof ) THEN 4073 !-- shf is used in diffusion_s and also 4074 !-- for calculation of surface layer fluxes 4075 !-- update for horizontal surfaces 4076 shf(j,i) = wshf_eb(l) / rho_cp 4077 ELSE 4078 !-- surface heat flux for vertical surfaces 4079 !-- used in diffusion_s 4080 wshf(l) = wshf_eb(l) / rho_cp 4081 ENDIF 4939 !-- calculate true tendency 4940 stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *& 4941 surf_usm_v(l)%tt_surface_m(m) ) / ( dt_3d * tsc(2) ) 4942 4943 !-- calculate t_surf tendencies for the next Runge-Kutta step 4944 IF ( timestep_scheme(1:5) == 'runge' ) THEN 4945 IF ( intermediate_timestep_count == 1 ) THEN 4946 surf_usm_v(l)%tt_surface_m(m) = stend 4947 ELSEIF ( intermediate_timestep_count < & 4948 intermediate_timestep_count_max ) THEN 4949 surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend + & 4950 5.3125_wp * surf_usm_h%tt_surface_m(m) 4951 ENDIF 4952 ENDIF 4953 4954 !-- in case of fast changes in the skin temperature, it is required to 4955 !-- update the radiative fluxes in order to keep the solution stable 4956 IF ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) THEN 4957 force_radiation_call_l = .TRUE. 4958 ENDIF 4959 4960 !-- calculate fluxes 4961 !-- rad_net_l is never used! 4962 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_net_l(m) + & 4963 3.0_wp * sigma_sb * & 4964 t_surf_v(l)%t(m)**4 - 4.0_wp * sigma_sb * & 4965 t_surf_v(l)%t(m)**3 * t_surf_v_p(l)%t(m) 4966 4967 surf_usm_v(l)%wghf_eb(m) = lambda_surface * & 4968 ( t_surf_v_p(l)%t(m) - t_wall_v(l)%t(nzb_wall,m) ) 4969 4970 !-- ground/wall/roof surface heat flux 4971 surf_usm_v(l)%wshf_eb(m) = - f_shf * ( pt1 - t_surf_v_p(l)%t(m) ) 4972 4973 ! 4974 !-- store kinematic surface heat fluxes for utilization in other processes 4975 !-- diffusion_s, surface_layer_fluxes,... 4976 surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / rho_cp 4977 4978 ENDDO 4082 4979 4083 4980 ENDDO 4084 4085 4981 ! 4982 !-- Add-up anthropogenic heat, for now only at upward-facing surfaces 4086 4983 IF ( usm_anthropogenic_heat .AND. & 4087 4984 intermediate_timestep_count == intermediate_timestep_count_max ) THEN 4088 !-- 4089 !-- 4090 !-- to the first layer, generalization would be worth4985 !-- application of the additional anthropogenic heat sources 4986 !-- we considere the traffic for now so all heat is absorbed 4987 !-- to the first layer, generalization would be worth. 4091 4988 4092 !-- calculation of actual profile coefficient 4093 !-- ??? check time_since_reference_point ??? 4094 dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp) 4095 dhour = INT(dtime/3600.0_wp) 4096 !-- linear interpolation of coeficient 4097 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1) 4098 DO i = nxl, nxr 4099 DO j = nys, nyn 4100 IF ( aheat(j,i) > 0.0_wp ) THEN 4101 !-- TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d 4102 !-- given to anthropogenic heat aheat*acoef (W*m-2) 4103 !-- k = nzb_s_inner(j,i)+1 4104 !-- pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz) 4105 !-- Instead of this, we can adjust shf in case AH only at surface 4106 shf(j,i) = shf(j,i) + aheat(j,i)*acoef * ddx * ddy / rho_cp 4107 ENDIF 4108 ENDDO 4109 ENDDO 4989 !-- calculation of actual profile coefficient 4990 !-- ??? check time_since_reference_point ??? 4991 dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp) 4992 dhour = INT(dtime/3600.0_wp) 4993 !-- linear interpolation of coeficient 4994 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1) 4995 4996 DO m = 1, surf_usm_h%ns 4997 ! 4998 !-- Get indices of respective grid point 4999 i = surf_usm_h%i(m) 5000 j = surf_usm_h%j(m) 5001 k = surf_usm_h%k(m) 5002 5003 IF ( aheat(j,i) > 0.0_wp ) THEN 5004 !-- TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d 5005 !-- given to anthropogenic heat aheat*acoef (W*m-2) 5006 !-- k = nzb_s_inner(j,i)+1 5007 !-- pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz) 5008 !-- Instead of this, we can adjust shf in case AH only at surface 5009 surf_usm_h%shf(m) = surf_usm_h%shf(m) + & 5010 aheat(j,i) * acoef * ddx * ddy / rho_cp 5011 ENDIF 5012 ENDDO 5013 4110 5014 ENDIF 4111 5015 … … 4113 5017 !-- get the borders from neighbours 4114 5018 CALL exchange_horiz( pt, nbgp ) 4115 CALL exchange_horiz_2d( shf ) 4116 4117 4118 !-- calculation of force_radiation_call: 4119 !-- Make logical OR for all processes. 4120 !-- Force radiation call if at least one processor forces it. 4121 IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 ) & 4122 THEN 5019 5020 5021 !-- calculation of force_radiation_call: 5022 !-- Make logical OR for all processes. 5023 !-- Force radiation call if at least one processor forces it. 5024 IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )& 5025 THEN 4123 5026 #if defined( __parallel ) 4124 5027 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 4125 CALL mpi_allreduce( force_radiation_call_l, force_radiation_call,&4126 5028 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 5029 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 4127 5030 #else 4128 5031 force_radiation_call = force_radiation_call_l … … 4130 5033 force_radiation_call_l = .FALSE. 4131 5034 ENDIF 4132 4133 5035 4134 5036 END SUBROUTINE usm_surface_energy_balance … … 4149 5051 4150 5052 #if defined( __nopointer ) 4151 t_surf = t_surf_p 4152 t_wall = t_wall_p 5053 t_surf_h = t_surf_h_p 5054 t_wall_h = t_wall_h_p 5055 t_surf_v = t_surf_v_p 5056 t_wall_v = t_wall_v_p 4153 5057 #else 4154 5058 SELECT CASE ( mod_count ) 4155 5059 CASE ( 0 ) 4156 t_surf => t_surf_1; t_surf_p => t_surf_2 4157 t_wall => t_wall_1; t_wall_p => t_wall_2 5060 ! 5061 !-- Horizontal surfaces 5062 t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2 5063 t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2 5064 ! 5065 !-- Vertical surfaces 5066 t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2 5067 t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2 4158 5068 CASE ( 1 ) 4159 t_surf => t_surf_2; t_surf_p => t_surf_1 4160 t_wall => t_wall_2; t_wall_p => t_wall_1 5069 ! 5070 !-- Horizontal surfaces 5071 t_surf_h => t_surf_h_2; t_surf_h_p => t_surf_h_1 5072 t_wall_h => t_wall_h_2; t_wall_h_p => t_wall_h_1 5073 ! 5074 !-- Vertical surfaces 5075 t_surf_v => t_surf_v_2; t_surf_v_p => t_surf_v_1 5076 t_wall_v => t_wall_v_2; t_wall_v_p => t_wall_v_1 4161 5077 END SELECT 4162 5078 #endif 4163 5079 4164 5080 END SUBROUTINE usm_swap_timelevel 4165 4166 4167 !------------------------------------------------------------------------------!4168 ! Description:4169 ! ------------4170 !4171 !> This function applies the kinematic wall heat fluxes4172 !> for walls in four directions for all gridboxes in urban layer.4173 !> It is called out from subroutine prognostic_equations.4174 !> TODO Compare performance with cycle runnig l=startwall,endwall...4175 !------------------------------------------------------------------------------!4176 SUBROUTINE usm_wall_heat_flux4177 4178 IMPLICIT NONE4179 4180 INTEGER(iwp) :: i,j,k,d,l !< running indices4181 4182 DO l = startenergy, endenergy4183 j = surfl(iy,l)4184 i = surfl(ix,l)4185 k = surfl(iz,l)4186 d = surfl(id,l)4187 tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)4188 ENDDO4189 4190 END SUBROUTINE usm_wall_heat_flux4191 4192 4193 !------------------------------------------------------------------------------!4194 ! Description:4195 ! ------------4196 !4197 !> This function applies the kinematic wall heat fluxes4198 !> for walls in four directions around the gridbox i,j.4199 !> It is called out from subroutine prognostic_equations.4200 !------------------------------------------------------------------------------!4201 SUBROUTINE usm_wall_heat_flux_ij(i,j)4202 4203 IMPLICIT NONE4204 4205 INTEGER(iwp), INTENT(in) :: i,j !< indices of grid box4206 INTEGER(iwp) :: ii,jj,k,d,l4207 4208 DO l = startenergy, endenergy4209 jj = surfl(iy,l)4210 ii = surfl(ix,l)4211 IF ( ii == i .AND. jj == j ) THEN4212 k = surfl(iz,l)4213 IF ( k >= nzb_s_inner(j,i)+1 .AND. k <= nzb_s_outer(j,i) ) THEN4214 d = surfl(id,l)4215 IF ( d >= 1 .and. d <= 4 ) THEN4216 tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)4217 ENDIF4218 ENDIF4219 ENDIF4220 ENDDO4221 4222 END SUBROUTINE usm_wall_heat_flux_ij4223 4224 5081 4225 5082 !------------------------------------------------------------------------------! … … 4242 5099 DO i = 0, io_blocks-1 4243 5100 IF ( i == io_group ) THEN 4244 WRITE ( 14 ) 't_surf ' 5101 5102 WRITE ( 14 ) 't_surf_h ' 4245 5103 #if defined( __nopointer ) 4246 WRITE ( 14 ) t_surf 5104 WRITE ( 14 ) t_surf_h 4247 5105 #else 4248 WRITE ( 14 ) t_surf_ 15106 WRITE ( 14 ) t_surf_h_1 4249 5107 #endif 4250 WRITE ( 14 ) 't_ wall'5108 WRITE ( 14 ) 't_surf_v(0) ' 4251 5109 #if defined( __nopointer ) 4252 WRITE ( 14 ) t_ wall5110 WRITE ( 14 ) t_surf_v(0)%t 4253 5111 #else 4254 WRITE ( 14 ) t_wall_1 5112 WRITE ( 14 ) t_surf_v_1(0)%t 5113 #endif 5114 WRITE ( 14 ) 't_surf_v(1) ' 5115 #if defined( __nopointer ) 5116 WRITE ( 14 ) t_surf_v(1)%t 5117 #else 5118 WRITE ( 14 ) t_surf_v_1(1)%t 5119 #endif 5120 WRITE ( 14 ) 't_surf_v(2) ' 5121 #if defined( __nopointer ) 5122 WRITE ( 14 ) t_surf_v(2)%t 5123 #else 5124 WRITE ( 14 ) t_surf_v_1(2)%t 5125 #endif 5126 WRITE ( 14 ) 't_surf_v(3) ' 5127 #if defined( __nopointer ) 5128 WRITE ( 14 ) t_surf_v(3)%t 5129 #else 5130 WRITE ( 14 ) t_surf_v_1(3)%t 5131 #endif 5132 WRITE ( 14 ) 't_wall_h ' 5133 #if defined( __nopointer ) 5134 WRITE ( 14 ) t_wall_h 5135 #else 5136 WRITE ( 14 ) t_wall_h_1 5137 #endif 5138 WRITE ( 14 ) 't_wall_v(0) ' 5139 #if defined( __nopointer ) 5140 WRITE ( 14 ) t_wall_v(0)%t 5141 #else 5142 WRITE ( 14 ) t_wall_v_1(0)%t 5143 #endif 5144 WRITE ( 14 ) 't_wall_v(1) ' 5145 #if defined( __nopointer ) 5146 WRITE ( 14 ) t_wall_v(1)%t 5147 #else 5148 WRITE ( 14 ) t_wall_v_1(1)%t 5149 #endif 5150 WRITE ( 14 ) 't_wall_v(2) ' 5151 #if defined( __nopointer ) 5152 WRITE ( 14 ) t_wall_v(2)%t 5153 #else 5154 WRITE ( 14 ) t_wall_v_1(2)%t 5155 #endif 5156 WRITE ( 14 ) 't_wall_v(3) ' 5157 #if defined( __nopointer ) 5158 WRITE ( 14 ) t_wall_v(3)%t 5159 #else 5160 WRITE ( 14 ) t_wall_v_1(3)%t 4255 5161 #endif 4256 5162 WRITE ( 14 ) '*** end usm *** '
Note: See TracChangeset
for help on using the changeset viewer.
