Changeset 3814 for palm/trunk/SOURCE
- Timestamp:
- Mar 26, 2019 8:40:31 AM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/data_output_3d.f90
r3766 r3814 298 298 299 299 USE radiation_model_mod, & 300 ONLY: nz ub, nzut300 ONLY: nz_urban_b, nz_urban_t 301 301 302 302 USE turbulence_closure_mod, & … … 723 723 ! 724 724 !-- For urban model quantities, it is required to re-allocate local_pf 725 nzb_do = nz ub726 nzt_do = nz ut725 nzb_do = nz_urban_b 726 nzt_do = nz_urban_t 727 727 728 728 DEALLOCATE ( local_pf ) -
palm/trunk/SOURCE/radiation_model_mod.f90
r3771 r3814 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Rename exported nzu, nzp and related variables due to name conflict 31 ! 32 ! 3771 2019-02-28 12:19:33Z raasch 30 33 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary 31 34 ! pointers to avoid compiler warnings about outlived pointer targets, … … 584 587 585 588 USE basic_constants_and_equations_mod, & 586 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, 589 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, & 587 590 barometric_formula 588 591 … … 708 711 /) 709 712 710 REAL(wp), PARAMETER :: sigma_sb = 5.67037321E-8_wp !< Stefan-Boltzmann constant 711 712 INTEGER(iwp) :: albedo_type = 9999999, & !< Albedo surface type 713 dots_rad = 0 !< starting index for timeseries output 713 INTEGER(iwp) :: albedo_type = 9999999_iwp, & !< Albedo surface type 714 dots_rad = 0_iwp !< starting index for timeseries output 714 715 715 716 LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed … … 912 913 ! 913 914 !-- Parameters of urban and land surface models 914 INTEGER(iwp) :: nzu !< number of layers of urban surface (will be calculated) 915 INTEGER(iwp) :: nzp !< number of layers of plant canopy (will be calculated) 916 INTEGER(iwp) :: nzub,nzut !< bottom and top layer of urban surface (will be calculated) 917 INTEGER(iwp) :: nzpt !< top layer of plant canopy (will be calculated) 915 INTEGER(iwp) :: nz_urban !< number of layers of urban surface (will be calculated) 916 INTEGER(iwp) :: nz_plant !< number of layers of plant canopy (will be calculated) 917 INTEGER(iwp) :: nz_urban_b !< bottom layer of urban surface (will be calculated) 918 INTEGER(iwp) :: nz_urban_t !< top layer of urban surface (will be calculated) 919 INTEGER(iwp) :: nz_plant_t !< top layer of plant canopy (will be calculated) 918 920 !-- parameters of urban and land surface models 919 921 INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers above top of of topography … … 1104 1106 INTEGER(iwp) :: nmrtfa !< dimmension of array allocated for storage of mrt 1105 1107 INTEGER(iwp) :: msvf, mcsf, mmrtf!< mod for swapping the growing array 1106 INTEGER(iwp), PARAMETER :: gasize = 100000 !< initial size of growing arrays1108 INTEGER(iwp), PARAMETER :: gasize = 100000_iwp !< initial size of growing arrays 1107 1109 REAL(wp), PARAMETER :: grow_factor = 1.4_wp !< growth factor of growing arrays 1108 1110 INTEGER(iwp) :: nsvfl !< number of svf for local processor … … 1298 1300 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 1299 1301 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 1300 nsurf_type, nz ub, nzut, nzu, pch, nsurf,&1302 nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf, & 1301 1303 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 1302 1304 radiation_interactions, startwall, startland, endland, endwall, & … … 2843 2845 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp 2844 2846 2845 pt1_l = SUM( pt(nz ut,nys:nyn,nxl:nxr) )2846 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz ut,nys:nyn,nxl:nxr) )2847 pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) ) 2848 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) ) 2847 2849 2848 2850 #if defined( __parallel ) … … 2866 2868 #endif 2867 2869 2868 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz ut) * ql12870 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1 2869 2871 ! 2870 2872 !-- Finally, divide by number of grid points … … 2905 2907 IF ( average_radiation ) THEN 2906 2908 2907 k = nz ut2909 k = nz_urban_t 2908 2910 2909 2911 surf%rad_sw_in = solar_constant * sky_trans * zenith(0) … … 3023 3025 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp 3024 3026 3025 pt1_l = SUM( pt(nz ut,nys:nyn,nxl:nxr) )3026 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz ut,nys:nyn,nxl:nxr) )3027 pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) ) 3028 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) ) 3027 3029 3028 3030 #if defined( __parallel ) … … 3044 3046 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l 3045 3047 #endif 3046 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz ut+1) * ql13048 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1 3047 3049 ! 3048 3050 !-- Finally, divide by number of grid points … … 3086 3088 surf%rad_net = net_radiation 3087 3089 3088 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(nz ut+1))**43090 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(nz_urban_t+1))**4 3089 3091 3090 3092 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 4979 4981 REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx) 4980 4982 REAL(wp), DIMENSION(0:nsurf_type) :: costheta !< direct irradiance factor of solar angle 4981 REAL(wp), DIMENSION(nz ub:nzut) :: pchf_prep !< precalculated factor for canopy temperature tendency4983 REAL(wp), DIMENSION(nz_urban_b:nz_urban_t) :: pchf_prep !< precalculated factor for canopy temperature tendency 4982 4984 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4983 4985 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff … … 5002 5004 5003 5005 IF ( plant_canopy ) THEN 5004 pchf_prep(:) = r_d * exner(nz ub:nzut) &5005 / (c_p * hyp(nz ub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)5006 pchf_prep(:) = r_d * exner(nz_urban_b:nz_urban_t) & 5007 / (c_p * hyp(nz_urban_b:nz_urban_t) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T) 5006 5008 ENDIF 5007 5009 … … 5923 5925 ENDDO 5924 5926 ! 5925 !-- Find nz ub, nzut, nzuvia wall_flag_0 array (nzb_s_inner will be5927 !-- Find nz_urban_b, nz_urban_t, nz_urban via wall_flag_0 array (nzb_s_inner will be 5926 5928 !-- removed later). The following contruct finds the lowest / largest index 5927 5929 !-- for any upward-facing wall (see bit 12). … … 5976 5978 5977 5979 #if defined( __parallel ) 5978 CALL MPI_AllReduce(nzubl, nz ub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )5980 CALL MPI_AllReduce(nzubl, nz_urban_b, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr ) 5979 5981 IF ( ierr /= 0 ) THEN 5980 WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nz ub5982 WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nz_urban_b 5981 5983 FLUSH(9) 5982 5984 ENDIF 5983 CALL MPI_AllReduce(nzutl, nz ut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )5985 CALL MPI_AllReduce(nzutl, nz_urban_t, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5984 5986 IF ( ierr /= 0 ) THEN 5985 WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nz ut5987 WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nz_urban_t 5986 5988 FLUSH(9) 5987 5989 ENDIF 5988 CALL MPI_AllReduce(nzptl, nz pt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )5990 CALL MPI_AllReduce(nzptl, nz_plant_t, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5989 5991 IF ( ierr /= 0 ) THEN 5990 WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nz pt5992 WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nz_plant_t 5991 5993 FLUSH(9) 5992 5994 ENDIF 5993 5995 #else 5994 nz ub = nzubl5995 nz ut = nzutl5996 nz pt = nzptl5996 nz_urban_b = nzubl 5997 nz_urban_t = nzutl 5998 nz_plant_t = nzptl 5997 5999 #endif 5998 6000 ! … … 6002 6004 !-- are active. ABS (...) is required because the default value of 6003 6005 !-- dz_stretch_level_start is -9999999.9_wp (negative). 6004 IF ( ABS( dz_stretch_level_start(1) ) <= zw(nz ut) ) THEN6006 IF ( ABS( dz_stretch_level_start(1) ) <= zw(nz_urban_t) ) THEN 6005 6007 WRITE( message_string, * ) 'The lowest level where vertical ', & 6006 6008 'stretching is applied have to be ', & 6007 'greater than ', zw(nz ut)6009 'greater than ', zw(nz_urban_t) 6008 6010 CALL message( 'radiation_interaction_init', 'PA0496', 1, 2, 0, 6, 0 ) 6009 6011 ENDIF 6010 6012 ! 6011 6013 !-- global number of urban and plant layers 6012 nz u = nzut - nzub + 16013 nz p = nzpt - nzub + 16014 nz_urban = nz_urban_t - nz_urban_b + 1 6015 nz_plant = nz_plant_t - nz_urban_b + 1 6014 6016 ! 6015 6017 !-- check max_raytracing_dist relative to urban surface layer height 6016 mrl = 2.0_wp * nz u* dz(1)6018 mrl = 2.0_wp * nz_urban * dz(1) 6017 6019 !-- set max_raytracing_dist to double the urban surface layer height, if not set 6018 6020 IF ( max_raytracing_dist == -999.0_wp ) THEN … … 6063 6065 IF ( npcbl > 0 ) THEN 6064 6066 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 6065 ALLOCATE( gridpcbl(nz ub:nzpt,nys:nyn,nxl:nxr) )6067 ALLOCATE( gridpcbl(nz_urban_b:nz_plant_t,nys:nyn,nxl:nxr) ) 6066 6068 pcbl = -1 6067 6069 gridpcbl(:,:,:) = 0 … … 6125 6127 ENDIF 6126 6128 6127 CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz u*nny*nnx, &6129 CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz_urban*nny*nnx, & 6128 6130 kind=MPI_ADDRESS_KIND), STORAGE_SIZE(1_iwp)/8, & 6129 6131 minfo, comm2d, gridsurf_rma_p, win_gridsurf, ierr) 6130 6132 IF ( ierr /= 0 ) THEN 6131 6133 WRITE(9,*) 'Error MPI_Win_allocate1:', ierr, & 6132 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz u*nny*nnx,kind=MPI_ADDRESS_KIND), &6134 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nz_urban*nny*nnx,kind=MPI_ADDRESS_KIND), & 6133 6135 STORAGE_SIZE(1_iwp)/8, win_gridsurf 6134 6136 FLUSH(9) … … 6147 6149 !-- pointer and then redirecting a multidimensional pointer to it works 6148 6150 !-- fine. 6149 CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nz u*nny*nnx /))6150 gridsurf(0:nsurf_type_u-1, nz ub:nzut, nys:nyn, nxl:nxr) => &6151 gridsurf_rma(1:nsurf_type_u*nz u*nny*nnx)6151 CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nz_urban*nny*nnx /)) 6152 gridsurf(0:nsurf_type_u-1, nz_urban_b:nz_urban_t, nys:nyn, nxl:nxr) => & 6153 gridsurf_rma(1:nsurf_type_u*nz_urban*nny*nnx) 6152 6154 #else 6153 ALLOCATE(gridsurf(0:nsurf_type_u-1,nz ub:nzut,nys:nyn,nxl:nxr) )6155 ALLOCATE(gridsurf(0:nsurf_type_u-1,nz_urban_b:nz_urban_t,nys:nyn,nxl:nxr) ) 6154 6156 #endif 6155 6157 gridsurf(:,:,:,:) = -999 … … 6511 6513 IF ( plant_canopy ) THEN 6512 6514 ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) ) 6513 maxboxesg = nx + ny + nz p+ 16515 maxboxesg = nx + ny + nz_plant + 1 6514 6516 max_track_len = nx + ny + 1 6515 6517 !-- temporary arrays storing values for csf calculation during raytracing … … 6564 6566 !-- optimization of memory should be done 6565 6567 !-- Argument X of function STORAGE_SIZE(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now 6566 size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nz p6568 size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nz_plant 6567 6569 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, & 6568 6570 lad_s_rma_p, win_lad, ierr) … … 6572 6574 FLUSH(9) 6573 6575 ENDIF 6574 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nz p*nny*nnx /))6575 sub_lad(nz ub:nzpt, nys:nyn, nxl:nxr) => lad_s_rma(1:nzp*nny*nnx)6576 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nz_plant*nny*nnx /)) 6577 sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr) => lad_s_rma(1:nz_plant*nny*nnx) 6576 6578 ELSE 6577 ALLOCATE(sub_lad(nz ub:nzpt, nys:nyn, nxl:nxr))6579 ALLOCATE(sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr)) 6578 6580 ENDIF 6579 6581 #else 6580 6582 plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) ) 6581 ALLOCATE(sub_lad(nz ub:nzpt, nys:nyn, nxl:nxr))6583 ALLOCATE(sub_lad(nz_urban_b:nz_plant_t, nys:nyn, nxl:nxr)) 6582 6584 #endif 6583 6585 plantt_max = MAXVAL(plantt) 6584 ALLOCATE( rt2_track(2, max_track_len), rt2_track_lad(nz ub:plantt_max, max_track_len), &6585 rt2_track_dist(0:max_track_len), rt2_dist(plantt_max-nz ub+2) )6586 ALLOCATE( rt2_track(2, max_track_len), rt2_track_lad(nz_urban_b:plantt_max, max_track_len), & 6587 rt2_track_dist(0:max_track_len), rt2_dist(plantt_max-nz_urban_b+2) ) 6586 6588 6587 6589 sub_lad(:,:,:) = 0._wp … … 6590 6592 k = get_topography_top_index_ji( j, i, 's' ) 6591 6593 6592 sub_lad(k:nz pt, j, i) = lad_s(0:nzpt-k, j, i)6594 sub_lad(k:nz_plant_t, j, i) = lad_s(0:nz_plant_t-k, j, i) 6593 6595 ENDDO 6594 6596 ENDDO … … 6608 6610 6609 6611 ELSE 6610 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz p-1) )6611 CALL MPI_AllGather( sub_lad, nnx*nny*nz p, MPI_REAL, &6612 sub_lad_g, nnx*nny*nz p, MPI_REAL, comm2d, ierr )6612 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz_plant-1) ) 6613 CALL MPI_AllGather( sub_lad, nnx*nny*nz_plant, MPI_REAL, & 6614 sub_lad_g, nnx*nny*nz_plant, MPI_REAL, comm2d, ierr ) 6613 6615 IF ( ierr /= 0 ) THEN 6614 6616 WRITE(9,*) 'Error MPI_AllGather3:', ierr, SIZE(sub_lad), & 6615 nnx*nny*nz p, SIZE(sub_lad_g), nnx*nny*nzp6617 nnx*nny*nz_plant, SIZE(sub_lad_g), nnx*nny*nz_plant 6616 6618 FLUSH(9) 6617 6619 ENDIF … … 7551 7553 #if defined( __parallel ) 7552 7554 lad_ip(ncsb) = ip 7553 lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz p) + (box(2)-py*nny)*nzp + box(1)-nzub7555 lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz_plant) + (box(2)-py*nny)*nz_plant + box(1)-nz_urban_b 7554 7556 #endif 7555 7557 ENDIF … … 7595 7597 lad_s_target = lad_s_ray(i) 7596 7598 ELSE 7597 lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nz p+ lad_disp(i))7599 lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nz_plant + lad_disp(i)) 7598 7600 ENDIF 7599 7601 #else … … 7736 7738 rt2_track_dist(0) = 0._wp 7737 7739 rt2_track_lad(:,:) = 0._wp 7738 nly = plantt_max - nz ub + 17740 nly = plantt_max - nz_urban_b + 1 7739 7741 ENDIF 7740 7742 … … 7871 7873 IF ( plantt(ig) < lowest_lad ) CYCLE 7872 7874 7873 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz p) + (rt2_track(1,i)-py*nny)*nzp + lowest_lad-nzub7875 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant + lowest_lad-nz_urban_b 7874 7876 wcount = plantt(ig)-lowest_lad+1 7875 7877 ! TODO send request ASAP - even during raytracing … … 7897 7899 py = rt2_track(1,i)/nny 7898 7900 ip = px*pdims(2)+py 7899 ig = ip*nnx*nny*nz p + (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp7900 rt2_track_lad(nz ub:plantt_max, i) = sub_lad_g(ig:ig+nly-1)7901 ig = ip*nnx*nny*nz_plant + (rt2_track(2,i)-px*nnx)*(nny*nz_plant) + (rt2_track(1,i)-py*nny)*nz_plant 7902 rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad_g(ig:ig+nly-1) 7901 7903 ENDDO 7902 7904 ENDIF 7903 7905 #else 7904 7906 DO i = 1, ntrack 7905 rt2_track_lad(nz ub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max)7907 rt2_track_lad(nz_urban_b:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nz_urban_b:plantt_max) 7906 7908 ENDDO 7907 7909 #endif … … 7942 7944 !-- Assert that we have space allocated for CSFs 7943 7945 !-- 7944 maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nz ub, &7945 nz ut - CEILING(origin(1)-.5_wp))) * nrays7946 maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nz_urban_b, & 7947 nz_urban_t - CEILING(origin(1)-.5_wp))) * nrays 7946 7948 IF ( ncsfl + maxboxes > ncsfla ) THEN 7947 7949 !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) … … 7955 7957 !-- Calculate transparencies and store new CSFs 7956 7958 !-- 7957 zbottom = REAL(nz ub, wp) - .5_wp7959 zbottom = REAL(nz_urban_b, wp) - .5_wp 7958 7960 ztop = REAL(plantt_max, wp) + .5_wp 7959 7961 … … 8116 8118 py = y / nny 8117 8119 iproc = px * pdims(2) + py 8118 target_displ = ((x-px*nnx) * nny + y - py*nny ) * nz u* nsurf_type_u +&8119 ( z-nz ub ) * nsurf_type_u + d8120 target_displ = ((x-px*nnx) * nny + y - py*nny ) * nz_urban * nsurf_type_u +& 8121 ( z-nz_urban_b ) * nsurf_type_u + d 8120 8122 ! 8121 8123 !-- Send MPI_Get request to obtain index target_surfl(i) -
palm/trunk/SOURCE/urban_surface_mod.f90
r3802 r3814 480 480 iz, iy, ix, nsurf, idsvf, ndsvf, & 481 481 idcsf, ndcsf, kdcsf, pct, & 482 nz ub, nzut, unscheduled_radiation_calls482 nz_urban_b, nz_urban_t, unscheduled_radiation_calls 483 483 484 484 USE statistics, & … … 8585 8585 DO i = nxl, nxr 8586 8586 DO j = nys, nyn 8587 DO k = nz ub, min(nzut,naheatlayers)8587 DO k = nz_urban_b, min(nz_urban_t,naheatlayers) 8588 8588 IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN 8589 8589 !
Note: See TracChangeset
for help on using the changeset viewer.