Changeset 3449 for palm/trunk/SOURCE
- Timestamp:
- Oct 29, 2018 7:36:56 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 5 edited
-
basic_constants_and_equations_mod.f90 (modified) (6 diffs)
-
chemistry_model_mod.f90 (modified) (13 diffs)
-
plant_canopy_model_mod.f90 (modified) (21 diffs)
-
radiation_model_mod.f90 (modified) (46 diffs)
-
urban_surface_mod.f90 (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/basic_constants_and_equations_mod.f90
r3361 r3449 25 25 ! ----------------- 26 26 ! $Id$ 27 ! +degc_to_k 28 ! 29 ! 3361 2018-10-16 20:39:37Z knoop 27 30 ! New module (introduced with modularization of bulk cloud physics model) 28 31 ! … … 43 46 IMPLICIT NONE 44 47 45 REAL(wp), PARAMETER :: pi = 3.141592654_wp !< PI46 47 48 REAL(wp), PARAMETER :: c_p = 1005.0_wp !< heat capacity of dry air (J kg-1 K-1) 49 REAL(wp), PARAMETER :: degc_to_k = 273.15_wp !< temperature (in K) of 0 deg C (K) 48 50 REAL(wp), PARAMETER :: g = 9.81_wp !< gravitational acceleration (m s-2) 49 51 REAL(wp), PARAMETER :: kappa = 0.4_wp !< von Karman constant … … 55 57 REAL(wp), PARAMETER :: molecular_weight_of_nh4no3 = 0.08004_wp !< mol. m. ammonium sulfate (kg mol-1) 56 58 REAL(wp), PARAMETER :: molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1) 59 REAL(wp), PARAMETER :: pi = 3.141592654_wp !< PI 57 60 REAL(wp), PARAMETER :: rho_l = 1.0E3_wp !< density of water (kg m-3) 58 61 REAL(wp), PARAMETER :: rho_nacl = 2165.0_wp !< density of NaCl (kg m-3) … … 61 64 REAL(wp), PARAMETER :: r_d = 287.0_wp !< sp. gas const. dry air (J kg-1 K-1) 62 65 REAL(wp), PARAMETER :: r_v = 461.51_wp !< sp. gas const. water vapor (J kg-1 K-1) 63 REAL(wp), PARAMETER :: sigma_sb = 5.67 E-08_wp!< Stefan-Boltzmann constant66 REAL(wp), PARAMETER :: sigma_sb = 5.67037E-08_wp !< Stefan-Boltzmann constant 64 67 REAL(wp), PARAMETER :: solar_constant = 1368.0_wp !< solar constant at top of atmosphere 65 68 REAL(wp), PARAMETER :: vanthoff_nacl = 2.0_wp !< van't Hoff factor for NaCl … … 143 146 ! 144 147 !-- Saturation vapor pressure for a specific temperature: 145 magnus_0d = 611.2_wp * EXP( 17.62_wp * ( t - 273.15_wp) / &148 magnus_0d = 611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) / & 146 149 ( t - 29.65_wp ) ) 147 150 … … 163 166 ! 164 167 !-- Saturation vapor pressure for a specific temperature: 165 magnus_1d = 611.2_wp * EXP( 17.62_wp * ( t - 273.15_wp) / &168 magnus_1d = 611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) / & 166 169 ( t - 29.65_wp ) ) 167 170 -
palm/trunk/SOURCE/chemistry_model_mod.f90
r3435 r3449 27 27 ! ----------------- 28 28 ! $Id$ 29 ! additional output - merged from branch resler 30 ! 31 ! 3438 2018-10-28 19:31:42Z pavelkrc 29 32 ! Add terrain-following masked output 30 33 ! … … 262 265 !< neumann = zero gradient 263 266 264 REAL(kind=wp), PUBLIC :: cs_time_step 267 REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp 265 268 CHARACTER(10), PUBLIC :: photolysis_scheme 266 269 ! 'constant', … … 410 413 IF ( mode == 'allocate' ) THEN 411 414 DO lsp = 1, nspec 412 IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 415 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 416 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 413 417 chem_species(lsp)%conc_av = 0.0_wp 414 415 418 ENDIF 416 419 ENDDO … … 418 421 ELSEIF ( mode == 'sum' ) THEN 419 422 420 421 423 DO lsp = 1, nspec 422 IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 424 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 425 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 423 426 DO i = nxlg, nxrg 424 427 DO j = nysg, nyng … … 462 465 ELSEIF ( mode == 'average' ) THEN 463 466 DO lsp = 1, nspec 464 IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 467 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 468 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 465 469 DO i = nxlg, nxrg 466 470 DO j = nysg, nyng … … 811 815 unit = 'illegal' 812 816 813 spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_'. 817 spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_' or 'em_'. 818 819 IF ( TRIM(var(1:3)) == 'em_' ) THEN 820 DO lsp=1,nspec 821 IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN 822 unit = 'mol m-2 s-1' 823 ENDIF 824 !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code 825 !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): 826 !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5) 827 IF (spec_name(1:2) == 'PM') THEN 828 unit = 'kg m-2 s-1' 829 ENDIF 830 ENDDO 831 832 ELSE 814 833 815 834 DO lsp=1,nspec … … 821 840 ! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5) 822 841 IF (spec_name(1:2) == 'PM') THEN 823 unit = 'kg m-3'842 unit = 'kg m-3' 824 843 ENDIF 825 844 ENDDO 826 845 827 DO lsp=1,nphot 846 DO lsp=1,nphot 828 847 IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN 829 848 unit = 'sec-1' 830 849 ENDIF 831 850 ENDDO 832 833 834 RETURN 851 ENDIF 852 853 854 RETURN 835 855 END SUBROUTINE chem_check_data_output 836 856 ! … … 933 953 CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 ) 934 954 ENDIF 955 956 !--------------------- 957 !-- chem_check_parameters is called before the array chem_species is allocated! 958 !-- temporary switch of this part of the check 959 RETURN 960 !--------------------- 935 961 936 962 !-- check for initial chem species input … … 1076 1102 USE kinds 1077 1103 1104 USE surface_mod 1105 1106 !USE chem_modules, ONLY: nvar 1078 1107 1079 1108 IMPLICIT NONE … … 1091 1120 1092 1121 !-- local variables 1093 1094 INTEGER :: i, j, k, lsp1095 1122 CHARACTER(len=16) :: spec_name 1123 INTEGER(iwp) :: i, j, k 1124 INTEGER(iwp) :: m, l !< running indices for surfaces 1125 INTEGER(iwp) :: lsp !< running index for chem spcs 1096 1126 1097 1127 1098 1128 found = .FALSE. 1129 IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN 1130 RETURN 1131 ENDIF 1099 1132 1100 1133 spec_name = TRIM(variable(4:)) 1101 1134 1102 DO lsp=1,nspec 1103 IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN 1135 IF ( variable(1:3) == 'em_' ) THEN 1136 1137 local_pf = 0.0_wp 1138 1139 DO lsp = 1, nvar !!! cssws - nvar species, chem_species - nspec species !!! 1140 IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) ) THEN 1141 ! no average for now 1142 DO m = 1, surf_usm_h%ns 1143 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = & 1144 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m) 1145 ENDDO 1146 DO m = 1, surf_lsm_h%ns 1147 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = & 1148 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m) 1149 ENDDO 1150 DO l = 0, 3 1151 DO m = 1, surf_usm_v(l)%ns 1152 local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = & 1153 local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m) 1154 ENDDO 1155 DO m = 1, surf_lsm_v(l)%ns 1156 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = & 1157 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m) 1158 ENDDO 1159 ENDDO 1160 found = .TRUE. 1161 ENDIF 1162 ENDDO 1163 ELSE 1164 DO lsp=1,nspec 1165 IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN 1104 1166 ! 1105 1167 !-- todo: remove or replace by "CALL message" mechanism (kanani) 1106 1168 ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // & 1107 1169 ! TRIM(chem_species(lsp)%name) 1108 1109 IF (av == 0) THEN 1110 DO i = nxl, nxr 1111 DO j = nys, nyn 1112 DO k = nzb_do, nzt_do 1113 local_pf(i,j,k) = MERGE( & 1114 chem_species(lsp)%conc(k,j,i), & 1115 REAL( fill_value, KIND = wp ), & 1116 BTEST( wall_flags_0(k,j,i), 0 ) ) 1117 ENDDO 1118 ENDDO 1119 ENDDO 1120 1121 ELSE 1122 DO i = nxl, nxr 1123 DO j = nys, nyn 1124 DO k = nzb_do, nzt_do 1125 local_pf(i,j,k) = MERGE( & 1126 chem_species(lsp)%conc_av(k,j,i),& 1127 REAL( fill_value, KIND = wp ), & 1128 BTEST( wall_flags_0(k,j,i), 0 ) ) 1129 ENDDO 1130 ENDDO 1131 ENDDO 1132 ENDIF 1133 1134 found = .TRUE. 1135 ENDIF 1136 ENDDO 1170 IF (av == 0) THEN 1171 DO i = nxl, nxr 1172 DO j = nys, nyn 1173 DO k = nzb_do, nzt_do 1174 local_pf(i,j,k) = MERGE( & 1175 chem_species(lsp)%conc(k,j,i), & 1176 REAL( fill_value, KIND = wp ), & 1177 BTEST( wall_flags_0(k,j,i), 0 ) ) 1178 ENDDO 1179 ENDDO 1180 ENDDO 1181 1182 ELSE 1183 DO i = nxl, nxr 1184 DO j = nys, nyn 1185 DO k = nzb_do, nzt_do 1186 local_pf(i,j,k) = MERGE( & 1187 chem_species(lsp)%conc_av(k,j,i),& 1188 REAL( fill_value, KIND = wp ), & 1189 BTEST( wall_flags_0(k,j,i), 0 ) ) 1190 ENDDO 1191 ENDDO 1192 ENDDO 1193 ENDIF 1194 found = .TRUE. 1195 ENDIF 1196 ENDDO 1197 ENDIF 1137 1198 1138 1199 RETURN 1200 1139 1201 END SUBROUTINE chem_data_output_3d 1140 1202 ! … … 1177 1239 1178 1240 spec_name = TRIM( variable(4:) ) 1179 !av == 01241 !av == 0 1180 1242 1181 1243 DO lsp=1,nspec … … 1292 1354 found = .TRUE. 1293 1355 1294 IF ( var(1:3) == 'kc_') THEN !< always the same grid for chemistry variables1356 IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' ) THEN !< always the same grid for chemistry variables 1295 1357 grid_x = 'x' 1296 1358 grid_y = 'y' … … 1873 1935 !-- Enable chemistry model 1874 1936 air_chemistry = .TRUE. 1875 GOTO 201876 1877 10 BACKSPACE( 11 )1878 READ( 11 , '(A)') line1879 CALL parin_fail_message( 'chemistry_parameters', line )1880 1881 20 CONTINUE1937 GOTO 20 1938 1939 10 BACKSPACE( 11 ) 1940 READ( 11 , '(A)') line 1941 CALL parin_fail_message( 'chemistry_parameters', line ) 1942 1943 20 CONTINUE 1882 1944 1883 1945 ! -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r3337 r3449 16 16 ! 17 17 ! Copyright 1997-2018 Leibniz Universitaet Hannover 18 ! Copyright 2018 Institute of Computer Science of the 19 ! Czech Academy of Sciences, Prague 18 20 !------------------------------------------------------------------------------! 19 21 ! … … 25 27 ! ----------------- 26 28 ! $Id$ 29 ! Add calculation of transpiration for resolved plant canopy (trees, shrubs) 30 ! (V. Fuka, MFF UK Prague, J.Resler, ICS AS, Prague) 31 ! 32 ! Fix reading plant canopy over buildings 33 ! 34 ! 3337 2018-10-12 15:17:09Z kanani 27 35 ! Fix reading plant canopy over buildings 28 36 ! … … 214 222 215 223 USE arrays_3d, & 216 ONLY: dzu, dzw, e, q, s, tend, u, v, w, zu, zw 224 ONLY: dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw 225 226 USE basic_constants_and_equations_mod, & 227 ONLY: c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv 228 229 USE control_parameters, & 230 ONLY: humidity 217 231 218 232 USE indices, & … … 222 236 USE kinds 223 237 238 USE pegrid 239 224 240 USE surface_mod, & 225 241 ONLY: get_topography_top_index_ji … … 229 245 230 246 231 CHARACTER (LEN=30) :: canopy_mode = 'block' !< canopy coverage 232 233 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 234 INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 235 236 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top 237 238 LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func. 247 CHARACTER (LEN=30) :: canopy_mode = 'block' !< canopy coverage 248 LOGICAL :: plant_canopy_transpiration = .FALSE. !< flag to switch calculation of transpiration and corresponding latent heat 249 !< for resolved plant canopy inside radiation model 250 !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod) 251 252 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 253 INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 254 255 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top 256 257 LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func. 239 258 240 259 REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation … … 261 280 REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad 262 281 263 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.264 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid265 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_heating_rate !< plant canopy heating rate282 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc. 283 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid 284 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_heating_rate !< plant canopy heating rate 266 285 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_transpiration_rate !< plant canopy transpiration rate 286 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_latent_rate !< plant canopy latent heating rate 267 287 268 288 SAVE … … 273 293 ! 274 294 !-- Public functions 275 PUBLIC pcm_check_data_output, pcm_check_parameters, pcm_data_output_3d, & 276 pcm_define_netcdf_grid, pcm_header, pcm_init, pcm_parin, pcm_tendency 295 PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output, & 296 pcm_check_parameters, pcm_data_output_3d, pcm_define_netcdf_grid, & 297 pcm_header, pcm_init, pcm_parin, pcm_tendency 277 298 278 299 ! 279 300 !-- Public variables and constants 280 PUBLIC pc_heating_rate, pc_transpiration_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, & 281 pch_index 282 301 PUBLIC pc_heating_rate, pc_transpiration_rate, pc_latent_rate, & 302 canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index, & 303 plant_canopy_transpiration 304 305 INTERFACE pcm_calc_transpiration_rate 306 MODULE PROCEDURE pcm_calc_transpiration_rate 307 END INTERFACE pcm_calc_transpiration_rate 283 308 284 309 INTERFACE pcm_check_data_output … … 323 348 324 349 350 351 !------------------------------------------------------------------------------! 352 ! Description: 353 ! ------------ 354 !> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart 355 !> with parametrizations described in Daudet et al. (1999; Agricult. and Forest 356 !> Meteorol. 97) and Ngao, Adam and Saudreau (2017; Agricult. and Forest Meteorol 357 !> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric. 358 !> and Forest. Meteorol. 43) instead, because they are valid for broader intervals 359 !> of values. Funcion f4 used in form present in van Wijk et al. (1998; 360 !> Tree Physiology 20). 361 !> 362 !> This subroutine is called from subroutine radiation_interaction 363 !> after the calculation of radiation in plant canopy boxes. 364 !> (arrays pcbinsw and pcbinlw). 365 !> 366 !------------------------------------------------------------------------------! 367 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh) 368 369 USE control_parameters, & 370 ONLY: dz 371 372 USE grid_variables, & 373 ONLY: dx, dy 374 375 IMPLICIT NONE 376 !-- input parameters 377 INTEGER(iwp), INTENT(IN) :: i, j, k, kk !< indices of the pc gridbox 378 REAL(wp), INTENT(IN) :: pcbsw !< sw radiation in gridbox (W) 379 REAL(wp), INTENT(IN) :: pcblw !< lw radiation in gridbox (W) 380 REAL(wp), INTENT(OUT) :: pcbtr !< transpiration rate dq/dt (kg/kg/s) 381 REAL(wp), INTENT(OUT) :: pcblh !< latent heat from transpiration dT/dt (K/s) 382 383 !-- variables and parameters for calculation of transpiration rate 384 REAL(wp) :: sat_press, sat_press_d, temp, v_lad 385 REAL(wp) :: d_fact, g_b, g_s, wind_speed, evapor_rate 386 REAL(wp) :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad 387 REAL(wp), PARAMETER :: gama_psychr = 66 !< psychrometric constant (Pa/K) 388 REAL(wp), PARAMETER :: g_s_max = 0.01 !< maximum stomatal conductivity (m/s) 389 REAL(wp), PARAMETER :: m_soil = 0.4_wp !< soil water content (needs to adjust or take from LSM) 390 REAL(wp), PARAMETER :: m_wilt = 0.01_wp !< wilting point soil water content (needs to adjust or take from LSM) 391 REAL(wp), PARAMETER :: m_sat = 0.51_wp !< saturation soil water content (needs to adjust or take from LSM) 392 REAL(wp), PARAMETER :: t2_min = 0.0_wp !< minimal temperature for calculation of f2 393 REAL(wp), PARAMETER :: t2_max = 40.0_wp !< maximal temperature for calculation of f2 394 395 396 !-- Temperature (deg C) 397 temp = pt(k,j,i) * exner(k) - degc_to_k 398 !-- Coefficient for conversion of radiation to grid to radiation to unit leaves surface 399 v_lad = 1.0_wp / ( MAX(lad_s(kk,j,i), 1.0e-10) * dx * dy * dz(1) ) 400 !-- Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1) 401 !-- There are updated formulas available, kept consistent with the rest of the parametrization 402 sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp)) 403 !-- Saturation pressure derivative (derivative of the above) 404 sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2 405 !-- Wind speed 406 wind_speed = SQRT( ((u(k,j,i) + u(k,j,i+1))/2.0_wp )**2 + & 407 ((v(k,j,i) + v(k,j,i+1))/2.0_wp )**2 + & 408 ((w(k,j,i) + w(k,j,i+1))/2.0_wp )**2 ) 409 !-- Aerodynamic conductivity (Daudet et al. (1999) eq. 14 410 g_b = 0.01_wp * wind_speed + 0.0071_wp 411 !-- Radiation flux per leaf surface unit 412 rad = pcbsw * v_lad 413 !-- First function for calculation of stomatal conductivity (radiation dependency) 414 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17 415 f1 = rad * (1000._wp+42.1_wp) / 1000._wp / (rad+42.1_wp) 416 !-- Second function for calculation of stomatal conductivity (temperature dependency) 417 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21 418 f2 = MAX(t2_min, (temp-t2_min) * (t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / & 419 ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) ) 420 !-- Water pressure deficit 421 !-- Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure 422 vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp ) 423 !-- Third function for calculation of stomatal conductivity (water pressure deficit dependency) 424 !-- Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988) 425 !-- The coefficients of the linear dependence should better correspond to broad-leaved trees 426 !-- than the coefficients from Stewart (1988) which correspond to conifer trees. 427 vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp) 428 f3 = -2e-4_wp * vpd + 1.154_wp 429 !-- Fourth function for calculation of stomatal conductivity (soil moisture dependency) 430 !-- Residual soil water content 431 !-- van Wijk et al. (1998; Tree Physiology 20) eq. 7 432 !-- TODO - over LSM surface might be calculated from LSM parameters 433 rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt ) 434 !-- van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988)) 435 f4 = MAX(0._wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp)) 436 !-- Stomatal conductivity 437 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12 438 !-- (notation according to Ngao, Adam and Saudreau (2017) and others) 439 g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0e-10_wp 440 !-- Decoupling factor 441 !-- Daudet et al. (1999) eq. 6 442 d_fact = (sat_press_d / gama_psychr + 2._wp ) / & 443 (sat_press_d / gama_psychr + 2._wp + 2 * g_b / g_s ) 444 !-- Equilibrium evaporation rate 445 !-- Daudet et al. (1999) eq. 4 446 e_eq = (pcbsw + pcblw) * v_lad * sat_press_d / & 447 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v 448 !-- Imposed evaporation rate 449 !-- Daudet et al. (1999) eq. 5 450 e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v 451 !-- Evaporation rate 452 !-- Daudet et al. (1999) eq. 3 453 !-- (evaporation rate is limited to non-negative values) 454 evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp) 455 !-- Conversion of evaporation rate to q tendency in gridbox 456 !-- dq/dt = E * LAD * V_g / (rho_air * V_g) 457 pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k) !-- = dq/dt 458 !-- latent heat from evaporation 459 pcblh = pcbtr * lv_d_cp !-- = - dT/dt 460 461 END SUBROUTINE pcm_calc_transpiration_rate 462 463 325 464 !------------------------------------------------------------------------------! 326 465 ! Description: … … 352 491 CASE ( 'pcm_transpirationrate' ) 353 492 unit = 'kg kg-1 s-1' 493 494 CASE ( 'pcm_latentrate' ) 495 unit = 'K s-1' 496 497 CASE ( 'pcm_bowenratio' ) 498 unit = 'K s-1' 354 499 355 500 CASE ( 'pcm_lad' ) … … 503 648 ENDDO 504 649 ENDIF 505 650 651 CASE ( 'pcm_latentrate' ) 652 IF ( av == 0 ) THEN 653 DO i = nxl, nxr 654 DO j = nys, nyn 655 IF ( pch_index_ji(j,i) /= 0 ) THEN 656 k_topo = get_topography_top_index_ji( j, i, 's' ) 657 DO k = k_topo, k_topo + pch_index_ji(j,i) 658 local_pf(i,j,k) = pc_latent_rate(k-k_topo,j,i) 659 ENDDO 660 ENDIF 661 ENDDO 662 ENDDO 663 ENDIF 664 665 CASE ( 'pcm_bowenratio' ) 666 IF ( av == 0 ) THEN 667 DO i = nxl, nxr 668 DO j = nys, nyn 669 IF ( pch_index_ji(j,i) /= 0 ) THEN 670 k_topo = get_topography_top_index_ji( j, i, 's' ) 671 DO k = k_topo, k_topo + pch_index_ji(j,i) 672 IF ( pc_latent_rate(k-k_topo,j,i) /= 0._wp ) THEN 673 local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) / & 674 pc_latent_rate(k-k_topo,j,i) 675 ENDIF 676 ENDDO 677 ENDIF 678 ENDDO 679 ENDDO 680 ENDIF 681 506 682 CASE ( 'pcm_lad' ) 507 683 IF ( av == 0 ) THEN … … 550 726 SELECT CASE ( TRIM( var ) ) 551 727 552 CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate' )728 CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate', 'pcm_bowenratio') 553 729 grid_x = 'x' 554 730 grid_y = 'y' … … 691 867 692 868 USE control_parameters, & 693 ONLY: humidity,message_string, ocean_mode, urban_surface869 ONLY: message_string, ocean_mode, urban_surface 694 870 695 871 USE netcdf_data_input_mod, & … … 705 881 INTEGER(iwp) :: k !< running index 706 882 INTEGER(iwp) :: m !< running index 883 INTEGER(iwp) :: pch_index_l 707 884 708 885 REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction … … 903 1080 CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp ) 904 1081 1082 !-- Calculate global pch_index value (index of top of plant canopy from ground) 1083 pch_index = MAXVAL(pch_index_ji) 1084 !-- Exchange pch_index from all processors 1085 #if defined( __parallel ) 1086 CALL MPI_Allreduce(MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr) 1087 IF ( ierr /= 0 ) THEN 1088 WRITE (9,*) 'Error MPI_Allreduce pch_index:', ierr 1089 FLUSH(9) 1090 ENDIF 1091 #endif 1092 1093 !-- Allocation of arrays pc_heating_rate, pc_transpiration_rate and pc_latent_rate 1094 ALLOCATE( pc_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1095 IF ( humidity ) THEN 1096 ALLOCATE( pc_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1097 pc_transpiration_rate = 0.0_wp 1098 ALLOCATE( pc_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1099 pc_latent_rate = 0.0_wp 1100 ENDIF 1101 905 1102 ! 906 1103 !-- Initialization of the canopy heat source distribution due to heating … … 915 1112 !-- 73â96). This approach has been applied e.g. by Shaw & Schumann (1992; 916 1113 !-- Bound.-Layer Meteorol. 61, 47â64). 917 !-- When using the urban surface model (USM), canopy heating (pc_heating_rate) 918 !-- by radiation is calculated in the USM. 919 IF ( cthf /= 0.0_wp .AND. .NOT. urban_surface ) THEN 920 921 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 922 pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1114 !-- When using the radiation_interactions, canopy heating (pc_heating_rate) 1115 !-- and plant canopy transpiration (pc_transpiration_rate, pc_latent_rate) 1116 !-- are calculated in the RTM after the calculation of radiation. 1117 !-- We cannot use variable radiation_interactions here to determine the situation 1118 !-- as it is assigned in init_3d_model after the call of pcm_init. 1119 IF ( cthf /= 0.0_wp ) THEN 1120 1121 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 923 1122 ! 924 1123 !-- Piecewise calculation of the cumulative leaf area index by vertical … … 996 1195 997 1196 ENDIF 998 !999 !-- Allocate transpiration rate1000 IF ( humidity ) &1001 ALLOCATE( pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )1002 1003 1197 1004 1198 … … 1028 1222 lai_beta, & 1029 1223 leaf_scalar_exch_coeff, & 1030 leaf_surface_conc, pch_index 1224 leaf_surface_conc, pch_index, & 1225 plant_canopy_transpiration 1031 1226 1032 1227 NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, & … … 1037 1232 lai_beta, & 1038 1233 leaf_scalar_exch_coeff, & 1039 leaf_surface_conc, pch_index 1234 leaf_surface_conc, pch_index, & 1235 plant_canopy_transpiration 1040 1236 1041 1237 line = ' ' … … 1434 1630 !-- potential temperature 1435 1631 CASE ( 4 ) 1436 DO i = nxl, nxr 1437 DO j = nys, nyn 1438 ! 1439 !-- Determine topography-top index on scalar-grid 1440 k_wall = get_topography_top_index_ji( j, i, 's' ) 1441 1442 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1443 1444 kk = k - k_wall !- lad arrays are defined flat 1445 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1632 IF ( humidity ) THEN 1633 DO i = nxl, nxr 1634 DO j = nys, nyn 1635 !-- Determine topography-top index on scalar-grid 1636 k_wall = get_topography_top_index_ji( j, i, 's' ) 1637 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1638 kk = k - k_wall !- lad arrays are defined flat 1639 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i) 1640 ENDDO 1446 1641 ENDDO 1447 1642 ENDDO 1448 ENDDO 1643 ELSE 1644 DO i = nxl, nxr 1645 DO j = nys, nyn 1646 !-- Determine topography-top index on scalar-grid 1647 k_wall = get_topography_top_index_ji( j, i, 's' ) 1648 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1649 kk = k - k_wall !- lad arrays are defined flat 1650 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1651 ENDDO 1652 ENDDO 1653 ENDDO 1654 ENDIF 1449 1655 1450 1656 ! … … 1460 1666 1461 1667 kk = k - k_wall !- lad arrays are defined flat 1462 pc_transpiration_rate(kk,j,i) = - lsec & 1463 * lad_s(kk,j,i) * & 1464 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1465 u(k,j,i+1) ) & 1466 )**2 + & 1467 ( 0.5_wp * ( v(k,j,i) + & 1468 v(k,j+1,i) ) & 1469 )**2 + & 1470 ( 0.5_wp * ( w(k-1,j,i) + & 1471 w(k,j,i) ) & 1472 )**2 & 1473 ) * & 1474 ( q(k,j,i) - lsc ) 1668 1669 IF ( .NOT. plant_canopy_transpiration ) THEN 1670 ! pc_transpiration_rate is calculated in radiation model 1671 ! in case of plant_canopy_transpiration = .T. 1672 ! to include also the dependecy to the radiation 1673 ! in the plant canopy box 1674 pc_transpiration_rate(kk,j,i) = - lsec & 1675 * lad_s(kk,j,i) * & 1676 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1677 u(k,j,i+1) ) & 1678 )**2 + & 1679 ( 0.5_wp * ( v(k,j,i) + & 1680 v(k,j+1,i) ) & 1681 )**2 + & 1682 ( 0.5_wp * ( w(k-1,j,i) + & 1683 w(k,j,i) ) & 1684 )**2 & 1685 ) * & 1686 ( q(k,j,i) - lsc ) 1687 ENDIF 1475 1688 1476 1689 tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i) … … 1778 1991 k_wall = get_topography_top_index_ji( j, i, 's' ) 1779 1992 1993 IF ( humidity ) THEN 1994 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1995 kk = k - k_wall !- lad arrays are defined flat 1996 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i) 1997 ENDDO 1998 ELSE 1999 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 2000 kk = k - k_wall !- lad arrays are defined flat 2001 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 2002 ENDDO 2003 ENDIF 2004 2005 ! 2006 !-- humidity 2007 CASE ( 5 ) 2008 ! 2009 !-- Determine topography-top index on scalar grid 2010 k_wall = get_topography_top_index_ji( j, i, 's' ) 2011 1780 2012 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1781 2013 kk = k - k_wall !- lad arrays are defined flat 1782 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1783 ENDDO 1784 1785 1786 ! 1787 !-- humidity 1788 CASE ( 5 ) 1789 ! 1790 !-- Determine topography-top index on scalar grid 1791 k_wall = get_topography_top_index_ji( j, i, 's' ) 1792 1793 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1794 kk = k - k_wall !- lad arrays are defined flat 1795 1796 pc_transpiration_rate(kk,j,i) = - lsec & 1797 * lad_s(kk,j,i) * & 1798 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1799 u(k,j,i+1) ) & 1800 )**2 + & 1801 ( 0.5_wp * ( v(k,j,i) + & 1802 v(k,j+1,i) ) & 1803 )**2 + & 1804 ( 0.5_wp * ( w(k-1,j,i) + & 1805 w(k,j,i) ) & 1806 )**2 & 1807 ) * & 1808 ( q(k,j,i) - lsc ) 2014 IF ( .NOT. plant_canopy_transpiration ) THEN 2015 ! pc_transpiration_rate is calculated in radiation model 2016 ! in case of plant_canopy_transpiration = .T. 2017 ! to include also the dependecy to the radiation 2018 ! in the plant canopy box 2019 pc_transpiration_rate(kk,j,i) = - lsec & 2020 * lad_s(kk,j,i) * & 2021 SQRT( ( 0.5_wp * ( u(k,j,i) + & 2022 u(k,j,i+1) ) & 2023 )**2 + & 2024 ( 0.5_wp * ( v(k,j,i) + & 2025 v(k,j+1,i) ) & 2026 )**2 + & 2027 ( 0.5_wp * ( w(k-1,j,i) + & 2028 w(k,j,i) ) & 2029 )**2 & 2030 ) * & 2031 ( q(k,j,i) - lsc ) 2032 ENDIF 1809 2033 1810 2034 tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i) -
palm/trunk/SOURCE/radiation_model_mod.f90
r3435 r3449 28 28 ! ----------------- 29 29 ! $Id$ 30 ! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague) 31 ! - Interaction of plant canopy with LW radiation 32 ! - Transpiration from resolved plant canopy dependent on radiation 33 ! called from RTM 34 ! 35 ! 36 ! 3435 2018-10-26 18:25:44Z gronemeier 30 37 ! - workaround: return unit=illegal in check_data_output for certain variables 31 38 ! when check called from init_masks … … 499 506 500 507 USE basic_constants_and_equations_mod, & 501 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, 508 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, & 502 509 barometric_formula 503 510 … … 544 551 545 552 USE plant_canopy_model_mod, & 546 ONLY: lad_s, pc_heating_rate, pc_transpiration_rate 553 ONLY: lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate, & 554 plant_canopy_transpiration, pcm_calc_transpiration_rate 547 555 548 556 USE pegrid … … 853 861 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices 854 862 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices 855 !< parameter but set in the code 863 REAL(wp), DIMENSION(0:nsurf_type) :: facearea !< area of single face in respective 864 !< direction (will be calc'd) 856 865 857 866 … … 887 896 !-- configuration parameters (they can be setup in PALM config) 888 897 LOGICAL :: raytrace_mpi_rma = .TRUE. !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing 889 LOGICAL :: read_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be read from a file at initialization 890 LOGICAL :: write_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be written out to a file 891 LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 898 LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 892 899 !< reflected radiation (as opposed to all mutually visible pairs) 900 LOGICAL :: plant_lw_interact = .TRUE. !< whether plant canopy interacts with LW radiation (in addition to SW) 893 901 INTEGER(iwp) :: mrt_nlevels = 0 !< number of vertical boxes above surface for which to calculate MRT 894 902 LOGICAL :: mrt_skip_roof = .TRUE. !< do not calculate MRT above roof surfaces 895 903 LOGICAL :: mrt_include_sw = .TRUE. !< should MRT calculation include SW radiation as well? 896 INTEGER(iwp) :: nrefsteps = 0!< number of reflection steps to perform904 INTEGER(iwp) :: nrefsteps = 3 !< number of reflection steps to perform 897 905 REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha) 898 906 INTEGER(iwp), PARAMETER :: rad_version_len = 10 !< length of identification string of rad version 899 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 1.1' !< identification of version of binary svf and restart files907 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 3.0' !< identification of version of binary svf and restart files 900 908 INTEGER(iwp) :: raytrace_discrete_elevs = 40 !< number of discretization steps for elevation (nadir to zenith) 901 909 INTEGER(iwp) :: raytrace_discrete_azims = 80 !< number of discretization steps for azimuth (out of 360 degrees) … … 929 937 INTEGER(iwp) :: ity !< 930 938 INTEGER(iwp) :: itz !< 931 INTEGER(iwp) :: isurfs !< 932 REAL(wp) :: rsvf !< 939 INTEGER(iwp) :: isurfs !< Idx of source face / -1 for sky 940 REAL(wp) :: rcvf !< Canopy view factor for faces / 941 !< canopy sink factor for sky (-1) 933 942 END TYPE 934 943 … … 976 985 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection 977 986 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection 987 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlg !< global array of incoming lw radiation from plant canopy 978 988 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection 979 989 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection … … 2928 2938 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 2929 2939 2930 NAMELIST /radiation_par/ albedo, albedo_ type, albedo_lw_dir,&2931 albedo_ lw_dif, albedo_sw_dir, albedo_sw_dif,&2932 constant_albedo, dt_radiation, emissivity, &2933 lw_radiation, m rt_nlevels, mrt_skip_roof,&2934 m rt_include_sw, net_radiation,&2935 radiation_scheme, skip_time_do_radiation,&2936 sw_radiation, unscheduled_radiation_calls,&2937 max_raytracing_dist, min_irrf_value,&2938 nrefsteps, raytrace_mpi_rma,&2939 surface_reflections, svfnorm_report_thresh,&2940 radiation_interactions_on,&2941 rad_angular_discretization,&2942 raytrace_discrete_azims, &2943 raytrace_discrete_elevs 2940 NAMELIST /radiation_par/ albedo, albedo_lw_dif, albedo_lw_dir, & 2941 albedo_sw_dif, albedo_sw_dir, albedo_type, & 2942 constant_albedo, dt_radiation, emissivity, & 2943 lw_radiation, max_raytracing_dist, & 2944 min_irrf_value, mrt_include_sw, mrt_nlevels, & 2945 mrt_skip_roof, net_radiation, nrefsteps, & 2946 plant_lw_interact, rad_angular_discretization,& 2947 radiation_interactions_on, radiation_scheme, & 2948 raytrace_discrete_azims, & 2949 raytrace_discrete_elevs, raytrace_mpi_rma, & 2950 skip_time_do_radiation, surface_reflections, & 2951 svfnorm_report_thresh, sw_radiation, & 2952 unscheduled_radiation_calls 2953 2944 2954 2945 NAMELIST /radiation_parameters/ albedo, albedo_type, albedo_lw_dir, & 2946 albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, & 2947 constant_albedo, dt_radiation, emissivity, & 2948 lw_radiation, mrt_nlevels, mrt_skip_roof, & 2949 mrt_include_sw, net_radiation, & 2950 radiation_scheme, skip_time_do_radiation, & 2951 sw_radiation, unscheduled_radiation_calls, & 2952 max_raytracing_dist, min_irrf_value, & 2953 nrefsteps, raytrace_mpi_rma, & 2954 surface_reflections, svfnorm_report_thresh, & 2955 radiation_interactions_on, & 2956 rad_angular_discretization, & 2957 raytrace_discrete_azims, & 2958 raytrace_discrete_elevs 2955 NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir, & 2956 albedo_sw_dif, albedo_sw_dir, albedo_type, & 2957 constant_albedo, dt_radiation, emissivity, & 2958 lw_radiation, max_raytracing_dist, & 2959 min_irrf_value, mrt_include_sw, mrt_nlevels, & 2960 mrt_skip_roof, net_radiation, nrefsteps, & 2961 plant_lw_interact, rad_angular_discretization,& 2962 radiation_interactions_on, radiation_scheme, & 2963 raytrace_discrete_azims, & 2964 raytrace_discrete_elevs, raytrace_mpi_rma, & 2965 skip_time_do_radiation, surface_reflections, & 2966 svfnorm_report_thresh, sw_radiation, & 2967 unscheduled_radiation_calls 2959 2968 2960 2969 line = ' ' … … 4648 4657 REAL(wp), DIMENSION(0:nsurf_type) :: costheta !< direct irradiance factor of solar angle 4649 4658 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temperature tendency 4650 REAL(wp), DIMENSION(nzub:nzut) :: pctf_prep !< precalculated factor for canopy transpiration tendency4651 4659 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4652 4660 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 4653 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 4661 REAL(wp) :: asrc !< area of source face 4662 REAL(wp) :: pcrad !< irradiance from plant canopy 4654 4663 REAL(wp) :: pabsswl = 0.0_wp !< total absorbed SW radiation energy in local processor (W) 4655 4664 REAL(wp) :: pabssw = 0.0_wp !< total absorbed SW radiation energy in all processors (W) … … 4672 4681 pchf_prep(:) = r_d * exner(nzub:nzut) & 4673 4682 / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T) 4674 pctf_prep(:) = r_d * exner(nzub:nzut) &4675 / (l_v * hyp(nzub:nzut) * dx*dy*dz(1))4676 4683 ENDIF 4677 4684 #endif … … 4729 4736 mrtinlw(:) = 0._wp 4730 4737 ENDIF 4738 surfinlg(:) = 0._wp !global 4731 4739 4732 4740 … … 4824 4832 ! 4825 4833 !-- For surface-to-surface factors we calculate thermal radiation in 1st pass 4826 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4834 IF ( plant_lw_interact ) THEN 4835 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc) 4836 ELSE 4837 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4838 ENDIF 4827 4839 ENDDO 4828 4840 ENDIF … … 4833 4845 i = surfl(ix, isurf) 4834 4846 surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf) 4835 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 4847 IF ( plant_lw_interact ) THEN 4848 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf) 4849 ELSE 4850 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 4851 ENDIF 4836 4852 ENDDO 4837 4853 ! … … 4880 4896 pcbinswdir(:) = 0._wp 4881 4897 pcbinswdif(:) = 0._wp 4882 pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore4883 ! 4884 !-- 4898 pcbinlw(:) = 0._wp 4899 ! 4900 !-- pcsf first pass 4885 4901 DO icsf = 1, ncsfl 4886 4902 ipcgb = csfsurf(1, icsf) … … 4891 4907 4892 4908 IF ( isurfsrc == -1 ) THEN 4893 !-- Diffuse rad from sky. 4894 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 4895 4896 !--Direct rad 4897 IF ( zenith(0) > 0 ) THEN 4898 !--Estimate directed box absorption 4899 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) 4900 4901 !--isd has already been established, see 1) 4902 pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area & 4903 * pc_abs_frac * dsitransc(ipcgb, isd) 4904 ENDIF 4909 ! 4910 !-- Diffuse rad from sky. 4911 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 4912 ! 4913 !-- Absorbed diffuse LW from sky minus emitted to sky 4914 IF ( plant_lw_interact ) THEN 4915 pcbinlw(ipcgb) = csf(1,icsf) & 4916 * (rad_lw_in_diff(j, i) & 4917 - sigma_sb * (pt(k,j,i)*exner(k))**4) 4918 ENDIF 4919 ! 4920 !-- Direct rad 4921 IF ( zenith(0) > 0 ) THEN 4922 !-- Estimate directed box absorption 4923 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) 4924 ! 4925 !-- isd has already been established, see 1) 4926 pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area & 4927 * pc_abs_frac * dsitransc(ipcgb, isd) 4928 ENDIF 4929 ELSE 4930 IF ( plant_lw_interact ) THEN 4931 ! 4932 !-- Thermal emission from plan canopy towards respective face 4933 pcrad = sigma_sb * (pt(k,j,i) * exner(k))**4 * csf(1,icsf) 4934 surfinlg(isurfsrc) = surfinlg(isurfsrc) + pcrad 4935 ! 4936 !-- Remove the flux above + absorb LW from first pass from surfaces 4937 asrc = facearea(surf(id, isurfsrc)) 4938 pcbinlw(ipcgb) = pcbinlw(ipcgb) & 4939 + (csf(1,icsf) * surfoutl(isurfsrc) & ! Absorb from first pass surf emit 4940 - pcrad) & ! Remove emitted heatflux 4941 * asrc 4942 ENDIF 4905 4943 ENDIF 4906 4944 ENDDO … … 4908 4946 pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:) 4909 4947 ENDIF 4948 4949 IF ( plant_lw_interact ) THEN 4950 ! 4951 !-- Exchange incoming lw radiation from plant canopy 4952 #if defined( __parallel ) 4953 CALL MPI_Allreduce(MPI_IN_PLACE, surfinlg, nsurf, MPI_REAL, MPI_SUM, comm2d, ierr) 4954 IF ( ierr /= 0 ) THEN 4955 WRITE (9,*) 'Error MPI_Allreduce:', ierr 4956 FLUSH(9) 4957 ENDIF 4958 surfinl(:) = surfinl(:) + surfinlg(surfstart(myid)+1:surfstart(myid+1)) 4959 #else 4960 surfinl(:) = surfinl(:) + surfinlg(:) 4961 #endif 4962 ENDIF 4963 4910 4964 surfins = surfinswdir + surfinswdif 4911 4965 surfinl = surfinl + surfinlwdif … … 4966 5020 isurfsrc = svfsurf(2, isvf) 4967 5021 surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 4968 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 5022 IF ( plant_lw_interact ) THEN 5023 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc) 5024 ELSE 5025 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 5026 ENDIF 4969 5027 ENDDO 4970 4971 !-- radiation absorbed by plant canopy 4972 DO icsf = 1, ncsfl 5028 ! 5029 !-- NOTE: PC absorbtion and MRT from reflected can both be done at once 5030 !-- after all reflections if we do one more MPI_ALLGATHERV on surfout. 5031 !-- Advantage: less local computation. Disadvantage: one more collective 5032 !-- MPI call. 5033 ! 5034 !-- Radiation absorbed by plant canopy 5035 DO icsf = 1, ncsfl 4973 5036 ipcgb = csfsurf(1, icsf) 4974 5037 isurfsrc = csfsurf(2, icsf) 4975 5038 IF ( isurfsrc == -1 ) CYCLE ! sky->face only in 1st pass, not here 4976 4977 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) 5039 ! 5040 !-- Calculate source surface area. If the `surf' array is removed 5041 !-- before timestepping starts (future version), then asrc must be 5042 !-- stored within `csf' 5043 asrc = facearea(surf(id, isurfsrc)) 5044 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) * asrc 5045 IF ( plant_lw_interact ) THEN 5046 pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc) * asrc 5047 ENDIF 4978 5048 ENDDO 4979 5049 ! … … 4996 5066 IF ( npcbl > 0 ) THEN 4997 5067 pc_heating_rate(:,:,:) = 0.0_wp 4998 pc_transpiration_rate(:,:,:) = 0.0_wp4999 5068 DO ipcgb = 1, npcbl 5000 5001 5069 j = pcbl(iy, ipcgb) 5002 5070 i = pcbl(ix, ipcgb) 5003 5071 k = pcbl(iz, ipcgb) 5004 5072 ! 5005 !-- 5073 !-- Following expression equals former kk = k - nzb_s_inner(j,i) 5006 5074 kk = k - get_topography_top_index_ji( j, i, 's' ) !- lad arrays are defined flat 5007 5075 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 5008 5076 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt 5009 5010 ! pc_transpiration_rate(kk,j,i) = 0.75_wp* (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &5011 ! * pctf_prep(k) * pt(k, j, i) !-- = dq/dt5012 5013 5077 ENDDO 5078 5079 IF ( humidity .AND. plant_canopy_transpiration ) THEN 5080 !-- Calculation of plant canopy transpiration rate and correspondidng latent heat rate 5081 pc_transpiration_rate(:,:,:) = 0.0_wp 5082 pc_latent_rate(:,:,:) = 0.0_wp 5083 DO ipcgb = 1, npcbl 5084 i = pcbl(ix, ipcgb) 5085 j = pcbl(iy, ipcgb) 5086 k = pcbl(iz, ipcgb) 5087 kk = k - get_topography_top_index_ji( j, i, 's' ) !- lad arrays are defined flat 5088 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), & 5089 pc_transpiration_rate(kk,j,i), pc_latent_rate(kk,j,i) ) 5090 ENDDO 5091 ENDIF 5014 5092 ENDIF 5015 5093 ! … … 5161 5239 !-- domain when using average_radiation in the respective radiation model 5162 5240 5163 !-- Precalculate face areas for all face directions using normal vector5164 DO d = 0, nsurf_type5165 facearea(d) = 1._wp5166 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx5167 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy5168 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)5169 ENDDO5170 5241 !-- calculate horizontal area 5171 5242 ! !!! ATTENTION!!! uniform grid is assumed here … … 5429 5500 IMPLICIT NONE 5430 5501 5431 INTEGER(iwp) :: i, j, k, l, m 5502 INTEGER(iwp) :: i, j, k, l, m, d 5432 5503 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 5433 5504 INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt … … 5439 5510 #endif 5440 5511 5512 ! 5513 !-- precalculate face areas for different face directions using normal vector 5514 DO d = 0, nsurf_type 5515 facearea(d) = 1._wp 5516 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 5517 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 5518 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1) 5519 ENDDO 5441 5520 ! 5442 5521 !-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be … … 5904 5983 ALLOCATE( surfouts(nsurf) ) 5905 5984 ALLOCATE( surfoutl(nsurf) ) 5985 ALLOCATE( surfinlg(nsurf) ) 5906 5986 ALLOCATE( skyvf(nsurfl) ) 5907 5987 ALLOCATE( skyvft(nsurfl) ) … … 5940 6020 REAL(wp) :: az1, az2 !< relative azimuth of section borders 5941 6021 REAL(wp) :: azmid !< ray (center) azimuth 6022 REAL(wp) :: yxlen !< |yxdir| 5942 6023 REAL(wp), DIMENSION(2) :: yxdir !< y,x *unit* vector of ray direction (in grid units) 5943 6024 REAL(wp), DIMENSION(:), ALLOCATABLE :: zdirs !< directions in z (tangent of elevation) 6025 REAL(wp), DIMENSION(:), ALLOCATABLE :: zcent !< zenith angle centers 5944 6026 REAL(wp), DIMENSION(:), ALLOCATABLE :: zbdry !< zenith angle boundaries 5945 6027 REAL(wp), DIMENSION(:), ALLOCATABLE :: vffrac !< view factor fractions for individual rays … … 5949 6031 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: itarget !< face indices of detected obstacles 5950 6032 INTEGER(iwp) :: itarg0, itarg1 5951 #if defined( __parallel ) 5952 #endif 5953 5954 5955 5956 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 6033 5957 6034 INTEGER(iwp) :: udim 5958 6035 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l … … 5966 6043 LOGICAL :: visible 5967 6044 REAL(wp), DIMENSION(3) :: sa, ta !< real coordinates z,y,x of source and target 6045 REAL(wp) :: difvf !< differential view factor 5968 6046 REAL(wp) :: transparency, rirrf, sqdist, svfsum 5969 6047 INTEGER(iwp) :: isurflt, isurfs, isurflt_prev … … 5983 6061 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 5984 6062 CALL radiation_write_debug_log('Start calculation of SVF and CSF') 5985 5986 !-- precalculate face areas for different face directions using normal vector5987 DO d = 0, nsurf_type5988 facearea(d) = 1._wp5989 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx5990 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy5991 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)5992 ENDDO5993 6063 5994 6064 !-- initialize variables and temporary arrays for calculation of svf and csf … … 6219 6289 END SELECT 6220 6290 6221 ALLOCATE ( zdirs(1:nzn), z bdry(0:nzn), vffrac(1:nzn*naz), &6291 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), & 6222 6292 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6223 6293 !in case of rad_angular_discretization … … 6225 6295 itarg0 = 1 6226 6296 itarg1 = nzn 6227 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6297 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6228 6298 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 6229 6299 IF ( td == iup_u .OR. td == iup_l ) THEN … … 6251 6321 !-- sum of whole vffrac equals 1, verified 6252 6322 ENDIF 6253 yxdir = (/ COS(azmid), SIN(azmid) /) 6323 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6324 yxlen = SQRT(SUM(yxdir(:)**2)) 6325 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6326 yxdir(:) = yxdir(:) / yxlen 6327 6254 6328 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6255 6329 surfstart(myid) + isurflt, facearea(td), & … … 6257 6331 .FALSE., lowest_free_ray, & 6258 6332 ztransp(itarg0:itarg1), & 6259 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6260 !FIXME itarget available only in 6261 !case of rad_angular_discretization 6333 itarget(itarg0:itarg1)) 6334 6262 6335 skyvf(isurflt) = skyvf(isurflt) + & 6263 6336 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) … … 6337 6410 ENDIF ! rad_angular_discretization 6338 6411 6339 DEALLOCATE ( zdirs, z bdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only6412 DEALLOCATE ( zdirs, zcent, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only 6340 6413 !in case of rad_angular_discretization 6341 6414 ! … … 6367 6440 ENDIF 6368 6441 6442 difvf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 6443 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 6444 / (pi * sqdist) ! square of distance between centers 6445 ! 6369 6446 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 6370 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 6371 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 6372 / (pi * sqdist) & ! square of distance between centers 6373 * facearea(sd) 6447 rirrf = difvf * facearea(sd) 6374 6448 6375 6449 !-- reject raytracing for potentially too small view factor values … … 6380 6454 6381 6455 !-- raytrace + process plant canopy sinks within 6382 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &6456 CALL raytrace(sa, ta, isurfs, difvf, facearea(td), .TRUE., & 6383 6457 visible, transparency) 6384 6458 … … 6428 6502 nzn = raytrace_discrete_elevs / 2 6429 6503 zns = pi / 2._wp / REAL(nzn, wp) 6430 ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn), &6504 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), vffrac(1:nzn), ztransp(1:nzn), & 6431 6505 itarget(1:nzn) ) 6432 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6506 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6433 6507 vffrac(:) = 0._wp 6434 6508 … … 6440 6514 DO iaz = 1, naz 6441 6515 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 6442 yxdir = (/ COS(azmid), SIN(azmid) /) 6516 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6517 yxlen = SQRT(SUM(yxdir(:)**2)) 6518 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6519 yxdir(:) = yxdir(:) / yxlen 6443 6520 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6444 6521 -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., & 6445 lowest_free_ray, ztransp, itarget) !FIXME unit vect in grid units + zdirs6522 lowest_free_ray, ztransp, itarget) 6446 6523 6447 6524 !-- Save direct solar transparency … … 6456 6533 ENDDO 6457 6534 ENDDO 6458 DEALLOCATE ( zdirs, vffrac, ztransp, itarget )6535 DEALLOCATE ( zdirs, zcent, vffrac, ztransp, itarget ) 6459 6536 !-- 6460 6537 !-- Raytrace to MRT boxes … … 6469 6546 nzn = raytrace_discrete_elevs 6470 6547 zns = pi / REAL(nzn, wp) 6471 ALLOCATE ( zdirs(1:nzn), z bdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &6548 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), & 6472 6549 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6473 6550 !in case of rad_angular_discretization 6474 6551 6475 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6552 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6476 6553 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 6477 6554 vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp) … … 6493 6570 DO iaz = 1, naz 6494 6571 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 6495 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), nzn, zdirs, & 6572 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6573 yxlen = SQRT(SUM(yxdir(:)**2)) 6574 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6575 yxdir(:) = yxdir(:) / yxlen 6576 6577 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6496 6578 -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., & 6497 6579 .FALSE., .TRUE., lowest_free_ray, & 6498 6580 ztransp(itarg0:itarg1), & 6499 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6500 !FIXME itarget available only in 6501 !case of rad_angular_discretization 6581 itarget(itarg0:itarg1)) 6502 6582 6503 6583 !-- Sky view factors for MRT … … 6574 6654 6575 6655 ENDDO ! imrt 6576 DEALLOCATE ( zdirs, z bdry, vffrac, vffrac0, ztransp, itarget )6656 DEALLOCATE ( zdirs, zcent, zbdry, vffrac, vffrac0, ztransp, itarget ) 6577 6657 ! 6578 6658 !-- Move MRT factors to final arrays … … 6778 6858 j = 0 6779 6859 ENDIF 6780 !-- fill out real values of rsvf, rtransp 6781 csflt(1,kcsf) = acsf(kcsf)%rsvf 6860 csflt(1,kcsf) = acsf(kcsf)%rcvf 6782 6861 !-- fill out integer values of itz,ity,itx,isurfs 6783 6862 kcsflt(1,kcsf) = acsf(kcsf)%itz … … 6946 7025 !> doesn't need to be the same in all three dimensions). 6947 7026 !------------------------------------------------------------------------------! 6948 SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency)7027 SUBROUTINE raytrace(src, targ, isrc, difvf, atarg, create_csf, visible, transparency) 6949 7028 IMPLICIT NONE 6950 7029 6951 7030 REAL(wp), DIMENSION(3), INTENT(in) :: src, targ !< real coordinates z,y,x 6952 7031 INTEGER(iwp), INTENT(in) :: isrc !< index of source face for csf 6953 REAL(wp), INTENT(in) :: rirrf !< irradiancefactor for csf7032 REAL(wp), INTENT(in) :: difvf !< differential view factor for csf 6954 7033 REAL(wp), INTENT(in) :: atarg !< target surface area for csf 6955 7034 LOGICAL, INTENT(in) :: create_csf !< whether to generate new CSFs during raytracing … … 7113 7192 acsf(ncsfl)%itz = boxes(1,i) 7114 7193 acsf(ncsfl)%isurfs = isrc 7115 acsf(ncsfl)%r svf = cursink*transparency*rirrf*atarg7194 acsf(ncsfl)%rcvf = cursink*transparency*difvf*atarg 7116 7195 ENDIF !< create_csf 7117 7196 … … 7514 7593 acsf(ncsfl)%itz = iz 7515 7594 acsf(ncsfl)%isurfs = iorig 7516 acsf(ncsfl)%r svf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)7595 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k) 7517 7596 ENDIF 7518 7597 … … 7571 7650 acsf(ncsfl)%ity = rt2_track(1,i) 7572 7651 acsf(ncsfl)%itz = iz 7573 acsf(ncsfl)%isurfs = itarget(k) !if above horizon, then itarget(k)==-1, which7574 !is also a special ID indicating sky7575 acsf(ncsfl)%r svf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)7652 IF ( itarget(k) /= -1 ) ERROR STOP !FIXME remove after test 7653 acsf(ncsfl)%isurfs = -1 7654 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k) 7576 7655 ENDIF ! create_csf 7577 7656 … … 7600 7679 INTEGER(iwp), TARGET, INTENT(out) :: isurfl 7601 7680 INTEGER(iwp), INTENT(out) :: iproc 7681 #if defined( __parallel ) 7682 #else 7683 INTEGER(iwp) :: target_displ !< index of the grid in the local gridsurf array 7684 #endif 7602 7685 INTEGER(iwp) :: px, py !< number of processors in x and y direction 7603 7686 !< before the processor in the question … … 8138 8221 .AND. acsfnew(iwrite)%isurfs == acsf(iread)%isurfs ) THEN 8139 8222 8140 acsfnew(iwrite)%r svf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf8223 acsfnew(iwrite)%rcvf = acsfnew(iwrite)%rcvf + acsf(iread)%rcvf 8141 8224 !-- advance reading index, keep writing index 8142 8225 ELSE -
palm/trunk/SOURCE/urban_surface_mod.f90
r3435 r3449 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Bugfix: Fix average arrays allocations in usm_average_3d_data (J.Resler) 31 ! Bugfix: Fix reading wall temperatures (J.Resler) 32 ! Bugfix: Fix treating of outputs for wall temperature and sky view factors (J.Resler) 33 ! 34 ! 35 ! 3435 2018-10-26 18:25:44Z gronemeier 30 36 ! Bugfix: allocate gamma_w_green_sat until nzt_wall+1 31 37 ! … … 411 417 412 418 USE plant_canopy_model_mod, & 413 ONLY: pc_heating_rate, pc_transpiration_rate 419 ONLY: pc_heating_rate, pc_transpiration_rate, pc_latent_rate 414 420 415 421 USE radiation_model_mod, & … … 1845 1851 CASE ( 'usm_rad_net' ) 1846 1852 !-- array of complete radiation balance 1847 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN 1848 ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) ) 1849 surf_usm_h%rad_net_av = 0.0_wp 1853 IF ( l == -1 ) THEN 1854 IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN 1855 ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) ) 1856 surf_usm_h%rad_net_av = 0.0_wp 1857 ENDIF 1850 1858 ELSE 1851 1859 IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN … … 1968 1976 CASE ( 'usm_rad_hf' ) 1969 1977 !-- array of heat flux from radiation for surfaces after i-th reflection 1970 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN 1971 ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) ) 1972 surf_usm_h%surfhf_av = 0.0_wp 1978 IF ( l == -1 ) THEN 1979 IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN 1980 ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) ) 1981 surf_usm_h%surfhf_av = 0.0_wp 1982 ENDIF 1973 1983 ELSE 1974 1984 IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN … … 1981 1991 !-- array of sensible heat flux from surfaces 1982 1992 !-- land surfaces 1983 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN 1984 ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) ) 1985 surf_usm_h%wshf_eb_av = 0.0_wp 1993 IF ( l == -1 ) THEN 1994 IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN 1995 ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) ) 1996 surf_usm_h%wshf_eb_av = 0.0_wp 1997 ENDIF 1986 1998 ELSE 1987 1999 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN … … 2036 2048 CASE ( 'usm_wghf' ) 2037 2049 !-- array of heat flux from ground (wall, roof, land) 2038 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN 2039 ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) ) 2040 surf_usm_h%wghf_eb_av = 0.0_wp 2050 IF ( l == -1 ) THEN 2051 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN 2052 ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) ) 2053 surf_usm_h%wghf_eb_av = 0.0_wp 2054 ENDIF 2041 2055 ELSE 2042 2056 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN … … 2048 2062 CASE ( 'usm_wghf_window' ) 2049 2063 !-- array of heat flux from window ground (wall, roof, land) 2050 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN 2051 ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) ) 2052 surf_usm_h%wghf_eb_window_av = 0.0_wp 2064 IF ( l == -1 ) THEN 2065 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN 2066 ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) ) 2067 surf_usm_h%wghf_eb_window_av = 0.0_wp 2068 ENDIF 2053 2069 ELSE 2054 2070 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN … … 2060 2076 CASE ( 'usm_wghf_green' ) 2061 2077 !-- array of heat flux from green ground (wall, roof, land) 2062 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN 2063 ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) ) 2064 surf_usm_h%wghf_eb_green_av = 0.0_wp 2078 IF ( l == -1 ) THEN 2079 IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN 2080 ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) ) 2081 surf_usm_h%wghf_eb_green_av = 0.0_wp 2082 ENDIF 2065 2083 ELSE 2066 2084 IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN … … 2072 2090 CASE ( 'usm_iwghf' ) 2073 2091 !-- array of heat flux from indoor ground (wall, roof, land) 2074 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN 2075 ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) ) 2076 surf_usm_h%iwghf_eb_av = 0.0_wp 2092 IF ( l == -1 ) THEN 2093 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN 2094 ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) ) 2095 surf_usm_h%iwghf_eb_av = 0.0_wp 2096 ENDIF 2077 2097 ELSE 2078 2098 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN … … 2084 2104 CASE ( 'usm_iwghf_window' ) 2085 2105 !-- array of heat flux from indoor window ground (wall, roof, land) 2086 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN 2087 ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) ) 2088 surf_usm_h%iwghf_eb_window_av = 0.0_wp 2106 IF ( l == -1 ) THEN 2107 IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN 2108 ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) ) 2109 surf_usm_h%iwghf_eb_window_av = 0.0_wp 2110 ENDIF 2089 2111 ELSE 2090 2112 IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN … … 2093 2115 ENDIF 2094 2116 ENDIF 2095 2117 2096 2118 CASE ( 'usm_t_surf_wall' ) 2097 2119 !-- surface temperature for surfaces 2098 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_wall_av) ) THEN 2099 ALLOCATE( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) ) 2100 surf_usm_h%t_surf_wall_av = 0.0_wp 2120 IF ( l == -1 ) THEN 2121 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_wall_av) ) THEN 2122 ALLOCATE( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) ) 2123 surf_usm_h%t_surf_wall_av = 0.0_wp 2124 ENDIF 2101 2125 ELSE 2102 2126 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_wall_av) ) THEN … … 2108 2132 CASE ( 'usm_t_surf_window' ) 2109 2133 !-- surface temperature for window surfaces 2110 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN 2111 ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) ) 2112 surf_usm_h%t_surf_window_av = 0.0_wp 2134 IF ( l == -1 ) THEN 2135 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN 2136 ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) ) 2137 surf_usm_h%t_surf_window_av = 0.0_wp 2138 ENDIF 2113 2139 ELSE 2114 2140 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN … … 2120 2146 CASE ( 'usm_t_surf_green' ) 2121 2147 !-- surface temperature for green surfaces 2122 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN 2123 ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) ) 2124 surf_usm_h%t_surf_green_av = 0.0_wp 2148 IF ( l == -1 ) THEN 2149 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN 2150 ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) ) 2151 surf_usm_h%t_surf_green_av = 0.0_wp 2152 ENDIF 2125 2153 ELSE 2126 2154 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN … … 2132 2160 CASE ( 'usm_t_surf_10cm' ) 2133 2161 !-- near surface temperature for whole surfaces 2134 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN 2135 ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) ) 2136 surf_usm_h%t_surf_10cm_av = 0.0_wp 2162 IF ( l == -1 ) THEN 2163 IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN 2164 ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) ) 2165 surf_usm_h%t_surf_10cm_av = 0.0_wp 2166 ENDIF 2137 2167 ELSE 2138 2168 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) ) THEN … … 2144 2174 CASE ( 'usm_t_wall' ) 2145 2175 !-- wall temperature for iwl layer of walls and land 2146 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN 2147 ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2148 surf_usm_h%t_wall_av = 0.0_wp 2176 IF ( l == -1 ) THEN 2177 IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN 2178 ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2179 surf_usm_h%t_wall_av = 0.0_wp 2180 ENDIF 2149 2181 ELSE 2150 2182 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN … … 2156 2188 CASE ( 'usm_t_window' ) 2157 2189 !-- window temperature for iwl layer of walls and land 2158 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN 2159 ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2160 surf_usm_h%t_window_av = 0.0_wp 2190 IF ( l == -1 ) THEN 2191 IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN 2192 ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2193 surf_usm_h%t_window_av = 0.0_wp 2194 ENDIF 2161 2195 ELSE 2162 2196 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN … … 2168 2202 CASE ( 'usm_t_green' ) 2169 2203 !-- green temperature for iwl layer of walls and land 2170 IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN 2171 ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2172 surf_usm_h%t_green_av = 0.0_wp 2204 IF ( l == -1 ) THEN 2205 IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN 2206 ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) 2207 surf_usm_h%t_green_av = 0.0_wp 2208 ENDIF 2173 2209 ELSE 2174 2210 IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN … … 3091 3127 CHARACTER(LEN=*),INTENT(OUT) :: unit !< 3092 3128 3093 INTEGER(iwp) :: i,j !< index 3129 INTEGER(iwp) :: i,j,l !< index 3130 CHARACTER(LEN=2) :: ls 3094 3131 CHARACTER(LEN=varnamelength) :: var !< TRIM(variable) 3095 INTEGER(iwp), PARAMETER :: nl1 = 3 2!< number of directional usm variables3132 INTEGER(iwp), PARAMETER :: nl1 = 31 !< number of directional usm variables 3096 3133 CHARACTER(LEN=varnamelength), DIMENSION(nl1) :: varlist1 = & !< list of directional usm variables 3097 3134 (/'usm_rad_net ', & … … 3119 3156 'usm_surfalb ', & 3120 3157 'usm_surfemis ', & 3121 'usm_t_surf 3158 'usm_t_surf_wall ', & 3122 3159 'usm_t_surf_window ', & 3123 3160 'usm_t_surf_green ', & 3124 3161 'usm_t_green ', & 3125 3162 'usm_t_surf_10cm ', & 3126 'usm_t_wall ', & 3127 'usm_t_window ', & 3128 'usm_t_green '/) 3129 3130 INTEGER(iwp), PARAMETER :: nl2 = 9 !< number of other variables 3163 'usm_skyvf ', & 3164 'usm_skyvft '/) 3165 3166 INTEGER(iwp), PARAMETER :: nl2 = 7 !< number of other variables 3131 3167 CHARACTER(LEN=varnamelength), DIMENSION(nl2) :: varlist2 = & !< list of other usm variables 3132 (/'usm_skyvf ', & 3133 'usm_skyvft ', & 3134 'usm_svf ', & 3168 (/'usm_svf ', & 3135 3169 'usm_dif ', & 3136 3170 'usm_rad_pc_inlw ', & … … 3139 3173 'usm_rad_pc_inswdif ', & 3140 3174 'usm_rad_pc_inswref '/) 3175 3176 INTEGER(iwp), PARAMETER :: nl3 = 3 !< number of directional layer usm variables 3177 CHARACTER(LEN=varnamelength), DIMENSION(nl3) :: varlist3 = & !< list of directional layer usm variables 3178 (/'usm_t_wall ', & 3179 'usm_t_window ', & 3180 'usm_t_green '/) 3141 3181 3142 3182 INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions … … 3158 3198 EXIT 3159 3199 ENDIF 3200 IF ( lfound ) EXIT 3201 ENDDO 3202 ENDDO 3203 IF ( lfound ) GOTO 10 3204 ! directional layer variables 3205 DO i = 1, nl3 3206 DO j = 1, nd 3207 DO l = nzb_wall, nzt_wall 3208 WRITE(ls,'(A1,I1)') '_',l 3209 IF ( TRIM(var) == TRIM(varlist3(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN 3210 lfound = .TRUE. 3211 EXIT 3212 ENDIF 3213 ENDDO 3160 3214 IF ( lfound ) EXIT 3161 3215 ENDDO … … 4503 4557 var(1:9) == 'usm_qsws_' .OR. var(1:13) == 'usm_qsws_veg_' .OR. & 4504 4558 var(1:13) == 'usm_qsws_liq_' .OR. & 4505 var(1:1 0) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. &4559 var(1:15) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. & 4506 4560 var(1:17) == 'usm_t_surf_window' .OR. var(1:12) == 'usm_t_window' .OR. & 4507 4561 var(1:16) == 'usm_t_surf_green' .OR. var(1:11) == 'usm_t_green' .OR. & … … 5910 5964 ENDIF 5911 5965 5912 IF ( plant_canopy ) THEN5913 5914 IF ( .NOT. ALLOCATED( pc_heating_rate) ) THEN5915 !-- then pc_heating_rate is allocated in init_plant_canopy5916 !-- in case of cthf /= 0 => we need to allocate it for our use here5917 ALLOCATE( pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )5918 5919 pc_heating_rate = 0.0_wp5920 5921 ENDIF5922 5923 IF ( .NOT. ALLOCATED( pc_transpiration_rate) ) THEN5924 !-- then pc_heating_rate is allocated in init_plant_canopy5925 !-- in case of cthf /= 0 => we need to allocate it for our use here5926 ALLOCATE( pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )5927 5928 pc_transpiration_rate = 0.0_wp5929 5930 5931 ENDIF5932 ENDIF5933 5966 ! 5934 5967 !-- Check for consistent initialization. … … 6092 6125 ENDDO 6093 6126 ENDDO 6094 ELSE6095 !-- If specified, replace constant wall temperatures with fully 3D values from file6096 IF ( read_wall_temp_3d ) CALL usm_read_wall_temperature()6097 !6098 6127 ENDIF 6099 6100 !-- 6128 6129 !-- If specified, replace constant wall temperatures with fully 3D values from file 6130 IF ( read_wall_temp_3d ) CALL usm_read_wall_temperature() 6131 6132 !-- 6101 6133 !-- Possibly DO user-defined actions (e.g. define heterogeneous wall surface) 6102 6134 CALL user_init_urban_surface … … 8870 8902 !> x, y, d, z coordinates 8871 8903 !------------------------------------------------------------------------------! 8872 PURE FUNCTION advance_surface(isurfl_start, isurfl_stop, x, y, z, d) & 8873 result(isurfl) 8874 8875 INTEGER(iwp), INTENT(in) :: isurfl_start, isurfl_stop 8904 PURE FUNCTION find_surface( x, y, z, d ) result(isurfl) 8905 8876 8906 INTEGER(iwp), INTENT(in) :: x, y, z, d 8877 INTEGER(iwp) :: isx, isy, isz, isd8878 8907 INTEGER(iwp) :: isurfl 8879 8880 DO isurfl = isurfl_start, isurfl_stop 8881 isx = surfl(ix, isurfl) 8882 isy = surfl(iy, isurfl) 8883 isz = surfl(iz, isurfl) 8884 isd = surfl(id, isurfl) 8885 IF ( isx==x .and. isy==y .and. isz==z .and. isd==d ) RETURN 8886 ENDDO 8908 INTEGER(iwp) :: isx, isy, isz 8909 8910 IF ( d == 0 ) THEN 8911 DO isurfl = 1, surf_usm_h%ns 8912 isx = surf_usm_h%i(isurfl) 8913 isy = surf_usm_h%j(isurfl) 8914 isz = surf_usm_h%k(isurfl) 8915 IF ( isx==x .and. isy==y .and. isz==z ) RETURN 8916 ENDDO 8917 ELSE 8918 DO isurfl = 1, surf_usm_v(d-1)%ns 8919 isx = surf_usm_v(d-1)%i(isurfl) 8920 isy = surf_usm_v(d-1)%j(isurfl) 8921 isz = surf_usm_v(d-1)%k(isurfl) 8922 IF ( isx==x .and. isy==y .and. isz==z ) RETURN 8923 ENDDO 8924 ENDIF 8887 8925 8888 8926 !-- coordinate not found … … 8927 8965 nys <= j .and. j <= nyn) THEN !< local processor 8928 8966 !-- identify surface id 8929 isurfl = advance_surface(isurfl+1, nsurfl, i, j, k, d)8967 isurfl = find_surface( i, j, k, d ) 8930 8968 IF ( isurfl == -1 ) THEN 8931 8969 WRITE(message_string, '(a,4i5,a,i5,a)') 'Coordinates (xyzd) ', i, j, k, d, & … … 8939 8977 t_surf_wall_h(isurfl) = rtsurf 8940 8978 t_wall_h(:,isurfl) = rtwall(:) 8979 t_window_h(:,isurfl) = rtwall(:) 8980 t_green_h(:,isurfl) = rtwall(:) 8941 8981 ELSE 8942 8982 t_surf_wall_v(d-1)%t(isurfl) = rtsurf 8943 8983 t_wall_v(d-1)%t(:,isurfl) = rtwall(:) 8984 t_window_v(d-1)%t(:,isurfl) = rtwall(:) 8985 t_green_v(d-1)%t(:,isurfl) = rtwall(:) 8944 8986 ENDIF 8945 8987 ENDIF … … 9663 9705 surf_usm_v(l)%rad_lw_out(m) 9664 9706 9665 if ((abs(t_surf_wall_v(l)%t(m)-276.).gt.100.).or.(abs(t_surf_window_v(l)%t(m)-276.).gt.300.) &9666 .or.(abs(t_surf_green_v(l)%t(m)-276.).gt.100)) then9667 print*, "tsurfvvvv",m,t_surf_wall_v(l)%t(m),t_surf_window_v(l)%t(m),t_surf_green_v(l)%t(m)9668 print*, "params", surf_usm_v(l)%emissivity(ind_veg_wall,m),lambda_surface, t_wall_v(l)%t(nzb_wall:nzt_wall,m), &9669 surf_usm_v(l)%emissivity(ind_wat_win,m),lambda_surface_window,t_window_v(l)%t(nzb_wall:nzt_wall,m),t_green_v(l)%t(nzb_wall:nzt_wall,m)9670 print*, "dicken",surf_usm_v(l)%zw(nzb_wall:nzt_wall,m),surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,m),surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,m)9671 print*, "c",surf_usm_v(l)%c_surface_window(m),surf_usm_v(l)%c_surface(m)9672 !if ((abs(t_surf_v(l)%t(m)-276.).gt.10.).or.(abs(t_surf_window_v(l)%t(m)-276.).gt.10.)) then9673 stop9674 endif9675 9707 !-- numerator of the prognostic equation 9676 9708 coef_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout included in calculation of radnet_l
Note: See TracChangeset
for help on using the changeset viewer.
