Changeset 4653
- Timestamp:
- Aug 27, 2020 8:54:43 AM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
init_grid.f90 (modified) (2 diffs)
-
radiation_model_mod.f90 (modified) (112 diffs)
-
urban_surface_mod.f90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r4648 r4653 2411 2411 !-- Pre-calculate topography top indices (former get_topography_top_index 2412 2412 !-- function) 2413 ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0: 4) )2413 ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:5) ) 2414 2414 ! 2415 2415 !-- Uppermost topography index on scalar grid … … 2437 2437 topo_top_ind(:,:,4) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,:,:), ibit ) ), DIM=1 ) & 2438 2438 - 1 2439 ! 2440 !-- Uppermost topography index including full-3D geometry 2441 ibit = 12 2442 DO k = nzb, nzt+1 2443 WHERE( BTEST( wall_flags_total_0(k,:,:), ibit ) ) topo_top_ind(:,:,5) = k 2444 ENDDO 2439 2445 2440 2446 END SUBROUTINE set_topo_flags -
palm/trunk/SOURCE/radiation_model_mod.f90
r4644 r4653 28 28 ! ----------------- 29 29 ! $Id$ 30 ! RTM version 4.0: 31 ! - Support full-3D geometry: 32 ! * Downward faces can only be default type (temporary PALM limitation) 33 ! * rad_angular_discretization = .FALSE. is not supported for full-3D cases 34 ! - Order faces in surf(l) by column 35 ! - Remove separate direction indices for urban and land faces 36 ! - Bugfix: skyvft output 37 ! 38 ! 4644 2020-08-22 15:05:46Z pavelkrc 30 39 ! Bugfix: remove invalid reordering of vertical surfaces in RTM 31 40 ! … … 426 435 427 436 USE surface_mod, & 428 ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, 437 ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, & 429 438 surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v, & 430 439 vertical_surfaces_exist … … 711 720 INTEGER(iwp), PARAMETER :: iy = 3 !< position of j-index in surfl and surf 712 721 INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf 713 INTEGER(iwp), PARAMETER :: im = 5 !< position of surface m-index in surfl and surf 714 INTEGER(iwp), PARAMETER :: nidx_surf = 5 !< number of indices in surfl and surf 715 716 INTEGER(iwp), PARAMETER :: nsurf_type = 10 !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1 717 718 INTEGER(iwp), PARAMETER :: iup_u = 0 !< 0 - index of urban upward surface (ground or roof) 719 INTEGER(iwp), PARAMETER :: idown_u = 1 !< 1 - index of urban downward surface (overhanging) 720 INTEGER(iwp), PARAMETER :: inorth_u = 2 !< 2 - index of urban northward facing wall 721 INTEGER(iwp), PARAMETER :: isouth_u = 3 !< 3 - index of urban southward facing wall 722 INTEGER(iwp), PARAMETER :: ieast_u = 4 !< 4 - index of urban eastward facing wall 723 INTEGER(iwp), PARAMETER :: iwest_u = 5 !< 5 - index of urban westward facing wall 724 725 INTEGER(iwp), PARAMETER :: iup_l = 6 !< 6 - index of land upward surface (ground or roof) 726 INTEGER(iwp), PARAMETER :: inorth_l = 7 !< 7 - index of land northward facing wall 727 INTEGER(iwp), PARAMETER :: isouth_l = 8 !< 8 - index of land southward facing wall 728 INTEGER(iwp), PARAMETER :: ieast_l = 9 !< 9 - index of land eastward facing wall 729 INTEGER(iwp), PARAMETER :: iwest_l = 10 !< 10- index of land westward facing wall 730 731 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/) !< surface normal direction x indices 732 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices 733 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices 722 INTEGER(iwp), PARAMETER :: nidx_surf = 4 !< number of indices in surfl and surf 723 724 INTEGER(iwp), PARAMETER :: nsurf_type = 5 !< number of surf types = surface directions 725 726 INTEGER(iwp), PARAMETER :: iup = 0 !< 0 - index of upward surface (ground or roof) 727 INTEGER(iwp), PARAMETER :: idown = 1 !< 1 - index of downward surface (overhanging) 728 INTEGER(iwp), PARAMETER :: inorth = 2 !< 2 - index of northward facing wall 729 INTEGER(iwp), PARAMETER :: isouth = 3 !< 3 - index of southward facing wall 730 INTEGER(iwp), PARAMETER :: ieast = 4 !< 4 - index of eastward facing wall 731 INTEGER(iwp), PARAMETER :: iwest = 5 !< 5 - index of westward facing wall 732 733 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0, 0, 0, 1,-1/) !< surface normal direction x indices 734 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0, 1,-1, 0, 0/) !< surface normal direction y indices 735 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1, 0, 0, 0, 0/) !< surface normal direction z indices 734 736 REAL(wp), DIMENSION(0:nsurf_type) :: facearea !< area of single face in respective 735 737 !< direction (will be calc'd) 736 738 737 738 !-- indices and sizes of urban and land surface models739 INTEGER(iwp) :: startland !< start index of block of land and roof surfaces740 INTEGER(iwp) :: endland !< end index of block of land and roof surfaces741 INTEGER(iwp) :: nlands !< number of land and roof surfaces in local processor742 INTEGER(iwp) :: startwall !< start index of block of wall surfaces743 INTEGER(iwp) :: endwall !< end index of block of wall surfaces744 INTEGER(iwp) :: nwalls !< number of wall surfaces in local processor745 739 746 740 !-- indices needed for RTM netcdf output subroutines 747 741 INTEGER(iwp), PARAMETER :: nd = 5 748 742 CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 749 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /) 750 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /) 751 INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart 752 INTEGER(iwp), DIMENSION(0:nd-1) :: dirend 743 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup, isouth, inorth, iwest, ieast /) 753 744 754 745 !-- indices and sizes of urban and land surface models … … 762 753 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfstart !< starts of blocks of surfaces for individual processors in array surf (indexed from 1) 763 754 !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1] 755 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfl_col_start !< start of surfaces in surfl for each x,y column (local surfaces) 756 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfg_col_start !< start of surfaces in surfl for each x,y column (all surfaces) 764 757 765 758 !-- block variables needed for calculation of the plant canopy model inside the urban surface model … … 786 779 REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha) 787 780 INTEGER(iwp), PARAMETER :: rad_version_len = 10 !< length of identification string of rad version 788 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 3.0' !< identification of version of binary svf and restart files781 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 4.0' !< identification of version of binary svf and restart files 789 782 INTEGER(iwp) :: raytrace_discrete_elevs = 40 !< number of discretization steps for elevation (nadir to zenith) 790 783 INTEGER(iwp) :: raytrace_discrete_azims = 80 !< number of discretization steps for azimuth (out of 360 degrees) … … 879 872 #endif 880 873 REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth 881 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing 874 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterrt !< temporary global arrays for raytracing 875 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterrb !< temporary global arrays for raytracing 876 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: plantt !< temporary global arrays for raytracing 882 877 INTEGER(iwp) :: plantt_max 883 878 … … 1093 1088 cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, & 1094 1089 idir, jdir, kdir, id, iz, iy, ix, & 1095 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 1096 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 1097 nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf, & 1090 iup, inorth, isouth, ieast, iwest, & 1091 nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf, & 1098 1092 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 1099 radiation_interactions, startwall, startland, endland, endwall, & 1100 skyvf, skyvft, radiation_interactions_on, average_radiation, & 1101 rad_sw_in_diff, rad_sw_in_dir 1093 radiation_interactions, skyvf, skyvft, radiation_interactions_on, & 1094 average_radiation, rad_sw_in_diff, rad_sw_in_dir 1102 1095 1103 1096 … … 5966 5959 ! 5967 5960 !-- Set up thermal radiation from surfaces 5968 !-- Horizontal walls5969 5961 mm = 1 5970 5962 DO i = nxl, nxr 5971 5963 DO j = nys, nyn 5972 5964 ! 5965 !-- Horizontal walls 5973 5966 !-- urban 5974 5967 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) … … 5996 5989 mm = mm + 1 5997 5990 ENDDO 5998 ENDDO 5999 ENDDO 6000 ! 6001 !-- Vertical walls 6002 DO i = nxl, nxr 6003 DO j = nys, nyn 5991 ! 5992 !-- Downward facing default surfaces 5993 !-- TODO: Downward-facing surfaces are only available as default surfaces which lack 5994 !-- temperature and heat balance, therefore they have to be simulated as neutral to 5995 !-- radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit 5996 !-- radiation and they have zero net heat flux. 5997 DO m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i) 5998 surfoutll(mm) = 0.0_wp 5999 albedo_surf(mm) = 1.0_wp 6000 emiss_surf(mm) = 0.0_wp 6001 mm = mm + 1 6002 ENDDO 6003 ! 6004 !-- Vertical walls 6004 6005 DO l = 0, 3 6005 6006 ! … … 6400 6401 !-- Transfer radiation arrays required for energy balance to the respective data types 6401 6402 ! and claculate relevant radiation model-RTM coupling terms 6403 mm = 1 6404 DO i = nxl, nxr 6405 DO j = nys, nyn 6406 ! 6407 !-- Horizontal walls 6408 !-- urban 6409 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 6410 surf_usm_h%rad_sw_in(m) = surfinsw(mm) 6411 surf_usm_h%rad_sw_out(m) = surfoutsw(mm) 6412 surf_usm_h%rad_sw_dir(m) = surfinswdir(mm) 6413 surf_usm_h%rad_sw_dif(m) = surfinswdif(mm) 6414 surf_usm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) - & 6415 surfinswdif(mm) 6416 surf_usm_h%rad_sw_res(m) = surfins(mm) 6417 surf_usm_h%rad_lw_in(m) = surfinlw(mm) 6418 surf_usm_h%rad_lw_out(m) = surfoutlw(mm) 6419 surf_usm_h%rad_net(m) = surfinsw(mm) - surfoutsw(mm) + & 6420 surfinlw(mm) - surfoutlw(mm) 6421 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net(m) 6422 surf_usm_h%rad_lw_dif(m) = surfinlwdif(mm) 6423 surf_usm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm) 6424 surf_usm_h%rad_lw_res(m) = surfinl(mm) 6425 mm = mm + 1 6426 ENDDO 6427 ! 6428 !-- land 6429 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 6430 surf_lsm_h%rad_sw_in(m) = surfinsw(mm) 6431 surf_lsm_h%rad_sw_out(m) = surfoutsw(mm) 6432 surf_lsm_h%rad_sw_dir(m) = surfinswdir(mm) 6433 surf_lsm_h%rad_sw_dif(m) = surfinswdif(mm) 6434 surf_lsm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) - & 6435 surfinswdif(mm) 6436 surf_lsm_h%rad_sw_res(m) = surfins(mm) 6437 surf_lsm_h%rad_lw_in(m) = surfinlw(mm) 6438 surf_lsm_h%rad_lw_out(m) = surfoutlw(mm) 6439 surf_lsm_h%rad_net(m) = surfinsw(mm) - surfoutsw(mm) + & 6440 surfinlw(mm) - surfoutlw(mm) 6441 surf_lsm_h%rad_lw_dif(m) = surfinlwdif(mm) 6442 surf_lsm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm) 6443 surf_lsm_h%rad_lw_res(m) = surfinl(mm) 6444 mm = mm + 1 6445 ENDDO 6446 ! 6447 !-- Downward facing default surfaces 6448 !-- TODO: Downward-facing surfaces are only available as default surfaces which lack 6449 !-- temperature and heat balance, therefore they have to be simulated as neutral to 6450 !-- radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit 6451 !-- radiation and they have zero net heat flux. 6452 DO m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i) 6453 mm = mm + 1 6454 ENDDO 6455 ! 6456 !-- Vertical walls 6457 DO l = 0, 3 6458 ! 6459 !-- urban 6460 DO m = surf_usm_v(l)%start_index(j,i), & 6461 surf_usm_v(l)%end_index(j,i) 6462 surf_usm_v(l)%rad_sw_in(m) = surfinsw(mm) 6463 surf_usm_v(l)%rad_sw_out(m) = surfoutsw(mm) 6464 surf_usm_v(l)%rad_sw_dir(m) = surfinswdir(mm) 6465 surf_usm_v(l)%rad_sw_dif(m) = surfinswdif(mm) 6466 surf_usm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) - & 6467 surfinswdif(mm) 6468 surf_usm_v(l)%rad_sw_res(m) = surfins(mm) 6469 surf_usm_v(l)%rad_lw_in(m) = surfinlw(mm) 6470 surf_usm_v(l)%rad_lw_out(m) = surfoutlw(mm) 6471 surf_usm_v(l)%rad_net(m) = surfinsw(mm) - surfoutsw(mm) + & 6472 surfinlw(mm) - surfoutlw(mm) 6473 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_net(m) 6474 surf_usm_v(l)%rad_lw_dif(m) = surfinlwdif(mm) 6475 surf_usm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm) 6476 surf_usm_v(l)%rad_lw_res(m) = surfinl(mm) 6477 mm = mm + 1 6478 ENDDO 6479 ! 6480 !-- land 6481 DO m = surf_lsm_v(l)%start_index(j,i), & 6482 surf_lsm_v(l)%end_index(j,i) 6483 surf_lsm_v(l)%rad_sw_in(m) = surfinsw(mm) 6484 surf_lsm_v(l)%rad_sw_out(m) = surfoutsw(mm) 6485 surf_lsm_v(l)%rad_sw_dir(m) = surfinswdir(mm) 6486 surf_lsm_v(l)%rad_sw_dif(m) = surfinswdif(mm) 6487 surf_lsm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) - & 6488 surfinswdif(mm) 6489 surf_lsm_v(l)%rad_sw_res(m) = surfins(mm) 6490 surf_lsm_v(l)%rad_lw_in(m) = surfinlw(mm) 6491 surf_lsm_v(l)%rad_lw_out(m) = surfoutlw(mm) 6492 surf_lsm_v(l)%rad_net(m) = surfinsw(mm) - surfoutsw(mm) + & 6493 surfinlw(mm) - surfoutlw(mm) 6494 surf_lsm_v(l)%rad_lw_dif(m) = surfinlwdif(mm) 6495 surf_lsm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm) 6496 surf_lsm_v(l)%rad_lw_res(m) = surfinl(mm) 6497 mm = mm + 1 6498 ENDDO 6499 ENDDO 6500 ENDDO 6501 ENDDO 6402 6502 6403 6503 DO i = 1, nsurfl 6404 m = surfl(im,i)6405 6504 d = surfl(id, i) 6406 !6407 !-- (1) Urban surfaces6408 !-- upward-facing6409 IF ( surfl(1,i) == iup_u ) THEN6410 surf_usm_h%rad_sw_in(m) = surfinsw(i)6411 surf_usm_h%rad_sw_out(m) = surfoutsw(i)6412 surf_usm_h%rad_sw_dir(m) = surfinswdir(i)6413 surf_usm_h%rad_sw_dif(m) = surfinswdif(i)6414 surf_usm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6415 surfinswdif(i)6416 surf_usm_h%rad_sw_res(m) = surfins(i)6417 surf_usm_h%rad_lw_in(m) = surfinlw(i)6418 surf_usm_h%rad_lw_out(m) = surfoutlw(i)6419 surf_usm_h%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6420 surfinlw(i) - surfoutlw(i)6421 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net(m)6422 surf_usm_h%rad_lw_dif(m) = surfinlwdif(i)6423 surf_usm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6424 surf_usm_h%rad_lw_res(m) = surfinl(i)6425 !6426 !-- northward-facding6427 ELSEIF ( surfl(1,i) == inorth_u ) THEN6428 surf_usm_v(0)%rad_sw_in(m) = surfinsw(i)6429 surf_usm_v(0)%rad_sw_out(m) = surfoutsw(i)6430 surf_usm_v(0)%rad_sw_dir(m) = surfinswdir(i)6431 surf_usm_v(0)%rad_sw_dif(m) = surfinswdif(i)6432 surf_usm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6433 surfinswdif(i)6434 surf_usm_v(0)%rad_sw_res(m) = surfins(i)6435 surf_usm_v(0)%rad_lw_in(m) = surfinlw(i)6436 surf_usm_v(0)%rad_lw_out(m) = surfoutlw(i)6437 surf_usm_v(0)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6438 surfinlw(i) - surfoutlw(i)6439 surf_usm_v(0)%rad_net_l(m) = surf_usm_v(0)%rad_net(m)6440 surf_usm_v(0)%rad_lw_dif(m) = surfinlwdif(i)6441 surf_usm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6442 surf_usm_v(0)%rad_lw_res(m) = surfinl(i)6443 !6444 !-- southward-facding6445 ELSEIF ( surfl(1,i) == isouth_u ) THEN6446 surf_usm_v(1)%rad_sw_in(m) = surfinsw(i)6447 surf_usm_v(1)%rad_sw_out(m) = surfoutsw(i)6448 surf_usm_v(1)%rad_sw_dir(m) = surfinswdir(i)6449 surf_usm_v(1)%rad_sw_dif(m) = surfinswdif(i)6450 surf_usm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6451 surfinswdif(i)6452 surf_usm_v(1)%rad_sw_res(m) = surfins(i)6453 surf_usm_v(1)%rad_lw_in(m) = surfinlw(i)6454 surf_usm_v(1)%rad_lw_out(m) = surfoutlw(i)6455 surf_usm_v(1)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6456 surfinlw(i) - surfoutlw(i)6457 surf_usm_v(1)%rad_net_l(m) = surf_usm_v(1)%rad_net(m)6458 surf_usm_v(1)%rad_lw_dif(m) = surfinlwdif(i)6459 surf_usm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6460 surf_usm_v(1)%rad_lw_res(m) = surfinl(i)6461 !6462 !-- eastward-facing6463 ELSEIF ( surfl(1,i) == ieast_u ) THEN6464 surf_usm_v(2)%rad_sw_in(m) = surfinsw(i)6465 surf_usm_v(2)%rad_sw_out(m) = surfoutsw(i)6466 surf_usm_v(2)%rad_sw_dir(m) = surfinswdir(i)6467 surf_usm_v(2)%rad_sw_dif(m) = surfinswdif(i)6468 surf_usm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6469 surfinswdif(i)6470 surf_usm_v(2)%rad_sw_res(m) = surfins(i)6471 surf_usm_v(2)%rad_lw_in(m) = surfinlw(i)6472 surf_usm_v(2)%rad_lw_out(m) = surfoutlw(i)6473 surf_usm_v(2)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6474 surfinlw(i) - surfoutlw(i)6475 surf_usm_v(2)%rad_net_l(m) = surf_usm_v(2)%rad_net(m)6476 surf_usm_v(2)%rad_lw_dif(m) = surfinlwdif(i)6477 surf_usm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6478 surf_usm_v(2)%rad_lw_res(m) = surfinl(i)6479 !6480 !-- westward-facding6481 ELSEIF ( surfl(1,i) == iwest_u ) THEN6482 surf_usm_v(3)%rad_sw_in(m) = surfinsw(i)6483 surf_usm_v(3)%rad_sw_out(m) = surfoutsw(i)6484 surf_usm_v(3)%rad_sw_dir(m) = surfinswdir(i)6485 surf_usm_v(3)%rad_sw_dif(m) = surfinswdif(i)6486 surf_usm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6487 surfinswdif(i)6488 surf_usm_v(3)%rad_sw_res(m) = surfins(i)6489 surf_usm_v(3)%rad_lw_in(m) = surfinlw(i)6490 surf_usm_v(3)%rad_lw_out(m) = surfoutlw(i)6491 surf_usm_v(3)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6492 surfinlw(i) - surfoutlw(i)6493 surf_usm_v(3)%rad_net_l(m) = surf_usm_v(3)%rad_net(m)6494 surf_usm_v(3)%rad_lw_dif(m) = surfinlwdif(i)6495 surf_usm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6496 surf_usm_v(3)%rad_lw_res(m) = surfinl(i)6497 !6498 !-- (2) land surfaces6499 !-- upward-facing6500 ELSEIF ( surfl(1,i) == iup_l ) THEN6501 surf_lsm_h%rad_sw_in(m) = surfinsw(i)6502 surf_lsm_h%rad_sw_out(m) = surfoutsw(i)6503 surf_lsm_h%rad_sw_dir(m) = surfinswdir(i)6504 surf_lsm_h%rad_sw_dif(m) = surfinswdif(i)6505 surf_lsm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6506 surfinswdif(i)6507 surf_lsm_h%rad_sw_res(m) = surfins(i)6508 surf_lsm_h%rad_lw_in(m) = surfinlw(i)6509 surf_lsm_h%rad_lw_out(m) = surfoutlw(i)6510 surf_lsm_h%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6511 surfinlw(i) - surfoutlw(i)6512 surf_lsm_h%rad_lw_dif(m) = surfinlwdif(i)6513 surf_lsm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6514 surf_lsm_h%rad_lw_res(m) = surfinl(i)6515 !6516 !-- northward-facding6517 ELSEIF ( surfl(1,i) == inorth_l ) THEN6518 surf_lsm_v(0)%rad_sw_in(m) = surfinsw(i)6519 surf_lsm_v(0)%rad_sw_out(m) = surfoutsw(i)6520 surf_lsm_v(0)%rad_sw_dir(m) = surfinswdir(i)6521 surf_lsm_v(0)%rad_sw_dif(m) = surfinswdif(i)6522 surf_lsm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6523 surfinswdif(i)6524 surf_lsm_v(0)%rad_sw_res(m) = surfins(i)6525 surf_lsm_v(0)%rad_lw_in(m) = surfinlw(i)6526 surf_lsm_v(0)%rad_lw_out(m) = surfoutlw(i)6527 surf_lsm_v(0)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6528 surfinlw(i) - surfoutlw(i)6529 surf_lsm_v(0)%rad_lw_dif(m) = surfinlwdif(i)6530 surf_lsm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6531 surf_lsm_v(0)%rad_lw_res(m) = surfinl(i)6532 !6533 !-- southward-facding6534 ELSEIF ( surfl(1,i) == isouth_l ) THEN6535 surf_lsm_v(1)%rad_sw_in(m) = surfinsw(i)6536 surf_lsm_v(1)%rad_sw_out(m) = surfoutsw(i)6537 surf_lsm_v(1)%rad_sw_dir(m) = surfinswdir(i)6538 surf_lsm_v(1)%rad_sw_dif(m) = surfinswdif(i)6539 surf_lsm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6540 surfinswdif(i)6541 surf_lsm_v(1)%rad_sw_res(m) = surfins(i)6542 surf_lsm_v(1)%rad_lw_in(m) = surfinlw(i)6543 surf_lsm_v(1)%rad_lw_out(m) = surfoutlw(i)6544 surf_lsm_v(1)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6545 surfinlw(i) - surfoutlw(i)6546 surf_lsm_v(1)%rad_lw_dif(m) = surfinlwdif(i)6547 surf_lsm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6548 surf_lsm_v(1)%rad_lw_res(m) = surfinl(i)6549 !6550 !-- eastward-facing6551 ELSEIF ( surfl(1,i) == ieast_l ) THEN6552 surf_lsm_v(2)%rad_sw_in(m) = surfinsw(i)6553 surf_lsm_v(2)%rad_sw_out(m) = surfoutsw(i)6554 surf_lsm_v(2)%rad_sw_dir(m) = surfinswdir(i)6555 surf_lsm_v(2)%rad_sw_dif(m) = surfinswdif(i)6556 surf_lsm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6557 surfinswdif(i)6558 surf_lsm_v(2)%rad_sw_res(m) = surfins(i)6559 surf_lsm_v(2)%rad_lw_in(m) = surfinlw(i)6560 surf_lsm_v(2)%rad_lw_out(m) = surfoutlw(i)6561 surf_lsm_v(2)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6562 surfinlw(i) - surfoutlw(i)6563 surf_lsm_v(2)%rad_lw_dif(m) = surfinlwdif(i)6564 surf_lsm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6565 surf_lsm_v(2)%rad_lw_res(m) = surfinl(i)6566 !6567 !-- westward-facing6568 ELSEIF ( surfl(1,i) == iwest_l ) THEN6569 surf_lsm_v(3)%rad_sw_in(m) = surfinsw(i)6570 surf_lsm_v(3)%rad_sw_out(m) = surfoutsw(i)6571 surf_lsm_v(3)%rad_sw_dir(m) = surfinswdir(i)6572 surf_lsm_v(3)%rad_sw_dif(m) = surfinswdif(i)6573 surf_lsm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) - &6574 surfinswdif(i)6575 surf_lsm_v(3)%rad_sw_res(m) = surfins(i)6576 surf_lsm_v(3)%rad_lw_in(m) = surfinlw(i)6577 surf_lsm_v(3)%rad_lw_out(m) = surfoutlw(i)6578 surf_lsm_v(3)%rad_net(m) = surfinsw(i) - surfoutsw(i) + &6579 surfinlw(i) - surfoutlw(i)6580 surf_lsm_v(3)%rad_lw_dif(m) = surfinlwdif(i)6581 surf_lsm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)6582 surf_lsm_v(3)%rad_lw_res(m) = surfinl(i)6583 ENDIF6584 6505 ! 6585 6506 !-- RTM coupling terms … … 6942 6863 INTEGER(iwp) :: k_topo !< vertical index indicating 6943 6864 !< topography top for given (j,i) 6865 INTEGER(iwp) :: icol !< flat column number (in (y,x) fortran order) 6944 6866 INTEGER(iwp) :: isurf, ipcgb, imrt 6945 6867 INTEGER(iwp) :: nzptl, nzubl, nzutl … … 6965 6887 !-- for any upward-facing wall (see bit 12). 6966 6888 nzubl = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ) 6967 nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr, 0) )6889 nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,5) ) 6968 6890 6969 6891 nzubl = MAX( nzubl, nzb ) … … 7077 6999 !-- Number of horizontal surfaces including land- and roof surfaces in both USM and LSM. Note that 7078 7000 !-- All horizontal surface elements are already counted in surface_mod. 7079 startland = 1 7080 nsurfl = surf_usm_h%ns + surf_lsm_h%ns 7081 endland = nsurfl 7082 nlands = endland - startland + 1 7083 7001 nsurfl = surf_usm_h%ns + surf_lsm_h%ns + surf_def_h(1)%ns 7084 7002 ! 7085 7003 !-- Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are 7086 7004 !-- already counted in surface_mod. 7087 startwall = nsurfl+17088 7005 DO i = 0,3 7089 7006 nsurfl = nsurfl + surf_usm_v(i)%ns + surf_lsm_v(i)%ns 7090 7007 ENDDO 7091 endwall = nsurfl7092 nwalls = endwall - startwall + 17093 dirstart = (/ startland, startwall, startwall, startwall, startwall /)7094 dirend = (/ endland, endwall, endwall, endwall, endwall /)7095 7008 7096 7009 !-- fill gridpcbl and pcbl … … 7196 7109 ENDIF 7197 7110 7198 !-- add horizontal surface elements (land and urban surfaces) 7199 !-- TODO: add urban overhanging surfaces (idown_u) 7111 ALLOCATE( surfl_col_start(0:nnx*nny-1), surfg_col_start(0:(nx+1)*(ny+1)) ) 7112 ! 7113 !-- Add horizontal and vertical surface elements (land and urban surfaces) 7114 !-- ordered by x,y column (y most varying) 7115 !-- TODO: remove the hard coding of l = 0 to l = idirection 7116 icol = 0 7200 7117 DO i = nxl, nxr 7201 7118 DO j = nys, nyn 7119 ! 7120 !-- Save column start 7121 surfl_col_start(icol) = isurf + 1 7122 icol = icol + 1 7123 ! 7124 !-- Horizontal surfaces 7202 7125 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 7203 7126 k = surf_usm_h%k(m) 7204 7127 isurf = isurf + 1 7205 surfl(:,isurf) = (/iup_u,k,j,i,m/) 7206 IF ( rad_angular_discretization ) THEN 7207 gridsurf(iup_u,k,j,i) = isurf 7208 ENDIF 7128 surfl(:,isurf) = (/iup,k,j,i/) 7209 7129 ENDDO 7210 7130 … … 7212 7132 k = surf_lsm_h%k(m) 7213 7133 isurf = isurf + 1 7214 surfl(:,isurf) = (/iup_l,k,j,i,m/) 7215 IF ( rad_angular_discretization ) THEN 7216 gridsurf(iup_u,k,j,i) = isurf 7217 ENDIF 7134 surfl(:,isurf) = (/iup,k,j,i/) 7218 7135 ENDDO 7219 7136 7220 ENDDO7221 ENDDO7222 7223 !-- add vertical surface elements (land and urban surfaces)7224 !-- TODO: remove the hard coding of l = 0 to l = idirection 7225 DO i = nxl, nxr 7226 DO j = nys, nyn 7137 DO m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i) 7138 k = surf_def_h(1)%k(m) 7139 isurf = isurf + 1 7140 surfl(:,isurf) = (/idown,k,j,i/) 7141 ENDDO 7142 ! 7143 !-- Vertical surfaces 7227 7144 l = 0 7228 7145 DO m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i) 7229 7146 k = surf_usm_v(l)%k(m) 7230 7147 isurf = isurf + 1 7231 surfl(:,isurf) = (/inorth_u,k,j,i,m/) 7232 IF ( rad_angular_discretization ) THEN 7233 gridsurf(inorth_u,k,j,i) = isurf 7234 ENDIF 7148 surfl(:,isurf) = (/inorth,k,j,i/) 7235 7149 ENDDO 7236 7150 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 7237 7151 k = surf_lsm_v(l)%k(m) 7238 7152 isurf = isurf + 1 7239 surfl(:,isurf) = (/inorth_l,k,j,i,m/) 7240 IF ( rad_angular_discretization ) THEN 7241 gridsurf(inorth_u,k,j,i) = isurf 7242 ENDIF 7153 surfl(:,isurf) = (/inorth,k,j,i/) 7243 7154 ENDDO 7244 7155 … … 7247 7158 k = surf_usm_v(l)%k(m) 7248 7159 isurf = isurf + 1 7249 surfl(:,isurf) = (/isouth_u,k,j,i,m/) 7250 IF ( rad_angular_discretization ) THEN 7251 gridsurf(isouth_u,k,j,i) = isurf 7252 ENDIF 7160 surfl(:,isurf) = (/isouth,k,j,i/) 7253 7161 ENDDO 7254 7162 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 7255 7163 k = surf_lsm_v(l)%k(m) 7256 7164 isurf = isurf + 1 7257 surfl(:,isurf) = (/isouth_l,k,j,i,m/) 7258 IF ( rad_angular_discretization ) THEN 7259 gridsurf(isouth_u,k,j,i) = isurf 7260 ENDIF 7165 surfl(:,isurf) = (/isouth,k,j,i/) 7261 7166 ENDDO 7262 7167 … … 7265 7170 k = surf_usm_v(l)%k(m) 7266 7171 isurf = isurf + 1 7267 surfl(:,isurf) = (/ieast_u,k,j,i,m/) 7268 IF ( rad_angular_discretization ) THEN 7269 gridsurf(ieast_u,k,j,i) = isurf 7270 ENDIF 7172 surfl(:,isurf) = (/ieast,k,j,i/) 7271 7173 ENDDO 7272 7174 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 7273 7175 k = surf_lsm_v(l)%k(m) 7274 7176 isurf = isurf + 1 7275 surfl(:,isurf) = (/ieast_l,k,j,i,m/) 7276 IF ( rad_angular_discretization ) THEN 7277 gridsurf(ieast_u,k,j,i) = isurf 7278 ENDIF 7177 surfl(:,isurf) = (/ieast,k,j,i/) 7279 7178 ENDDO 7280 7179 … … 7283 7182 k = surf_usm_v(l)%k(m) 7284 7183 isurf = isurf + 1 7285 surfl(:,isurf) = (/iwest_u,k,j,i,m/) 7286 IF ( rad_angular_discretization ) THEN 7287 gridsurf(iwest_u,k,j,i) = isurf 7288 ENDIF 7184 surfl(:,isurf) = (/iwest,k,j,i/) 7289 7185 ENDDO 7290 7186 DO m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i) 7291 7187 k = surf_lsm_v(l)%k(m) 7292 7188 isurf = isurf + 1 7293 surfl(:,isurf) = (/iwest_l,k,j,i,m/) 7294 IF ( rad_angular_discretization ) THEN 7295 gridsurf(iwest_u,k,j,i) = isurf 7296 ENDIF 7189 surfl(:,isurf) = (/iwest,k,j,i/) 7297 7190 ENDDO 7298 7191 ENDDO … … 7406 7299 surf = surfl 7407 7300 #endif 7408 7409 !-- 7301 ! 7302 !-- Populate gridsurf with reverse global indices (->surf) 7303 IF ( rad_angular_discretization ) THEN 7304 DO isurf = 1, nsurfl 7305 gridsurf(surfl(id,isurf),surfl(iz,isurf), & 7306 surfl(iy,isurf),surfl(ix,isurf)) = isurf + surfstart(myid) 7307 ENDDO 7308 ENDIF 7309 ! 7310 !-- Gather column start indices 7311 #if defined( __parallel ) 7312 CALL MPI_Allgather(surfl_col_start, nnx*nny, MPI_INTEGER, & 7313 surfg_col_start, nnx*nny, MPI_INTEGER, comm2d, ierr) 7314 IF ( ierr /= 0 ) THEN 7315 WRITE(9,*) 'Error MPI_AllGather1b:', ierr, surfl_col_start, surfg_col_start 7316 FLUSH(9) 7317 ENDIF 7318 ! 7319 !-- Convert local indices (->surfl) to global (->surf) 7320 DO i = 0, numprocs-1 7321 surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) = & 7322 surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) + surfstart(i) 7323 ENDDO 7324 #else 7325 surfg_col_start(0:(nx+1)*(ny+1)-1) = surfl_col_start(0:(nx+1)*(ny+1)-1) 7326 #endif 7327 surfg_col_start((nx+1)*(ny+1)) = nsurf+1 7328 ! 7410 7329 !-- allocation of the arrays for direct and diffusion radiation 7411 7330 IF ( debug_output ) CALL debug_message( 'allocation of radiation arrays', 'info' ) … … 7483 7402 REAL(wp), DIMENSION(:), ALLOCATABLE :: vffrac0 !< dtto (original values) 7484 7403 REAL(wp), DIMENSION(:), ALLOCATABLE :: ztransp !< array of transparency in z steps 7485 INTEGER(iwp) :: lowest_free_ray !< index into zdirs7486 7404 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: itarget !< face indices of detected obstacles 7487 7405 INTEGER(iwp) :: itarg0, itarg1 … … 7502 7420 INTEGER(iwp) :: max_track_len !< maximum 2d track length 7503 7421 #if defined( __parallel ) 7504 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l 7505 INTEGER(iwp), DIMENSION(:,:), POINTER :: nzterrl 7422 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrtl_l 7423 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrbl_l 7424 INTEGER(iwp), DIMENSION(:,:), POINTER :: nzterrtl 7425 INTEGER(iwp), DIMENSION(:,:), POINTER :: nzterrbl 7506 7426 INTEGER(iwp) :: minfo 7507 7427 REAL(wp), DIMENSION(:), POINTER, SAVE :: lad_s_rma !< fortran 1D pointer … … 7540 7460 7541 7461 !-- initialize temporary terrain and plant canopy height arrays (global 2D array!) 7542 ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) ) 7462 ALLOCATE( nzterrt(0:(nx+1)*(ny+1)-1) ) 7463 ALLOCATE( nzterrb(0:(nx+1)*(ny+1)-1) ) 7543 7464 #if defined( __parallel ) 7544 !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) ) 7545 ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) ) 7546 nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1)) 7547 nzterrl = topo_top_ind(nys:nyn,nxl:nxr,0) 7548 CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, & 7549 nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr ) 7465 ALLOCATE( nzterrtl_l((nyn-nys+1)*(nxr-nxl+1)) ) 7466 ALLOCATE( nzterrbl_l((nyn-nys+1)*(nxr-nxl+1)) ) 7467 nzterrtl(nys:nyn,nxl:nxr) => nzterrtl_l(1:(nyn-nys+1)*(nxr-nxl+1)) 7468 nzterrbl(nys:nyn,nxl:nxr) => nzterrbl_l(1:(nyn-nys+1)*(nxr-nxl+1)) 7469 nzterrtl = topo_top_ind(nys:nyn,nxl:nxr,5) 7470 nzterrbl = topo_top_ind(nys:nyn,nxl:nxr,0) 7471 CALL MPI_AllGather( nzterrtl_l, nnx*nny, MPI_INTEGER, & 7472 nzterrt, nnx*nny, MPI_INTEGER, comm2d, ierr ) 7550 7473 IF ( ierr /= 0 ) THEN 7551 WRITE(9,*) 'Error MPI_AllGather1 :', ierr, SIZE(nzterrl_l), nnx*nny, &7552 SIZE(nzterr ), nnx*nny7474 WRITE(9,*) 'Error MPI_AllGather1t:', ierr, SIZE(nzterrtl_l), nnx*nny, & 7475 SIZE(nzterrt), nnx*nny 7553 7476 FLUSH(9) 7554 7477 ENDIF 7555 DEALLOCATE(nzterrl_l) 7478 CALL MPI_AllGather( nzterrbl_l, nnx*nny, MPI_INTEGER, & 7479 nzterrb, nnx*nny, MPI_INTEGER, comm2d, ierr ) 7480 IF ( ierr /= 0 ) THEN 7481 WRITE(9,*) 'Error MPI_AllGather1b:', ierr, SIZE(nzterrbl_l), nnx*nny, & 7482 SIZE(nzterrb), nnx*nny 7483 FLUSH(9) 7484 ENDIF 7485 DEALLOCATE(nzterrtl_l) 7486 DEALLOCATE(nzterrbl_l) 7556 7487 #else 7557 nzterr = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) ) 7488 nzterrt = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,5), (/(nx+1)*(ny+1)/) ) 7489 nzterrb = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) ) 7558 7490 #endif 7559 7491 IF ( plant_canopy ) THEN … … 7701 7633 !--Select a proper half-sphere for 2D raytracing 7702 7634 SELECT CASE ( td ) 7703 CASE ( iup _u, iup_l)7635 CASE ( iup ) 7704 7636 az0 = 0._wp 7705 7637 naz = raytrace_discrete_azims … … 7708 7640 nzn = raytrace_discrete_elevs / 2 7709 7641 zns = pi / 2._wp / REAL(nzn, wp) 7710 CASE ( isouth_u, isouth_l ) 7642 CASE ( idown ) 7643 az0 = 0._wp 7644 naz = raytrace_discrete_azims 7645 azs = 2._wp * pi / REAL(naz, wp) 7646 zn0 = pi / 2._wp 7647 nzn = raytrace_discrete_elevs / 2 7648 zns = pi / 2._wp / REAL(nzn, wp) 7649 CASE ( isouth ) 7711 7650 az0 = pi / 2._wp 7712 7651 naz = raytrace_discrete_azims / 2 … … 7715 7654 nzn = raytrace_discrete_elevs 7716 7655 zns = pi / REAL(nzn, wp) 7717 CASE ( inorth _u, inorth_l)7656 CASE ( inorth ) 7718 7657 az0 = - pi / 2._wp 7719 7658 naz = raytrace_discrete_azims / 2 … … 7722 7661 nzn = raytrace_discrete_elevs 7723 7662 zns = pi / REAL(nzn, wp) 7724 CASE ( iwest _u, iwest_l)7663 CASE ( iwest ) 7725 7664 az0 = pi 7726 7665 naz = raytrace_discrete_azims / 2 … … 7729 7668 nzn = raytrace_discrete_elevs 7730 7669 zns = pi / REAL(nzn, wp) 7731 CASE ( ieast _u, ieast_l)7670 CASE ( ieast ) 7732 7671 az0 = 0._wp 7733 7672 naz = raytrace_discrete_azims / 2 … … 7744 7683 7745 7684 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), & 7746 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 7747 !in case of rad_angular_discretization 7685 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) 7748 7686 7749 7687 itarg0 = 1 … … 7751 7689 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 7752 7690 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 7753 IF ( td == iup _u .OR. td == iup_l) THEN7691 IF ( td == iup ) THEN 7754 7692 vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp) 7755 7693 ! … … 7759 7697 ENDDO 7760 7698 !-- sum of whole vffrac equals 1, verified 7699 ELSEIF ( td == idown ) THEN 7700 vffrac(1:nzn) = -(COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp) 7701 DO iaz = 1, naz-1 7702 vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn) 7703 ENDDO 7704 !-- sum of whole vffrac equals 1, verified 7761 7705 ENDIF 7762 7706 ! … … 7764 7708 DO iaz = 1, naz 7765 7709 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 7766 IF ( td /= iup _u .AND. td /= iup_l) THEN7710 IF ( td /= iup .AND. td /= idown ) THEN 7767 7711 az2 = REAL(iaz, wp) * azs - pi/2._wp 7768 7712 az1 = az2 - azs … … 7783 7727 surfstart(myid) + isurflt, facearea(td), & 7784 7728 vffrac(itarg0:itarg1), .TRUE., .TRUE., & 7785 .FALSE., lowest_free_ray, & 7786 ztransp(itarg0:itarg1), & 7729 .FALSE., ztransp(itarg0:itarg1), & 7787 7730 itarget(itarg0:itarg1)) 7788 7731 7789 skyvf(isurflt) = skyvf(isurflt) + & 7790 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) 7791 skyvft(isurflt) = skyvft(isurflt) + & 7792 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) & 7793 * vffrac(itarg0:itarg0+lowest_free_ray-1)) 7794 7732 skyvf(isurflt) = skyvf(isurflt) + & 7733 SUM(vffrac(itarg0:itarg1), & 7734 MASK=(itarget(itarg0:itarg1)<0)) 7735 skyvft(isurflt) = skyvft(isurflt) + & 7736 SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), & 7737 MASK=(itarget(itarg0:itarg1)<0)) 7738 ! 7795 7739 !-- Save direct solar transparency 7796 7740 j = MODULO(NINT(azmid/ & 7797 7741 (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), & 7798 7742 raytrace_discrete_azims) 7799 7800 DO k = 1, raytrace_discrete_elevs/2 7801 i = dsidir_rev(k-1, j) 7802 IF ( i /= -1 .AND. k <= lowest_free_ray ) & 7803 dsitrans(isurflt, i) = ztransp(itarg0+k-1) 7804 ENDDO 7805 7743 ! 7744 !-- For down direction there is no direct irradiance, otherwise az0=0 7745 IF ( td /= idown ) THEN 7746 DO k = 1, raytrace_discrete_elevs/2 7747 i = dsidir_rev(k-1, j) 7748 IF ( i /= -1 .AND. itarget(itarg0+k-1) < 0 ) & 7749 dsitrans(isurflt, i) = ztransp(itarg0+k-1) 7750 ENDDO 7751 ENDIF 7806 7752 ! 7807 7753 !-- Advance itarget indices … … 7820 7766 itarg0 = 1 7821 7767 DO WHILE ( itarg0 <= nzn*naz ) 7822 IF ( itarget(itarg0) /= -1) EXIT7768 IF ( itarget(itarg0) >= 0 ) EXIT 7823 7769 itarg0 = itarg0 + 1 7824 7770 ENDDO … … 7983 7929 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 7984 7930 -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., & 7985 lowest_free_ray,ztransp, itarget)7931 ztransp, itarget) 7986 7932 7987 7933 !-- Save direct solar transparency … … 7991 7937 DO k = 1, raytrace_discrete_elevs/2 7992 7938 i = dsidir_rev(k-1, j) 7993 IF ( i /= -1 .AND. k <= lowest_free_ray) &7939 IF ( i /= -1 .AND. itarget(k) < 0 ) & 7994 7940 dsitransc(ipcgb, i) = ztransp(k) 7995 7941 ENDDO … … 8010 7956 zns = pi / REAL(nzn, wp) 8011 7957 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), & 8012 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 8013 !in case of rad_angular_discretization 7958 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) 8014 7959 8015 7960 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) … … 8059 8004 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 8060 8005 -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., & 8061 .FALSE., .TRUE., lowest_free_ray, & 8062 ztransp(itarg0:itarg1), & 8006 .FALSE., .TRUE., ztransp(itarg0:itarg1), & 8063 8007 itarget(itarg0:itarg1)) 8064 8008 8065 8009 !-- Sky view factors for MRT 8066 8010 mrtsky(imrt) = mrtsky(imrt) + & 8067 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) 8011 SUM(vffrac(itarg0:itarg1), & 8012 MASK=(itarget(itarg0:itarg1)<0)) 8068 8013 mrtskyt(imrt) = mrtskyt(imrt) + & 8069 SUM(ztransp(itarg0:itarg 0+lowest_free_ray-1)&8070 * vffrac(itarg0:itarg0+lowest_free_ray-1))8014 SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), & 8015 MASK=(itarget(itarg0:itarg1)<0)) 8071 8016 !-- Direct solar transparency for MRT 8072 8017 j = MODULO(NINT(azmid/ & … … 8075 8020 DO k = 1, raytrace_discrete_elevs/2 8076 8021 i = dsidir_rev(k-1, j) 8077 IF ( i /= -1 .AND. k <= lowest_free_ray) &8022 IF ( i /= -1 .AND. itarget(itarg0+k-1) < 0 ) & 8078 8023 mrtdsit(imrt, i) = ztransp(itarg0+k-1) 8079 8024 ENDDO … … 8092 8037 itarg0 = 1 8093 8038 DO WHILE ( itarg0 <= nzn*naz ) 8094 IF ( itarget(itarg0) /= -1) EXIT8039 IF ( itarget(itarg0) >= 0 ) EXIT 8095 8040 itarg0 = itarg0 + 1 8096 8041 ENDDO … … 8183 8128 8184 8129 !-- deallocate temporary global arrays 8185 DEALLOCATE( nzterr)8130 DEALLOCATE( nzterrt, nzterrb ) 8186 8131 8187 8132 IF ( plant_canopy ) THEN … … 8607 8552 !-- in the array nzterr and plantt and id of the coresponding processor 8608 8553 CALL radiation_calc_global_offset( box(3), box(2), 0, 1, offs_glob=ig ) 8609 IF ( box(1) <= nzterr (ig) ) THEN8554 IF ( box(1) <= nzterrb(ig) ) THEN 8610 8555 visible = .FALSE. 8611 8556 RETURN … … 8703 8648 SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac, & 8704 8649 calc_svf, create_csf, skip_1st_pcb, & 8705 lowest_free_ray,transparency, itarget)8650 transparency, itarget) 8706 8651 IMPLICIT NONE 8707 8652 … … 8716 8661 LOGICAL, INTENT(in) :: create_csf !< whether to create canopy sink factors 8717 8662 LOGICAL, INTENT(in) :: skip_1st_pcb !< whether to skip first plant canopy box during raytracing 8718 INTEGER(iwp), INTENT(out) :: lowest_free_ray !< index into zdirs8719 8663 REAL(wp), DIMENSION(nrays), INTENT(OUT) :: transparency !< transparencies of zdirs paths 8720 INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) :: itarget !< global indices of target faces for zdirs 8721 8722 INTEGER(iwp), DIMENSION(nrays) :: target_procs 8723 REAL(wp) :: horizon !< highest horizon found after raytracing (z/hdist) 8664 INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) :: itarget !< global indices of target faces for zdirs or <0 for sky 8665 8666 LOGICAL :: is_mixed_col !< whether current column contains full-3D geometry 8667 INTEGER(iwp) :: lowest_free_ray !< index into zdirs 8668 INTEGER(iwp) :: lowest_mixed_ray !< index into zdirs 8669 REAL(wp) :: full_horizon !< highest full horizon found after raytracing (z/hdist) 8724 8670 INTEGER(iwp) :: i, k, l, d 8671 INTEGER(iwp) :: iray !< index into zdirs 8672 INTEGER(iwp) :: isurf !< index into surf(l) 8725 8673 INTEGER(iwp) :: seldim !< dimension to be incremented 8726 8674 REAL(wp), DIMENSION(2) :: yxorigin !< horizontal copy of origin (y,x) … … 8729 8677 REAL(wp) :: nextdist !< end of current crossing 8730 8678 REAL(wp) :: crmid !< midpoint of crossing 8731 REAL(wp) :: horz_entry !< horizon at entry to column 8732 REAL(wp) :: horz_exit !< horizon at exit from column 8679 REAL(wp) :: horz_full_entry !< full horizon at entry to column 8680 REAL(wp) :: horz_full_exit !< full horizon at exit from column 8681 REAL(wp) :: horz_mixed_entry !< mixed horizon at entry to column 8682 REAL(wp) :: horz_mixed_exit !< mixed horizon at exit from column 8733 8683 REAL(wp) :: bdydim !< boundary for current dimension 8734 8684 REAL(wp), DIMENSION(2) :: crossdist !< distances to boundary for dimensions … … 8738 8688 INTEGER(iwp), DIMENSION(2) :: dimdelta !< dimension direction = +- 1 8739 8689 INTEGER(iwp) :: ip !< number of processor where gridbox reside 8740 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array 8690 INTEGER(iwp) :: ip_last !< previous number of processor where gridbox reside 8691 INTEGER(iwp) :: ig !< 1D index of grid column in global 2D array 8692 INTEGER(iwp) :: ig_last !< 1D index of previous column in global 2D array 8741 8693 INTEGER(iwp) :: maxboxes !< max no of CSF created 8742 8694 INTEGER(iwp) :: nly !< maximum plant canopy height … … 8746 8698 INTEGER(iwp) :: zb1 8747 8699 INTEGER(iwp) :: nz 8748 INTEGER(iwp) :: iz8700 INTEGER(iwp) :: kz 8749 8701 INTEGER(iwp) :: zsgn 8750 8702 INTEGER(iwp) :: lastdir !< wall direction before hitting this column … … 8755 8707 INTEGER(iwp) :: wcount !< RMA window item count 8756 8708 INTEGER(MPI_ADDRESS_KIND) :: wdisp !< RMA window displacement 8709 INTEGER(iwp), DIMENSION(:), ALLOCATABLE:: vert_col !< vertical faces from previous to current column 8710 INTEGER(iwp), DIMENSION(:), ALLOCATABLE:: up_col !< upward oriented surfaces in current column 8711 INTEGER(iwp), DIMENSION(:), ALLOCATABLE:: down_col !< downward oriented surfaces in current column 8757 8712 #endif 8758 8713 … … 8769 8724 yxorigin(:) = origin(2:3) 8770 8725 transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing 8771 horizon = -HUGE(1._wp) 8726 full_horizon = -HUGE(1._wp) 8727 lowest_mixed_ray = nrays 8772 8728 lowest_free_ray = nrays 8773 IF ( rad_angular_discretization .AND. calc_svf) THEN8729 IF ( rad_angular_discretization ) THEN 8774 8730 ALLOCATE(target_surfl(nrays)) 8775 8731 target_surfl(:) = -1 8776 8732 lastdir = -999 8777 8733 lastcolumn(:) = -999 8734 ALLOCATE( vert_col(nz_urban_b:nz_urban_t) ) 8735 ALLOCATE( up_col(nz_urban_b:nz_urban_t) ) 8736 ALLOCATE( down_col(nz_urban_b:nz_urban_t) ) 8778 8737 ENDIF 8779 8738 … … 8806 8765 ENDIF 8807 8766 8767 ip_last = -1 8768 ig_last = -1 8808 8769 lastdist = 0._wp 8809 8770 … … 8839 8800 8840 8801 !-- calculate index of the grid with global indices (column(1),column(2)) 8841 !-- in the array nzterr and plantt and id of the coresponding processor 8842 CALL radiation_calc_global_offset( column(2), column(1), 0, 1, offs_glob=ig ) 8843 8844 IF ( lastdist == 0._wp ) THEN 8845 horz_entry = -HUGE(1._wp) 8802 !-- in the array nzterrt/b and plantt and id of the coresponding processor 8803 CALL radiation_calc_global_offset( column(2), column(1), 0, 1, & 8804 iproc=ip, offs_glob=ig ) 8805 8806 IF ( ip_last < 0 ) THEN 8807 horz_full_entry = -HUGE(1._wp) 8808 horz_mixed_entry = -HUGE(1._wp) 8846 8809 ELSE 8847 horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist 8810 horz_full_entry = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / lastdist 8811 horz_mixed_entry = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / lastdist 8848 8812 ENDIF 8849 horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist 8850 8851 IF ( rad_angular_discretization .AND. calc_svf ) THEN 8852 ! 8853 !-- Identify vertical obstacles hit by rays in current column 8854 DO WHILE ( lowest_free_ray > 0 ) 8855 IF ( zdirs(lowest_free_ray) > horz_entry ) EXIT 8813 horz_full_exit = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / nextdist 8814 horz_mixed_exit = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / nextdist 8815 is_mixed_col = ( nzterrt(ig) /= nzterrb(ig) ) 8816 8817 IF ( rad_angular_discretization ) THEN 8818 ! 8819 !-- Identify vertical full obstacles hit by rays in current column, 8820 !-- mixed rays need to be checked if they are already obstructed. 8821 DO WHILE ( lowest_mixed_ray > lowest_free_ray ) 8822 IF ( zdirs(lowest_mixed_ray) > horz_full_entry ) EXIT 8856 8823 ! 8857 8824 !-- This may only happen after 1st column, so lastdir and lastcolumn are valid 8858 CALL request_itarget(lastdir, & 8859 CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), & 8860 lastcolumn(1), lastcolumn(2), & 8861 target_surfl(lowest_free_ray), target_procs(lowest_free_ray)) 8862 lowest_free_ray = lowest_free_ray - 1 8825 IF ( target_surfl(lowest_mixed_ray) < 0 ) THEN 8826 CALL request_itarget(lastdir, & 8827 CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist), & 8828 lastcolumn(1), lastcolumn(2), & 8829 target_surfl(lowest_mixed_ray)) 8830 ENDIF 8831 lowest_mixed_ray = lowest_mixed_ray - 1 8863 8832 ENDDO 8864 8833 ! 8865 !-- Identify horizontal obstacles hit by rays in current column 8866 DO WHILE ( lowest_free_ray > 0 ) 8867 IF ( zdirs(lowest_free_ray) > horz_exit ) EXIT 8868 CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), & 8869 target_surfl(lowest_free_ray), & 8870 target_procs(lowest_free_ray)) 8871 lowest_free_ray = lowest_free_ray - 1 8834 !-- Identify vertical full obstacles hit by rays in current column, 8835 !-- free rays need no individual checks. 8836 DO WHILE ( lowest_mixed_ray > 0 ) 8837 IF ( zdirs(lowest_mixed_ray) > horz_full_entry ) EXIT 8838 ! 8839 !-- This may only happen after 1st column, so lastdir and lastcolumn are valid 8840 CALL request_itarget(lastdir, & 8841 CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist), & 8842 lastcolumn(1), lastcolumn(2), & 8843 target_surfl(lowest_mixed_ray)) 8844 lowest_mixed_ray = lowest_mixed_ray - 1 8872 8845 ENDDO 8873 ENDIF 8874 8875 horizon = MAX(horizon, horz_entry, horz_exit) 8846 IF ( lowest_free_ray > lowest_mixed_ray ) lowest_free_ray = lowest_mixed_ray 8847 ! 8848 !-- Identify targets for vertical mixed obstacles. 8849 !-- lowest_mixed_ray now points to bottom of vertical mixed obstacles. 8850 IF ( is_mixed_col .AND. ip_last >= 0 ) THEN 8851 ! 8852 !-- Load vertical surfaces belonging to previous column 8853 vert_col(:) = -999 8854 DO isurf = surfg_col_start(ig_last), surfg_col_start(ig_last+1)-1 8855 IF ( surf(id, isurf) == lastdir ) THEN 8856 vert_col(surf(iz, isurf)) = isurf 8857 ENDIF 8858 ENDDO 8859 ! 8860 !-- Previously mixed rays need to be checked whether they are obstructed 8861 DO iray = lowest_mixed_ray, lowest_free_ray-1, -1 8862 IF ( zdirs(iray) > horz_mixed_entry ) EXIT 8863 IF ( target_surfl(iray) >= 0 ) CYCLE 8864 target_surfl(iray) = vert_col(CEILING(-0.5_wp + & 8865 origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface 8866 ENDDO 8867 ! 8868 !-- Previously free rays cannot be obstructed yet 8869 DO iray = lowest_free_ray, 1, -1 8870 IF ( zdirs(iray) > horz_mixed_entry ) THEN 8871 ! 8872 !-- Extend mixed rays by raising the lowest_free ray 8873 !-- (remains unchanged if this is hit in 1st iteration) 8874 lowest_free_ray = iray 8875 EXIT 8876 ENDIF 8877 target_surfl(iray) = vert_col(CEILING(-0.5_wp + & 8878 origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface 8879 ENDDO 8880 ENDIF ! end of mixed horizon 8881 ! 8882 !-- Identify horizontal full obstacles hit by rays in current column, 8883 !-- mixed rays need to be checked if they are already obstructed. 8884 DO WHILE ( lowest_mixed_ray > lowest_free_ray ) 8885 IF ( zdirs(lowest_mixed_ray) > horz_full_exit ) EXIT 8886 IF ( target_surfl(lowest_mixed_ray) < 0 ) THEN 8887 CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), & 8888 target_surfl(lowest_mixed_ray)) 8889 ENDIF 8890 lowest_mixed_ray = lowest_mixed_ray - 1 8891 ENDDO 8892 ! 8893 !-- Identify horizontal full obstacles hit by rays in current column, 8894 !-- free rays need no individual checks. 8895 DO WHILE ( lowest_mixed_ray > 0 ) 8896 IF ( zdirs(lowest_mixed_ray) > horz_full_exit ) EXIT 8897 CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), & 8898 target_surfl(lowest_mixed_ray)) 8899 lowest_mixed_ray = lowest_mixed_ray - 1 8900 ENDDO 8901 IF ( lowest_free_ray > lowest_mixed_ray ) lowest_free_ray = lowest_mixed_ray 8902 ! 8903 !-- Identify targets for horizontal mixed obstacles. 8904 !-- lowest_mixed_ray now points _above_ horizontal full obstacles. 8905 IF ( is_mixed_col ) THEN 8906 ! 8907 !-- Load horizontal surfaces corresponding to current column 8908 up_col(:) = -999 8909 down_col(:) = -999 8910 DO isurf = surfg_col_start(ig), surfg_col_start(ig+1)-1 8911 SELECT CASE ( surf(id, isurf) ) 8912 CASE ( iup ) 8913 up_col(surf(iz, isurf)) = isurf 8914 CASE ( idown ) 8915 down_col(surf(iz, isurf)) = isurf 8916 ENDSELECT 8917 ENDDO 8918 ! 8919 !-- Previously mixed rays need to be checked whether they are obstructed 8920 DO iray = lowest_mixed_ray, lowest_free_ray-1, -1 8921 IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) ) EXIT 8922 IF ( target_surfl(iray) >= 0 ) CYCLE 8923 IF ( zdirs(iray) <= 0._wp ) THEN 8924 ! 8925 !-- Downward pointed ray, cycle k down from entry to exit, 8926 !-- search for upward oriented faces 8927 DO k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), & 8928 CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1 8929 target_surfl(iray) = up_col(k) ! contains -999 if missing surface 8930 IF ( target_surfl(iray) >= 0 ) EXIT 8931 ENDDO 8932 ELSE 8933 ! 8934 !-- Upward pointed ray, cycle k up from entry to exit, 8935 !-- search for downward oriented faces 8936 DO k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), & 8937 FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1 8938 target_surfl(iray) = down_col(k) ! contains -999 if missing surface 8939 IF ( target_surfl(iray) >= 0 ) EXIT 8940 ENDDO 8941 ENDIF 8942 ENDDO 8943 ! 8944 !-- Previously free rays cannot be obstructed yet 8945 DO iray = lowest_free_ray, 1, -1 8946 IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) ) THEN 8947 ! 8948 !-- Extend mixed rays by raising the lowest_free ray 8949 !-- (remains unchanged if this is hit in 1st iteration) 8950 lowest_free_ray = iray 8951 EXIT 8952 ENDIF 8953 IF ( zdirs(iray) <= 0._wp ) THEN 8954 ! 8955 !-- Downward pointed ray, cycle k down from entry to exit, 8956 !-- search upward oriented faces 8957 DO k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), & 8958 CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1 8959 target_surfl(iray) = up_col(k) ! contains -999 if missing surface 8960 IF ( target_surfl(iray) >= 0 ) EXIT 8961 ENDDO 8962 ELSE 8963 ! 8964 !-- Upward pointed ray, cycle k up from entry to exit, 8965 !-- search downward oriented faces 8966 DO k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), & 8967 FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1 8968 target_surfl(iray) = down_col(k) ! contains -999 if missing surface 8969 IF ( target_surfl(iray) >= 0 ) EXIT 8970 ENDDO 8971 ENDIF 8972 ENDDO 8973 ENDIF ! end of mixed horizon 8974 ENDIF ! rad_angular_discretization 8975 full_horizon = MAX(full_horizon, horz_full_entry, horz_full_exit) 8876 8976 8877 8977 IF ( plant_canopy ) THEN … … 8883 8983 IF ( nextdist + eps >= distance ) EXIT 8884 8984 8885 IF ( rad_angular_discretization .AND. calc_svf) THEN8985 IF ( rad_angular_discretization ) THEN 8886 8986 ! 8887 8987 !-- Save wall direction of coming building column (= this air column) 8888 8988 IF ( seldim == 1 ) THEN 8889 8989 IF ( dimdelta(seldim) == 1 ) THEN 8890 lastdir = isouth _u8990 lastdir = isouth 8891 8991 ELSE 8892 lastdir = inorth _u8992 lastdir = inorth 8893 8993 ENDIF 8894 8994 ELSE 8895 8995 IF ( dimdelta(seldim) == 1 ) THEN 8896 lastdir = iwest _u8996 lastdir = iwest 8897 8997 ELSE 8898 lastdir = ieast _u8998 lastdir = ieast 8899 8999 ENDIF 8900 9000 ENDIF 8901 9001 lastcolumn = column 8902 9002 ENDIF 9003 ip_last = ip 9004 ig_last = ig 8903 9005 lastdist = nextdist 8904 9006 dimnext(seldim) = dimnext(seldim) + dimdelta(seldim) … … 8921 9023 !-- For fixed view resolution, we need plant canopy even for rays 8922 9024 !-- to opposing surfaces 8923 lowest_lad = nzterr (ig) + 19025 lowest_lad = nzterrb(ig) + 1 8924 9026 ELSE 8925 9027 ! 8926 !-- We only need LAD for rays directed above horizon (to sky)9028 !-- We only need LAD for rays directed above full horizon (to sky) 8927 9029 lowest_lad = CEILING( -0.5_wp + origin(1) + & 8928 MIN( horizon * rt2_track_dist(i-1), & ! entry8929 horizon * rt2_track_dist(i) ) ) ! exit9030 MIN( full_horizon * rt2_track_dist(i-1), & ! entry 9031 full_horizon * rt2_track_dist(i) ) ) ! exit 8930 9032 ENDIF 8931 9033 ! … … 8970 9072 ENDIF ! plant_canopy 8971 9073 8972 IF ( rad_angular_discretization .AND. calc_svf) THEN9074 IF ( rad_angular_discretization ) THEN 8973 9075 #if defined( __parallel ) 8974 9076 !-- wait for all gridsurf requests to complete … … 8979 9081 ENDIF 8980 9082 #endif 8981 ! 8982 !-- recalculate local surf indices into global ones 8983 DO i = 1, nrays 8984 IF ( target_surfl(i) == -1 ) THEN 8985 itarget(i) = -1 8986 ELSE 8987 itarget(i) = target_surfl(i) + surfstart(target_procs(i)) 8988 ENDIF 9083 itarget(:) = target_surfl(:) 9084 DEALLOCATE( target_surfl ) 9085 ELSE 9086 iray = nrays 9087 DO WHILE ( iray > 0 .AND. zdirs(iray) <= full_horizon ) 9088 itarget(iray) = HUGE(1_iwp) 9089 iray = iray - 1 8989 9090 ENDDO 8990 8991 DEALLOCATE( target_surfl ) 8992 8993 ELSE 8994 itarget(:) = -1 9091 itarget(1:iray) = -1 8995 9092 ENDIF ! rad_angular_discretization 8996 9093 … … 9041 9138 !-- cause less moiree. 9042 9139 IF ( .NOT. rad_angular_discretization ) THEN 9043 IF ( zdirs(k) <= horizon ) CYCLE 9140 ! TODO: remove tradiational discretization (incompatible with full 3D) 9141 IF ( zdirs(k) <= full_horizon ) CYCLE 9044 9142 ENDIF 9045 9143 … … 9052 9150 nz = 2 9053 9151 rt2_dist(nz) = SQRT(dxxyy) 9054 iz = CEILING(-.5_wp + zorig, iwp)9152 kz = CEILING(-.5_wp + zorig, iwp) 9055 9153 ELSE 9056 9154 zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop) … … 9062 9160 qdist = rt2_dist(nz) / (zexit-zorig) 9063 9161 rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/) 9064 iz = zb0 * zsgn9162 kz = zb0 * zsgn 9065 9163 ENDIF 9066 9164 9067 9165 DO l = 2, nz 9068 IF ( rt2_track_lad( iz, i) > 0._wp ) THEN9069 curtrans = exp(-ext_coef * rt2_track_lad( iz, i) * (rt2_dist(l)-rt2_dist(l-1)))9166 IF ( rt2_track_lad(kz, i) > 0._wp ) THEN 9167 curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1))) 9070 9168 9071 9169 IF ( create_csf ) THEN … … 9074 9172 acsf(ncsfl)%itx = rt2_track(2,i) 9075 9173 acsf(ncsfl)%ity = rt2_track(1,i) 9076 acsf(ncsfl)%itz = iz9174 acsf(ncsfl)%itz = kz 9077 9175 acsf(ncsfl)%isurfs = iorig 9078 9176 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k) … … 9081 9179 transparency(k) = transparency(k) * curtrans 9082 9180 ENDIF 9083 iz = iz + zsgn9181 kz = kz + zsgn 9084 9182 ENDDO ! l = 1, nz - 1 9085 9183 ENDDO ! k = 1, nrays 9086 9184 ENDDO ! i = 1, ntrack 9087 9088 transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778) 9185 ! 9186 !-- Reset rays above horizon to transparent (see NOTE 6778) 9187 WHERE( itarget < 0 ) transparency = 1.0_wp 9089 9188 ENDIF 9090 9189 … … 9098 9197 ! 9099 9198 !-- See NOTE 6778 above 9100 IF ( zdirs(k) <= horizon) CYCLE9199 IF ( itarget(k) >= 0 ) CYCLE 9101 9200 9102 9201 zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1) … … 9108 9207 nz = 2 9109 9208 rt2_dist(nz) = SQRT(dxxyy) 9110 iz = NINT(zexit, iwp)9209 kz = NINT(zexit, iwp) 9111 9210 ELSE 9112 9211 zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop) … … 9118 9217 qdist = rt2_dist(nz) / (zexit-zorig) 9119 9218 rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/) 9120 iz = zb0 * zsgn9219 kz = zb0 * zsgn 9121 9220 ENDIF 9122 9221 9123 9222 DO l = 2, nz 9124 IF ( rt2_track_lad( iz, i) > 0._wp ) THEN9125 curtrans = exp(-ext_coef * rt2_track_lad( iz, i) * (rt2_dist(l)-rt2_dist(l-1)))9223 IF ( rt2_track_lad(kz, i) > 0._wp ) THEN 9224 curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1))) 9126 9225 9127 9226 IF ( create_csf ) THEN … … 9130 9229 acsf(ncsfl)%itx = rt2_track(2,i) 9131 9230 acsf(ncsfl)%ity = rt2_track(1,i) 9132 acsf(ncsfl)%itz = iz 9133 IF ( itarget(k) /= -1 ) STOP 1 !FIXME remove after test 9231 acsf(ncsfl)%itz = kz 9134 9232 acsf(ncsfl)%isurfs = -1 9135 9233 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k) … … 9138 9236 transparency(k) = transparency(k) * curtrans 9139 9237 ENDIF 9140 iz = iz + zsgn9238 kz = kz + zsgn 9141 9239 ENDDO ! l = 1, nz - 1 9142 9240 ENDDO ! k = 1, nrays … … 9144 9242 ENDIF ! plant_canopy 9145 9243 9146 IF ( .NOT. (rad_angular_discretization .AND. calc_svf) ) THEN9147 !9148 !-- Just update lowest_free_ray according to horizon9149 DO WHILE ( lowest_free_ray > 0 )9150 IF ( zdirs(lowest_free_ray) > horizon ) EXIT9151 lowest_free_ray = lowest_free_ray - 19152 ENDDO9153 ENDIF9154 9155 9244 CONTAINS 9156 9245 9157 SUBROUTINE request_itarget( d, z, y, x, isurfl , iproc)9246 SUBROUTINE request_itarget( d, z, y, x, isurfl ) 9158 9247 9159 9248 INTEGER(iwp), INTENT(in) :: d, z, y, x 9160 9249 INTEGER(iwp), TARGET, INTENT(out) :: isurfl 9161 INTEGER(iwp) , INTENT(out):: iproc9250 INTEGER(iwp) :: iproc 9162 9251 9163 9252 #if defined( __parallel ) … … 9181 9270 !-- set index target_surfl(i) 9182 9271 isurfl = gridsurf(d,z,y,x) 9183 iproc = 0 ! required to avoid compile error about unused variable in serial mode9184 9272 #endif 9185 9273 … … 9299 9387 9300 9388 !-- first check: are the two surfaces directed in the same direction 9301 IF ( (d==iup_u .OR. d==iup_l ) & 9302 .AND. (d2==iup_u .OR. d2==iup_l) ) RETURN 9303 IF ( (d==isouth_u .OR. d==isouth_l ) & 9304 .AND. (d2==isouth_u .OR. d2==isouth_l) ) RETURN 9305 IF ( (d==inorth_u .OR. d==inorth_l ) & 9306 .AND. (d2==inorth_u .OR. d2==inorth_l) ) RETURN 9307 IF ( (d==iwest_u .OR. d==iwest_l ) & 9308 .AND. (d2==iwest_u .OR. d2==iwest_l ) ) RETURN 9309 IF ( (d==ieast_u .OR. d==ieast_l ) & 9310 .AND. (d2==ieast_u .OR. d2==ieast_l ) ) RETURN 9389 IF ( d==iup .AND. d2==iup ) RETURN 9390 IF ( d==isouth .AND. d2==isouth ) RETURN 9391 IF ( d==inorth .AND. d2==inorth ) RETURN 9392 IF ( d==iwest .AND. d2==iwest ) RETURN 9393 IF ( d==ieast .AND. d2==ieast ) RETURN 9311 9394 9312 9395 !-- second check: are surfaces facing away from each other 9313 9396 SELECT CASE (d) 9314 CASE (iup _u, iup_l)!< upward facing surfaces9397 CASE (iup) !< upward facing surfaces 9315 9398 IF ( z2 < z ) RETURN 9316 CASE (isouth _u, isouth_l) !< southward facing surfaces9399 CASE (isouth) !< southward facing surfaces 9317 9400 IF ( y2 > y ) RETURN 9318 CASE (inorth _u, inorth_l) !< northward facing surfaces9401 CASE (inorth) !< northward facing surfaces 9319 9402 IF ( y2 < y ) RETURN 9320 CASE (iwest _u, iwest_l)!< westward facing surfaces9403 CASE (iwest) !< westward facing surfaces 9321 9404 IF ( x2 > x ) RETURN 9322 CASE (ieast _u, ieast_l)!< eastward facing surfaces9405 CASE (ieast) !< eastward facing surfaces 9323 9406 IF ( x2 < x ) RETURN 9324 9407 END SELECT 9325 9408 9326 9409 SELECT CASE (d2) 9327 CASE (iup _u)!< ground, roof9410 CASE (iup) !< ground, roof 9328 9411 IF ( z < z2 ) RETURN 9329 CASE (isouth _u, isouth_l) !< south facing9412 CASE (isouth) !< south facing 9330 9413 IF ( y > y2 ) RETURN 9331 CASE (inorth _u, inorth_l) !< north facing9414 CASE (inorth) !< north facing 9332 9415 IF ( y < y2 ) RETURN 9333 CASE (iwest _u, iwest_l)!< west facing9416 CASE (iwest) !< west facing 9334 9417 IF ( x > x2 ) RETURN 9335 CASE (ieast _u, ieast_l)!< east facing9418 CASE (ieast) !< east facing 9336 9419 IF ( x < x2 ) RETURN 9337 9420 CASE (-1) … … 9938 10021 INTEGER(iwp) :: l, m !< index of current surface element 9939 10022 9940 INTEGER(iwp) :: ids, idsint _u, idsint_l, isurf10023 INTEGER(iwp) :: ids, idsint, isurf 9941 10024 CHARACTER(LEN=varnamelength) :: var 9942 10025 … … 9951 10034 IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN 9952 10035 ids = i 9953 idsint_u = dirint_u(ids) 9954 idsint_l = dirint_l(ids) 10036 idsint = dirint(ids) 9955 10037 var = var(:k-j) 9956 10038 EXIT … … 10399 10481 CASE ( 'rtm_rad_net' ) 10400 10482 !-- array of complete radiation balance 10401 DO isurf = dirstart(ids), dirend(ids)10402 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10483 DO isurf = 1, nsurfl 10484 IF ( surfl(id,isurf) == idsint ) THEN 10403 10485 surfradnet_av(isurf) = surfradnet_av(isurf) + & 10404 10486 surfinsw(isurf) - surfoutsw(isurf) + & … … 10409 10491 CASE ( 'rtm_rad_insw' ) 10410 10492 !-- array of sw radiation falling to surface after i-th reflection 10411 DO isurf = dirstart(ids), dirend(ids)10412 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10493 DO isurf = 1, nsurfl 10494 IF ( surfl(id,isurf) == idsint ) THEN 10413 10495 surfinsw_av(isurf) = surfinsw_av(isurf) + surfinsw(isurf) 10414 10496 ENDIF … … 10417 10499 CASE ( 'rtm_rad_inlw' ) 10418 10500 !-- array of lw radiation falling to surface after i-th reflection 10419 DO isurf = dirstart(ids), dirend(ids)10420 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10501 DO isurf = 1, nsurfl 10502 IF ( surfl(id,isurf) == idsint ) THEN 10421 10503 surfinlw_av(isurf) = surfinlw_av(isurf) + surfinlw(isurf) 10422 10504 ENDIF … … 10425 10507 CASE ( 'rtm_rad_inswdir' ) 10426 10508 !-- array of direct sw radiation falling to surface from sun 10427 DO isurf = dirstart(ids), dirend(ids)10428 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10509 DO isurf = 1, nsurfl 10510 IF ( surfl(id,isurf) == idsint ) THEN 10429 10511 surfinswdir_av(isurf) = surfinswdir_av(isurf) + surfinswdir(isurf) 10430 10512 ENDIF … … 10433 10515 CASE ( 'rtm_rad_inswdif' ) 10434 10516 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 10435 DO isurf = dirstart(ids), dirend(ids)10436 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10517 DO isurf = 1, nsurfl 10518 IF ( surfl(id,isurf) == idsint ) THEN 10437 10519 surfinswdif_av(isurf) = surfinswdif_av(isurf) + surfinswdif(isurf) 10438 10520 ENDIF … … 10441 10523 CASE ( 'rtm_rad_inswref' ) 10442 10524 !-- array of sw radiation falling to surface from reflections 10443 DO isurf = dirstart(ids), dirend(ids)10444 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10525 DO isurf = 1, nsurfl 10526 IF ( surfl(id,isurf) == idsint ) THEN 10445 10527 surfinswref_av(isurf) = surfinswref_av(isurf) + surfinsw(isurf) - & 10446 10528 surfinswdir(isurf) - surfinswdif(isurf) … … 10451 10533 CASE ( 'rtm_rad_inlwdif' ) 10452 10534 !-- array of sw radiation falling to surface after i-th reflection 10453 DO isurf = dirstart(ids), dirend(ids)10454 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10535 DO isurf = 1, nsurfl 10536 IF ( surfl(id,isurf) == idsint ) THEN 10455 10537 surfinlwdif_av(isurf) = surfinlwdif_av(isurf) + surfinlwdif(isurf) 10456 10538 ENDIF … … 10459 10541 CASE ( 'rtm_rad_inlwref' ) 10460 10542 !-- array of lw radiation falling to surface from reflections 10461 DO isurf = dirstart(ids), dirend(ids)10462 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10543 DO isurf = 1, nsurfl 10544 IF ( surfl(id,isurf) == idsint ) THEN 10463 10545 surfinlwref_av(isurf) = surfinlwref_av(isurf) + & 10464 10546 surfinlw(isurf) - surfinlwdif(isurf) … … 10468 10550 CASE ( 'rtm_rad_outsw' ) 10469 10551 !-- array of sw radiation emitted from surface after i-th reflection 10470 DO isurf = dirstart(ids), dirend(ids)10471 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10552 DO isurf = 1, nsurfl 10553 IF ( surfl(id,isurf) == idsint ) THEN 10472 10554 surfoutsw_av(isurf) = surfoutsw_av(isurf) + surfoutsw(isurf) 10473 10555 ENDIF … … 10476 10558 CASE ( 'rtm_rad_outlw' ) 10477 10559 !-- array of lw radiation emitted from surface after i-th reflection 10478 DO isurf = dirstart(ids), dirend(ids)10479 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10560 DO isurf = 1, nsurfl 10561 IF ( surfl(id,isurf) == idsint ) THEN 10480 10562 surfoutlw_av(isurf) = surfoutlw_av(isurf) + surfoutlw(isurf) 10481 10563 ENDIF … … 10484 10566 CASE ( 'rtm_rad_ressw' ) 10485 10567 !-- array of residua of sw radiation absorbed in surface after last reflection 10486 DO isurf = dirstart(ids), dirend(ids)10487 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10568 DO isurf = 1, nsurfl 10569 IF ( surfl(id,isurf) == idsint ) THEN 10488 10570 surfins_av(isurf) = surfins_av(isurf) + surfins(isurf) 10489 10571 ENDIF … … 10492 10574 CASE ( 'rtm_rad_reslw' ) 10493 10575 !-- array of residua of lw radiation absorbed in surface after last reflection 10494 DO isurf = dirstart(ids), dirend(ids)10495 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10576 DO isurf = 1, nsurfl 10577 IF ( surfl(id,isurf) == idsint ) THEN 10496 10578 surfinl_av(isurf) = surfinl_av(isurf) + surfinl(isurf) 10497 10579 ENDIF … … 10696 10778 CASE ( 'rtm_rad_net' ) 10697 10779 !-- array of complete radiation balance 10698 DO isurf = dirstart(ids), dirend(ids)10699 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10780 DO isurf = 1, nsurfl 10781 IF ( surfl(id,isurf) == idsint ) THEN 10700 10782 surfradnet_av(isurf) = surfradnet_av(isurf) / REAL( average_count_3d, kind=wp ) 10701 10783 ENDIF … … 10704 10786 CASE ( 'rtm_rad_insw' ) 10705 10787 !-- array of sw radiation falling to surface after i-th reflection 10706 DO isurf = dirstart(ids), dirend(ids)10707 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10788 DO isurf = 1, nsurfl 10789 IF ( surfl(id,isurf) == idsint ) THEN 10708 10790 surfinsw_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp ) 10709 10791 ENDIF … … 10712 10794 CASE ( 'rtm_rad_inlw' ) 10713 10795 !-- array of lw radiation falling to surface after i-th reflection 10714 DO isurf = dirstart(ids), dirend(ids)10715 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10796 DO isurf = 1, nsurfl 10797 IF ( surfl(id,isurf) == idsint ) THEN 10716 10798 surfinlw_av(isurf) = surfinlw_av(isurf) / REAL( average_count_3d, kind=wp ) 10717 10799 ENDIF … … 10720 10802 CASE ( 'rtm_rad_inswdir' ) 10721 10803 !-- array of direct sw radiation falling to surface from sun 10722 DO isurf = dirstart(ids), dirend(ids)10723 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10804 DO isurf = 1, nsurfl 10805 IF ( surfl(id,isurf) == idsint ) THEN 10724 10806 surfinswdir_av(isurf) = surfinswdir_av(isurf) / REAL( average_count_3d, kind=wp ) 10725 10807 ENDIF … … 10728 10810 CASE ( 'rtm_rad_inswdif' ) 10729 10811 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 10730 DO isurf = dirstart(ids), dirend(ids)10731 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10812 DO isurf = 1, nsurfl 10813 IF ( surfl(id,isurf) == idsint ) THEN 10732 10814 surfinswdif_av(isurf) = surfinswdif_av(isurf) / REAL( average_count_3d, kind=wp ) 10733 10815 ENDIF … … 10736 10818 CASE ( 'rtm_rad_inswref' ) 10737 10819 !-- array of sw radiation falling to surface from reflections 10738 DO isurf = dirstart(ids), dirend(ids)10739 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10820 DO isurf = 1, nsurfl 10821 IF ( surfl(id,isurf) == idsint ) THEN 10740 10822 surfinswref_av(isurf) = surfinswref_av(isurf) / REAL( average_count_3d, kind=wp ) 10741 10823 ENDIF … … 10744 10826 CASE ( 'rtm_rad_inlwdif' ) 10745 10827 !-- array of sw radiation falling to surface after i-th reflection 10746 DO isurf = dirstart(ids), dirend(ids)10747 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10828 DO isurf = 1, nsurfl 10829 IF ( surfl(id,isurf) == idsint ) THEN 10748 10830 surfinlwdif_av(isurf) = surfinlwdif_av(isurf) / REAL( average_count_3d, kind=wp ) 10749 10831 ENDIF … … 10752 10834 CASE ( 'rtm_rad_inlwref' ) 10753 10835 !-- array of lw radiation falling to surface from reflections 10754 DO isurf = dirstart(ids), dirend(ids)10755 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10836 DO isurf = 1, nsurfl 10837 IF ( surfl(id,isurf) == idsint ) THEN 10756 10838 surfinlwref_av(isurf) = surfinlwref_av(isurf) / REAL( average_count_3d, kind=wp ) 10757 10839 ENDIF … … 10760 10842 CASE ( 'rtm_rad_outsw' ) 10761 10843 !-- array of sw radiation emitted from surface after i-th reflection 10762 DO isurf = dirstart(ids), dirend(ids)10763 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10844 DO isurf = 1, nsurfl 10845 IF ( surfl(id,isurf) == idsint ) THEN 10764 10846 surfoutsw_av(isurf) = surfoutsw_av(isurf) / REAL( average_count_3d, kind=wp ) 10765 10847 ENDIF … … 10768 10850 CASE ( 'rtm_rad_outlw' ) 10769 10851 !-- array of lw radiation emitted from surface after i-th reflection 10770 DO isurf = dirstart(ids), dirend(ids)10771 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10852 DO isurf = 1, nsurfl 10853 IF ( surfl(id,isurf) == idsint ) THEN 10772 10854 surfoutlw_av(isurf) = surfoutlw_av(isurf) / REAL( average_count_3d, kind=wp ) 10773 10855 ENDIF … … 10776 10858 CASE ( 'rtm_rad_ressw' ) 10777 10859 !-- array of residua of sw radiation absorbed in surface after last reflection 10778 DO isurf = dirstart(ids), dirend(ids)10779 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10860 DO isurf = 1, nsurfl 10861 IF ( surfl(id,isurf) == idsint ) THEN 10780 10862 surfins_av(isurf) = surfins_av(isurf) / REAL( average_count_3d, kind=wp ) 10781 10863 ENDIF … … 10784 10866 CASE ( 'rtm_rad_reslw' ) 10785 10867 !-- array of residua of lw radiation absorbed in surface after last reflection 10786 DO isurf = dirstart(ids), dirend(ids)10787 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN10868 DO isurf = 1, nsurfl 10869 IF ( surfl(id,isurf) == idsint ) THEN 10788 10870 surfinl_av(isurf) = surfinl_av(isurf) / REAL( average_count_3d, kind=wp ) 10789 10871 ENDIF … … 11354 11436 11355 11437 CHARACTER (len=varnamelength) :: var, surfid 11356 INTEGER(iwp) :: ids,idsint _u,idsint_l,isurf,isvf,isurfs,isurflt,ipcgb11438 INTEGER(iwp) :: ids,idsint,isurf,isvf,isurfs,isurflt,ipcgb 11357 11439 INTEGER(iwp) :: is, js, ks, istat 11358 11440 … … 11378 11460 IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN 11379 11461 ids = i 11380 idsint_u = dirint_u(ids) 11381 idsint_l = dirint_l(ids) 11462 idsint = dirint(ids) 11382 11463 var = var(:k-j) 11383 11464 EXIT … … 11596 11677 CASE ( 'rtm_rad_net' ) 11597 11678 !-- array of complete radiation balance 11598 DO isurf = dirstart(ids), dirend(ids)11599 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11679 DO isurf = 1, nsurfl 11680 IF ( surfl(id,isurf) == idsint ) THEN 11600 11681 IF ( av == 0 ) THEN 11601 11682 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = & … … 11609 11690 CASE ( 'rtm_rad_insw' ) 11610 11691 !-- array of sw radiation falling to surface after i-th reflection 11611 DO isurf = dirstart(ids), dirend(ids)11612 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11692 DO isurf = 1, nsurfl 11693 IF ( surfl(id,isurf) == idsint ) THEN 11613 11694 IF ( av == 0 ) THEN 11614 11695 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw(isurf) … … 11621 11702 CASE ( 'rtm_rad_inlw' ) 11622 11703 !-- array of lw radiation falling to surface after i-th reflection 11623 DO isurf = dirstart(ids), dirend(ids)11624 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11704 DO isurf = 1, nsurfl 11705 IF ( surfl(id,isurf) == idsint ) THEN 11625 11706 IF ( av == 0 ) THEN 11626 11707 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) … … 11633 11714 CASE ( 'rtm_rad_inswdir' ) 11634 11715 !-- array of direct sw radiation falling to surface from sun 11635 DO isurf = dirstart(ids), dirend(ids)11636 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11716 DO isurf = 1, nsurfl 11717 IF ( surfl(id,isurf) == idsint ) THEN 11637 11718 IF ( av == 0 ) THEN 11638 11719 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir(isurf) … … 11645 11726 CASE ( 'rtm_rad_inswdif' ) 11646 11727 !-- array of difusion sw radiation falling to surface from sky and borders of the domain 11647 DO isurf = dirstart(ids), dirend(ids)11648 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11728 DO isurf = 1, nsurfl 11729 IF ( surfl(id,isurf) == idsint ) THEN 11649 11730 IF ( av == 0 ) THEN 11650 11731 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif(isurf) … … 11657 11738 CASE ( 'rtm_rad_inswref' ) 11658 11739 !-- array of sw radiation falling to surface from reflections 11659 DO isurf = dirstart(ids), dirend(ids)11660 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11740 DO isurf = 1, nsurfl 11741 IF ( surfl(id,isurf) == idsint ) THEN 11661 11742 IF ( av == 0 ) THEN 11662 11743 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = & … … 11670 11751 CASE ( 'rtm_rad_inlwdif' ) 11671 11752 !-- array of difusion lw radiation falling to surface from sky and borders of the domain 11672 DO isurf = dirstart(ids), dirend(ids)11673 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11753 DO isurf = 1, nsurfl 11754 IF ( surfl(id,isurf) == idsint ) THEN 11674 11755 IF ( av == 0 ) THEN 11675 11756 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif(isurf) … … 11682 11763 CASE ( 'rtm_rad_inlwref' ) 11683 11764 !-- array of lw radiation falling to surface from reflections 11684 DO isurf = dirstart(ids), dirend(ids)11685 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11765 DO isurf = 1, nsurfl 11766 IF ( surfl(id,isurf) == idsint ) THEN 11686 11767 IF ( av == 0 ) THEN 11687 11768 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) - surfinlwdif(isurf) … … 11694 11775 CASE ( 'rtm_rad_outsw' ) 11695 11776 !-- array of sw radiation emitted from surface after i-th reflection 11696 DO isurf = dirstart(ids), dirend(ids)11697 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11777 DO isurf = 1, nsurfl 11778 IF ( surfl(id,isurf) == idsint ) THEN 11698 11779 IF ( av == 0 ) THEN 11699 11780 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw(isurf) … … 11706 11787 CASE ( 'rtm_rad_outlw' ) 11707 11788 !-- array of lw radiation emitted from surface after i-th reflection 11708 DO isurf = dirstart(ids), dirend(ids)11709 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11789 DO isurf = 1, nsurfl 11790 IF ( surfl(id,isurf) == idsint ) THEN 11710 11791 IF ( av == 0 ) THEN 11711 11792 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw(isurf) … … 11718 11799 CASE ( 'rtm_rad_ressw' ) 11719 11800 !-- average of array of residua of sw radiation absorbed in surface after last reflection 11720 DO isurf = dirstart(ids), dirend(ids)11721 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11801 DO isurf = 1, nsurfl 11802 IF ( surfl(id,isurf) == idsint ) THEN 11722 11803 IF ( av == 0 ) THEN 11723 11804 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins(isurf) … … 11730 11811 CASE ( 'rtm_rad_reslw' ) 11731 11812 !-- average of array of residua of lw radiation absorbed in surface after last reflection 11732 DO isurf = dirstart(ids), dirend(ids)11733 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11813 DO isurf = 1, nsurfl 11814 IF ( surfl(id,isurf) == idsint ) THEN 11734 11815 IF ( av == 0 ) THEN 11735 11816 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl(isurf) … … 11838 11919 ! 11839 11920 !-- sky view factor 11840 DO isurf = dirstart(ids), dirend(ids)11841 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11921 DO isurf = 1, nsurfl 11922 IF ( surfl(id,isurf) == idsint ) THEN 11842 11923 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvf(isurf) 11843 11924 ENDIF … … 11847 11928 ! 11848 11929 !-- sky view factor 11849 DO isurf = dirstart(ids), dirend(ids)11850 IF ( surfl(id,isurf) == ids ) THEN11930 DO isurf = 1, nsurfl 11931 IF ( surfl(id,isurf) == idsint ) THEN 11851 11932 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvft(isurf) 11852 11933 ENDIF … … 11866 11947 11867 11948 IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. surf(iz,isurfs) == ks .AND. & 11868 (surfl(id,isurflt) == idsint_u .OR. surfl(id,isurflt) == idsint_l )) THEN11949 surfl(id,isurflt) == idsint ) THEN 11869 11950 ! 11870 11951 !-- correct source surface … … 11876 11957 ! 11877 11958 !-- surface albedo 11878 DO isurf = dirstart(ids), dirend(ids)11879 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11959 DO isurf = 1, nsurfl 11960 IF ( surfl(id,isurf) == idsint ) THEN 11880 11961 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = albedo_surf(isurf) 11881 11962 ENDIF … … 11885 11966 ! 11886 11967 !-- surface emissivity, weighted average 11887 DO isurf = dirstart(ids), dirend(ids)11888 IF ( surfl(id,isurf) == idsint _u .OR. surfl(id,isurf) == idsint_l) THEN11968 DO isurf = 1, nsurfl 11969 IF ( surfl(id,isurf) == idsint ) THEN 11889 11970 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = emiss_surf(isurf) 11890 11971 ENDIF -
palm/trunk/SOURCE/urban_surface_mod.f90
r4630 r4653 362 362 force_radiation_call, & 363 363 id, & 364 ieast_l, & 365 ieast_u, & 366 inorth_l, & 367 inorth_u, & 368 isouth_l, & 369 isouth_u, & 370 iup_l, & 371 iup_u, & 372 iwest_l, & 373 iwest_u, & 364 ieast, & 365 inorth, & 366 isouth, & 367 iup, & 368 iwest, & 374 369 nz_urban_b, & 375 370 nz_urban_t, & … … 1343 1338 INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions 1344 1339 CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /) 1345 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup _u, isouth_u, inorth_u, iwest_u, ieast_u/)1340 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup, isouth, inorth, iwest, ieast /) 1346 1341 1347 1342 … … 2460 2455 2461 2456 INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions 2462 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u/) !<2457 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup, isouth, inorth, iwest, ieast /) !< 2463 2458 INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: diridx = (/ -1, 1, 0, 3, 2 /) 2464 2459 !< index for surf_*_v: 0:3 = (North, South, East, West) … … 2531 2526 ! 2532 2527 !-- Array of surface height (z) 2533 IF ( idsint == iup _u) THEN2528 IF ( idsint == iup ) THEN 2534 2529 DO m = 1, surf_usm_h%ns 2535 2530 i = surf_usm_h%i(m) … … 2551 2546 ! 2552 2547 !-- Surface category 2553 IF ( idsint == iup _u) THEN2548 IF ( idsint == iup ) THEN 2554 2549 DO m = 1, surf_usm_h%ns 2555 2550 i = surf_usm_h%i(m) … … 2571 2566 ! 2572 2567 !-- Transmissivity window tiles 2573 IF ( idsint == iup _u) THEN2568 IF ( idsint == iup ) THEN 2574 2569 DO m = 1, surf_usm_h%ns 2575 2570 i = surf_usm_h%i(m) … … 2592 2587 !-- Array of sensible heat flux from surfaces 2593 2588 IF ( av == 0 ) THEN 2594 IF ( idsint == iup _u) THEN2589 IF ( idsint == iup ) THEN 2595 2590 DO m = 1, surf_usm_h%ns 2596 2591 i = surf_usm_h%i(m) … … 2609 2604 ENDIF 2610 2605 ELSE 2611 IF ( idsint == iup _u) THEN2606 IF ( idsint == iup ) THEN 2612 2607 DO m = 1, surf_usm_h%ns 2613 2608 i = surf_usm_h%i(m) … … 2632 2627 !-- Array of latent heat flux from surfaces 2633 2628 IF ( av == 0 ) THEN 2634 IF ( idsint == iup _u) THEN2629 IF ( idsint == iup ) THEN 2635 2630 DO m = 1, surf_usm_h%ns 2636 2631 i = surf_usm_h%i(m) … … 2649 2644 ENDIF 2650 2645 ELSE 2651 IF ( idsint == iup _u) THEN2646 IF ( idsint == iup ) THEN 2652 2647 DO m = 1, surf_usm_h%ns 2653 2648 i = surf_usm_h%i(m) … … 2671 2666 !-- Array of latent heat flux from vegetation surfaces 2672 2667 IF ( av == 0 ) THEN 2673 IF ( idsint == iup _u) THEN2668 IF ( idsint == iup ) THEN 2674 2669 DO m = 1, surf_usm_h%ns 2675 2670 i = surf_usm_h%i(m) … … 2688 2683 ENDIF 2689 2684 ELSE 2690 IF ( idsint == iup _u) THEN2685 IF ( idsint == iup ) THEN 2691 2686 DO m = 1, surf_usm_h%ns 2692 2687 i = surf_usm_h%i(m) … … 2710 2705 !-- Array of latent heat flux from surfaces with liquid 2711 2706 IF ( av == 0 ) THEN 2712 IF ( idsint == iup _u) THEN2707 IF ( idsint == iup ) THEN 2713 2708 DO m = 1, surf_usm_h%ns 2714 2709 i = surf_usm_h%i(m) … … 2727 2722 ENDIF 2728 2723 ELSE 2729 IF ( idsint == iup _u) THEN2724 IF ( idsint == iup ) THEN 2730 2725 DO m = 1, surf_usm_h%ns 2731 2726 i = surf_usm_h%i(m) … … 2749 2744 !-- Array of heat flux from ground (land, wall, roof) 2750 2745 IF ( av == 0 ) THEN 2751 IF ( idsint == iup _u) THEN2746 IF ( idsint == iup ) THEN 2752 2747 DO m = 1, surf_usm_h%ns 2753 2748 i = surf_usm_h%i(m) … … 2766 2761 ENDIF 2767 2762 ELSE 2768 IF ( idsint == iup _u) THEN2763 IF ( idsint == iup ) THEN 2769 2764 DO m = 1, surf_usm_h%ns 2770 2765 i = surf_usm_h%i(m) … … 2788 2783 !-- Array of heat flux from window ground (land, wall, roof) 2789 2784 IF ( av == 0 ) THEN 2790 IF ( idsint == iup _u) THEN2785 IF ( idsint == iup ) THEN 2791 2786 DO m = 1, surf_usm_h%ns 2792 2787 i = surf_usm_h%i(m) … … 2805 2800 ENDIF 2806 2801 ELSE 2807 IF ( idsint == iup _u) THEN2802 IF ( idsint == iup ) THEN 2808 2803 DO m = 1, surf_usm_h%ns 2809 2804 i = surf_usm_h%i(m) … … 2827 2822 !-- Array of heat flux from green ground (land, wall, roof) 2828 2823 IF ( av == 0 ) THEN 2829 IF ( idsint == iup _u) THEN2824 IF ( idsint == iup ) THEN 2830 2825 DO m = 1, surf_usm_h%ns 2831 2826 i = surf_usm_h%i(m) … … 2844 2839 ENDIF 2845 2840 ELSE 2846 IF ( idsint == iup _u) THEN2841 IF ( idsint == iup ) THEN 2847 2842 DO m = 1, surf_usm_h%ns 2848 2843 i = surf_usm_h%i(m) … … 2866 2861 !-- Array of heat flux from indoor ground (land, wall, roof) 2867 2862 IF ( av == 0 ) THEN 2868 IF ( idsint == iup _u) THEN2863 IF ( idsint == iup ) THEN 2869 2864 DO m = 1, surf_usm_h%ns 2870 2865 i = surf_usm_h%i(m) … … 2883 2878 ENDIF 2884 2879 ELSE 2885 IF ( idsint == iup _u) THEN2880 IF ( idsint == iup ) THEN 2886 2881 DO m = 1, surf_usm_h%ns 2887 2882 i = surf_usm_h%i(m) … … 2905 2900 !-- Array of heat flux from indoor window ground (land, wall, roof) 2906 2901 IF ( av == 0 ) THEN 2907 IF ( idsint == iup _u) THEN2902 IF ( idsint == iup ) THEN 2908 2903 DO m = 1, surf_usm_h%ns 2909 2904 i = surf_usm_h%i(m) … … 2922 2917 ENDIF 2923 2918 ELSE 2924 IF ( idsint == iup _u) THEN2919 IF ( idsint == iup ) THEN 2925 2920 DO m = 1, surf_usm_h%ns 2926 2921 i = surf_usm_h%i(m) … … 2944 2939 !-- Surface temperature for surfaces 2945 2940 IF ( av == 0 ) THEN 2946 IF ( idsint == iup _u) THEN2941 IF ( idsint == iup ) THEN 2947 2942 DO m = 1, surf_usm_h%ns 2948 2943 i = surf_usm_h%i(m) … … 2961 2956 ENDIF 2962 2957 ELSE 2963 IF ( idsint == iup _u) THEN2958 IF ( idsint == iup ) THEN 2964 2959 DO m = 1, surf_usm_h%ns 2965 2960 i = surf_usm_h%i(m) … … 2983 2978 !-- Surface temperature for window surfaces 2984 2979 IF ( av == 0 ) THEN 2985 IF ( idsint == iup _u) THEN2980 IF ( idsint == iup ) THEN 2986 2981 DO m = 1, surf_usm_h%ns 2987 2982 i = surf_usm_h%i(m) … … 3001 2996 3002 2997 ELSE 3003 IF ( idsint == iup _u) THEN2998 IF ( idsint == iup ) THEN 3004 2999 DO m = 1, surf_usm_h%ns 3005 3000 i = surf_usm_h%i(m) … … 3025 3020 !-- Surface temperature for green surfaces 3026 3021 IF ( av == 0 ) THEN 3027 IF ( idsint == iup _u) THEN3022 IF ( idsint == iup ) THEN 3028 3023 DO m = 1, surf_usm_h%ns 3029 3024 i = surf_usm_h%i(m) … … 3043 3038 3044 3039 ELSE 3045 IF ( idsint == iup _u) THEN3040 IF ( idsint == iup ) THEN 3046 3041 DO m = 1, surf_usm_h%ns 3047 3042 i = surf_usm_h%i(m) … … 3067 3062 !-- Near surface temperature for whole surfaces 3068 3063 IF ( av == 0 ) THEN 3069 IF ( idsint == iup _u) THEN3064 IF ( idsint == iup ) THEN 3070 3065 DO m = 1, surf_usm_h%ns 3071 3066 i = surf_usm_h%i(m) … … 3086 3081 3087 3082 ELSE 3088 IF ( idsint == iup _u) THEN3083 IF ( idsint == iup ) THEN 3089 3084 DO m = 1, surf_usm_h%ns 3090 3085 i = surf_usm_h%i(m) … … 3109 3104 !-- Wall temperature for iwl layer of walls and land 3110 3105 IF ( av == 0 ) THEN 3111 IF ( idsint == iup _u) THEN3106 IF ( idsint == iup ) THEN 3112 3107 DO m = 1, surf_usm_h%ns 3113 3108 i = surf_usm_h%i(m) … … 3126 3121 ENDIF 3127 3122 ELSE 3128 IF ( idsint == iup _u) THEN3123 IF ( idsint == iup ) THEN 3129 3124 DO m = 1, surf_usm_h%ns 3130 3125 i = surf_usm_h%i(m) … … 3148 3143 !-- Window temperature for iwl layer of walls and land 3149 3144 IF ( av == 0 ) THEN 3150 IF ( idsint == iup _u) THEN3145 IF ( idsint == iup ) THEN 3151 3146 DO m = 1, surf_usm_h%ns 3152 3147 i = surf_usm_h%i(m) … … 3165 3160 ENDIF 3166 3161 ELSE 3167 IF ( idsint == iup _u) THEN3162 IF ( idsint == iup ) THEN 3168 3163 DO m = 1, surf_usm_h%ns 3169 3164 i = surf_usm_h%i(m) … … 3187 3182 !-- Green temperature for iwl layer of walls and land 3188 3183 IF ( av == 0 ) THEN 3189 IF ( idsint == iup _u) THEN3184 IF ( idsint == iup ) THEN 3190 3185 DO m = 1, surf_usm_h%ns 3191 3186 i = surf_usm_h%i(m) … … 3204 3199 ENDIF 3205 3200 ELSE 3206 IF ( idsint == iup _u) THEN3201 IF ( idsint == iup ) THEN 3207 3202 DO m = 1, surf_usm_h%ns 3208 3203 i = surf_usm_h%i(m) … … 3226 3221 !-- Soil water content for iwl layer of walls and land 3227 3222 IF ( av == 0 ) THEN 3228 IF ( idsint == iup _u) THEN3223 IF ( idsint == iup ) THEN 3229 3224 DO m = 1, surf_usm_h%ns 3230 3225 i = surf_usm_h%i(m) … … 3237 3232 ENDIF 3238 3233 ELSE 3239 IF ( idsint == iup _u) THEN3234 IF ( idsint == iup ) THEN 3240 3235 DO m = 1, surf_usm_h%ns 3241 3236 i = surf_usm_h%i(m)
Note: See TracChangeset
for help on using the changeset viewer.
