Changeset 4717 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- Sep 30, 2020 10:27:40 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/urban_surface_mod.f90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/urban_surface_mod.f90
r4713 r4717 27 27 ! ----------------- 28 28 ! $Id$ 29 ! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP 30 ! directives (J. Resler) 31 ! 32 ! 4713 2020-09-29 12:02:05Z pavelkrc 29 33 ! - Do not change original fractions in USM energy balance 30 34 ! - Correct OpenMP parallelization … … 5159 5163 ENDIF 5160 5164 5161 !$OMP PARALLEL PRIVATE (m, i, j, k, kw, wtend, wintend, win_absorp, wall_mod)5162 5165 wall_mod=1.0_wp 5163 5166 IF ( usm_wall_mod .AND. during_spinup ) THEN … … 5184 5187 ! 5185 5188 !-- Cycle for all surfaces in given direction 5186 !$OMP DOSCHEDULE (STATIC)5189 !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw, wtend, wintend, win_absorp, wall_mod) SCHEDULE (STATIC) 5187 5190 DO m = 1, surf%ns 5188 5191 ! … … 5194 5197 !-- Prognostic equation for ground/roof temperature t_wall 5195 5198 wtend(:) = 0.0_wp 5196 wtend(nzb_wall) = ( 1.0_wp / surf%rho_c_wall(nzb_wall,m) ) &5197 * ( surf%lambda_h(nzb_wall,m) * wall_mod(nzb_wall) &5198 * ( t_wall%val(nzb_wall+1,m) - t_wall%val(nzb_wall,m) ) &5199 * surf%ddz_wall(nzb_wall+1,m) &5200 + surf%frac(m,ind_veg_wall) &5201 / ( surf%frac(m,ind_veg_wall) &5202 + surf%frac(m,ind_pav_green) ) &5203 * surf%wghf_eb(m) &5204 - surf%frac(m,ind_pav_green) &5205 / ( surf%frac(m,ind_veg_wall) &5206 + surf%frac(m,ind_pav_green) ) &5207 * ( surf%lambda_h_green(nzt_wall,m) &5208 * wall_mod(nzt_wall) &5209 * surf%ddz_green(nzt_wall,m) &5210 + surf%lambda_h(nzb_wall,m) &5211 * wall_mod(nzb_wall) &5212 * surf%ddz_wall(nzb_wall,m) ) &5213 / ( surf%ddz_green(nzt_wall,m) &5214 + surf%ddz_wall(nzb_wall,m) ) &5215 * ( t_wall%val(nzb_wall,m) - t_green%val(nzt_wall,m) ) &5199 wtend(nzb_wall) = ( 1.0_wp / surf%rho_c_wall(nzb_wall,m) ) & 5200 * ( surf%lambda_h(nzb_wall,m) * wall_mod(nzb_wall) & 5201 * ( t_wall%val(nzb_wall+1,m) - t_wall%val(nzb_wall,m) ) & 5202 * surf%ddz_wall(nzb_wall+1,m) & 5203 + surf%frac(m,ind_veg_wall) & 5204 / ( surf%frac(m,ind_veg_wall) & 5205 + surf%frac(m,ind_pav_green) ) & 5206 * surf%wghf_eb(m) & 5207 - surf%frac(m,ind_pav_green) & 5208 / ( surf%frac(m,ind_veg_wall) & 5209 + surf%frac(m,ind_pav_green) ) & 5210 * ( surf%lambda_h_green(nzt_wall,m) & 5211 * wall_mod(nzt_wall) & 5212 * surf%ddz_green(nzt_wall,m) & 5213 + surf%lambda_h(nzb_wall,m) & 5214 * wall_mod(nzb_wall) & 5215 * surf%ddz_wall(nzb_wall,m) ) & 5216 / ( surf%ddz_green(nzt_wall,m) & 5217 + surf%ddz_wall(nzb_wall,m) ) & 5218 * ( t_wall%val(nzb_wall,m) - t_green%val(nzt_wall,m) ) & 5216 5219 ) * surf%ddz_wall_stag(nzb_wall,m) 5217 5220 … … 5226 5229 5227 5230 DO kw = nzb_wall+1, nzt_wall-1 5228 wtend(kw) = ( 1.0_wp / surf%rho_c_wall(kw,m) ) &5229 * ( surf%lambda_h(kw,m) * wall_mod(kw)&5230 * ( t_wall%val(kw+1,m) - t_wall%val(kw,m) )&5231 * surf%ddz_wall(kw+1,m)&5232 - surf%lambda_h(kw-1,m) * wall_mod(kw-1)&5233 * ( t_wall%val(kw,m) - t_wall%val(kw-1,m) )&5234 * surf%ddz_wall(kw,m)&5235 5231 wtend(kw) = ( 1.0_wp / surf%rho_c_wall(kw,m) ) & 5232 * ( surf%lambda_h(kw,m) * wall_mod(kw) & 5233 * ( t_wall%val(kw+1,m) - t_wall%val(kw,m) ) & 5234 * surf%ddz_wall(kw+1,m) & 5235 - surf%lambda_h(kw-1,m) * wall_mod(kw-1) & 5236 * ( t_wall%val(kw,m) - t_wall%val(kw-1,m) ) & 5237 * surf%ddz_wall(kw,m) & 5238 ) * surf%ddz_wall_stag(kw,m) 5236 5239 ENDDO 5237 5238 wtend(nzt_wall) = ( 1.0_wp / surf%rho_c_wall(nzt_wall,m) ) & 5239 * ( -surf%lambda_h(nzt_wall-1,m) * wall_mod(nzt_wall-1) & 5240 * ( t_wall%val(nzt_wall,m) - t_wall%val(nzt_wall-1,m) ) & 5241 * surf%ddz_wall(nzt_wall,m) & 5242 + surf%iwghf_eb(m) & 5243 ) * surf%ddz_wall_stag(nzt_wall,m) 5244 5245 t_wall_p%val(nzb_wall:nzt_wall,m) = t_wall%val(nzb_wall:nzt_wall,m) + dt_3d & 5246 * ( tsc(2) * wtend(nzb_wall:nzt_wall) + tsc(3) & 5240 wtend(nzt_wall) = ( 1.0_wp / surf%rho_c_wall(nzt_wall,m) ) & 5241 * ( -surf%lambda_h(nzt_wall-1,m) * wall_mod(nzt_wall-1) & 5242 * ( t_wall%val(nzt_wall,m) - t_wall%val(nzt_wall-1,m) ) & 5243 * surf%ddz_wall(nzt_wall,m) + surf%iwghf_eb(m) & 5244 ) * surf%ddz_wall_stag(nzt_wall,m) 5245 5246 t_wall_p%val(nzb_wall:nzt_wall,m) = t_wall%val(nzb_wall:nzt_wall,m) + dt_3d & 5247 * ( tsc(2) * wtend(nzb_wall:nzt_wall) + tsc(3) & 5247 5248 * surf%tt_wall_m(nzb_wall:nzt_wall,m) ) 5248 5249 … … 5255 5256 !-- of shortwave radiation into account 5256 5257 wintend(:) = 0.0_wp 5257 wintend(nzb_wall) = ( 1.0_wp / surf%rho_c_window(nzb_wall,m) ) &5258 * ( surf%lambda_h_window(nzb_wall,m) &5259 * ( t_window%val(nzb_wall+1,m) - t_window%val(nzb_wall,m) ) &5260 * surf%ddz_window(nzb_wall+1,m) &5261 + surf%wghf_eb_window(m) &5262 + surf%rad_sw_in(m) &5263 * ( 1.0_wp - exp( -win_absorp &5264 * surf%zw_window(nzb_wall,m) ) ) &5258 wintend(nzb_wall) = ( 1.0_wp / surf%rho_c_window(nzb_wall,m) ) & 5259 * ( surf%lambda_h_window(nzb_wall,m) & 5260 * ( t_window%val(nzb_wall+1,m) - t_window%val(nzb_wall,m) ) & 5261 * surf%ddz_window(nzb_wall+1,m) & 5262 + surf%wghf_eb_window(m) & 5263 + surf%rad_sw_in(m) & 5264 * ( 1.0_wp - exp( -win_absorp & 5265 * surf%zw_window(nzb_wall,m) ) ) & 5265 5266 ) * surf%ddz_window_stag(nzb_wall,m) 5266 5267 … … 5272 5273 ENDIF 5273 5274 5275 5274 5276 DO kw = nzb_wall+1, nzt_wall-1 5275 wintend(kw) = ( 1.0_wp / surf%rho_c_window(kw,m) ) &5276 * ( surf%lambda_h_window(kw,m)&5277 * ( t_window%val(kw+1,m) - t_window%val(kw,m) )&5278 * surf%ddz_window(kw+1,m)&5279 - surf%lambda_h_window(kw-1,m)&5280 * ( t_window%val(kw,m) - t_window%val(kw-1,m) )&5281 * surf%ddz_window(kw,m)&5282 + surf%rad_sw_in(m)&5283 * ( exp( -win_absorp * surf%zw_window(kw-1,m) )&5284 - exp(-win_absorp * surf%zw_window(kw,m) )&5285 )&5286 5277 wintend(kw) = ( 1.0_wp / surf%rho_c_window(kw,m) ) & 5278 * ( surf%lambda_h_window(kw,m) & 5279 * ( t_window%val(kw+1,m) - t_window%val(kw,m) ) & 5280 * surf%ddz_window(kw+1,m) & 5281 - surf%lambda_h_window(kw-1,m) & 5282 * ( t_window%val(kw,m) - t_window%val(kw-1,m) ) & 5283 * surf%ddz_window(kw,m) & 5284 + surf%rad_sw_in(m) & 5285 * ( exp( -win_absorp * surf%zw_window(kw-1,m) ) & 5286 - exp(-win_absorp * surf%zw_window(kw,m) ) & 5287 ) & 5288 ) * surf%ddz_window_stag(kw,m) 5287 5289 5288 5290 ENDDO 5289 5291 5290 wintend(nzt_wall) = ( 1.0_wp / surf%rho_c_window(nzt_wall,m) ) &5291 * ( -surf%lambda_h_window(nzt_wall-1,m)&5292 * ( t_window%val(nzt_wall,m) - t_window%val(nzt_wall-1,m) ) &5293 * surf%ddz_window(nzt_wall,m) &5294 + surf%iwghf_eb_window(m) &5295 + surf%rad_sw_in(m) &5296 * ( exp( -win_absorp &5297 * surf%zw_window(nzt_wall-1,m) )&5298 - exp( -win_absorp&5299 * surf%zw_window(nzt_wall,m) )&5300 ) &5301 ) * surf%ddz_window_stag(nzt_wall,m)5302 5303 t_window_p%val(nzb_wall:nzt_wall,m) = t_window%val(nzb_wall:nzt_wall,m) + dt_3d &5304 * ( tsc(2) * wintend(nzb_wall:nzt_wall) + tsc(3) &5292 wintend(nzt_wall) = ( 1.0_wp / surf%rho_c_window(nzt_wall,m) ) & 5293 * ( -surf%lambda_h_window(nzt_wall-1,m) & 5294 * ( t_window%val(nzt_wall,m) - t_window%val(nzt_wall-1,m) ) & 5295 * surf%ddz_window(nzt_wall,m) & 5296 + surf%iwghf_eb_window(m) & 5297 + surf%rad_sw_in(m) & 5298 * ( exp( -win_absorp & 5299 * surf%zw_window(nzt_wall-1,m) ) & 5300 - exp( -win_absorp & 5301 * surf%zw_window(nzt_wall,m) ) & 5302 ) & 5303 ) * surf%ddz_window_stag(nzt_wall,m) 5304 5305 t_window_p%val(nzb_wall:nzt_wall,m) = t_window%val(nzb_wall:nzt_wall,m) + dt_3d & 5306 * ( tsc(2) * wintend(nzb_wall:nzt_wall) + tsc(3) & 5305 5307 * surf%tt_window_m(nzb_wall:nzt_wall,m) ) 5306 5308 5307 5309 ENDIF 5308 5309 5310 ! 5310 5311 !-- Calculate t_wall tendencies for the next Runge-Kutta step … … 5316 5317 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 5317 5318 DO kw = nzb_wall, nzt_wall 5318 surf%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) + &5319 surf%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) + & 5319 5320 5.3125_wp * surf%tt_wall_m(kw,m) 5320 5321 ENDDO … … 5332 5333 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 5333 5334 DO kw = nzb_wall, nzt_wall 5334 surf%tt_window_m(kw,m) = -9.5625_wp * wintend(kw) + &5335 surf%tt_window_m(kw,m) = -9.5625_wp * wintend(kw) + & 5335 5336 5.3125_wp * surf%tt_window_m(kw,m) 5336 5337 ENDDO … … 5339 5340 ENDIF 5340 5341 ENDDO 5341 !$OMP END PARALLEL5342 5342 5343 5343 IF ( debug_output_timestep ) THEN … … 5403 5403 ENDIF 5404 5404 5405 !$OMP PARALLEL PRIVATE (m, i, j, k, kw, lambda_h_green_sat, ke, lambda_green_temp, gtend, & 5406 !$OMP& tend, h_vg, gamma_green_temp, m_total, root_extr_green) 5407 !$OMP DO SCHEDULE (STATIC) 5405 !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw, lambda_h_green_sat, ke, lambda_green_temp, & 5406 !$OMP& gtend, tend, h_vg, gamma_green_temp, m_total, root_extr_green) SCHEDULE (STATIC) 5408 5407 DO m = 1, surf%ns 5409 5408 IF (surf%frac(m,ind_pav_green) > 0.0_wp) THEN … … 5417 5416 ! 5418 5417 !-- Calculate volumetric heat capacity of the soil, taking into account water content 5419 surf%rho_c_total_green(kw,m) = (surf%rho_c_green(kw,m) &5420 * (1.0_wp - swc_sat_h(l)%val(kw,m)) 5418 surf%rho_c_total_green(kw,m) = (surf%rho_c_green(kw,m) & 5419 * (1.0_wp - swc_sat_h(l)%val(kw,m)) & 5421 5420 + rho_c_water * swc_h(l)%val(kw,m)) 5422 5421 5423 5422 ! 5424 5423 !-- Calculate soil heat conductivity at the center of the soil layers 5425 lambda_h_green_sat = lambda_h_green_sm ** ( 1.0_wp - swc_sat_h(l)%val(kw,m) ) 5424 lambda_h_green_sat = lambda_h_green_sm ** ( 1.0_wp - swc_sat_h(l)%val(kw,m) ) & 5426 5425 * lambda_h_water ** swc_h(l)%val(kw,m) 5427 5426 5428 5427 ke = 1.0_wp + LOG10( MAX( 0.1_wp,swc_h(l)%val(kw,m) / swc_sat_h(l)%val(kw,m) ) ) 5429 5428 5430 lambda_green_temp(kw) = ke * (lambda_h_green_sat - lambda_h_green_dry) 5429 lambda_green_temp(kw) = ke * (lambda_h_green_sat - lambda_h_green_dry) & 5431 5430 + lambda_h_green_dry 5432 5431 … … 5439 5438 !-- For pavement surface, the true pavement depth is considered 5440 5439 DO kw = nzb_wall, nzt_wall 5441 surf%lambda_h_green(kw,m) = ( lambda_green_temp(kw+1) + lambda_green_temp(kw) ) &5440 surf%lambda_h_green(kw,m) = ( lambda_green_temp(kw+1) + lambda_green_temp(kw) ) & 5442 5441 * 0.5_wp 5443 5442 ENDDO … … 5447 5446 !-- Prognostic equation for ground/roof temperature t_green_h 5448 5447 gtend(:) = 0.0_wp 5449 gtend(nzb_wall) = ( 1.0_wp / surf%rho_c_total_green(nzb_wall,m) ) &5450 * ( surf%lambda_h_green(nzb_wall,m) &5451 * ( t_green_h(l)%val(nzb_wall+1,m) 5452 - t_green_h(l)%val(nzb_wall,m) ) 5453 * surf%ddz_green(nzb_wall+1,m) &5454 + surf%wghf_eb_green(m) &5448 gtend(nzb_wall) = ( 1.0_wp / surf%rho_c_total_green(nzb_wall,m) ) & 5449 * ( surf%lambda_h_green(nzb_wall,m) & 5450 * ( t_green_h(l)%val(nzb_wall+1,m) & 5451 - t_green_h(l)%val(nzb_wall,m) ) & 5452 * surf%ddz_green(nzb_wall+1,m) & 5453 + surf%wghf_eb_green(m) & 5455 5454 ) * surf%ddz_green_stag(nzb_wall,m) 5456 5455 5457 5456 DO kw = nzb_wall+1, nzt_wall 5458 gtend(kw) = ( 1.0_wp / surf%rho_c_total_green(kw,m) ) &5459 * ( surf%lambda_h_green(kw,m) &5460 * ( t_green_h(l)%val(kw+1,m) - t_green_h(l)%val(kw,m) ) 5461 * surf%ddz_green(kw+1,m) &5462 - surf%lambda_h_green(kw-1,m) &5463 * ( t_green_h(l)%val(kw,m) - t_green_h(l)%val(kw-1,m) ) 5464 * surf%ddz_green(kw,m) &5457 gtend(kw) = ( 1.0_wp / surf%rho_c_total_green(kw,m) ) & 5458 * ( surf%lambda_h_green(kw,m) & 5459 * ( t_green_h(l)%val(kw+1,m) - t_green_h(l)%val(kw,m) ) & 5460 * surf%ddz_green(kw+1,m) & 5461 - surf%lambda_h_green(kw-1,m) & 5462 * ( t_green_h(l)%val(kw,m) - t_green_h(l)%val(kw-1,m) ) & 5463 * surf%ddz_green(kw,m) & 5465 5464 ) * surf%ddz_green_stag(kw,m) 5466 5465 ENDDO 5467 5466 5468 t_green_h_p(l)%val(nzb_wall:nzt_wall,m) = t_green_h(l)%val(nzb_wall:nzt_wall,m) &5469 + dt_3d * ( tsc(2) * gtend(nzb_wall:nzt_wall) + tsc(3) &5467 t_green_h_p(l)%val(nzb_wall:nzt_wall,m) = t_green_h(l)%val(nzb_wall:nzt_wall,m) & 5468 + dt_3d * ( tsc(2) * gtend(nzb_wall:nzt_wall) + tsc(3) & 5470 5469 * surf%tt_green_m(nzb_wall:nzt_wall,m) ) 5471 5470 … … 5480 5479 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 5481 5480 DO kw = nzb_wall, nzt_wall 5482 surf%tt_green_m(kw,m) = -9.5625_wp * gtend(kw) + 5.3125_wp &5481 surf%tt_green_m(kw,m) = -9.5625_wp * gtend(kw) + 5.3125_wp & 5483 5482 * surf%tt_green_m(kw,m) 5484 5483 ENDDO … … 5490 5489 ! 5491 5490 !-- Calculate soil diffusivity at the center of the soil layers 5492 lambda_green_temp(kw) = ( - b_ch * surf%gamma_w_green_sat(kw,m) * psi_sat &5493 / swc_sat_h(l)%val(kw,m) ) &5494 * ( MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) &5491 lambda_green_temp(kw) = ( - b_ch * surf%gamma_w_green_sat(kw,m) * psi_sat & 5492 / swc_sat_h(l)%val(kw,m) ) & 5493 * ( MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) & 5495 5494 / swc_sat_h(l)%val(kw,m) )**( b_ch + 2.0_wp ) 5496 5495 … … 5500 5499 ! 5501 5500 !-- Calculate the hydraulic conductivity after Van Genuchten (1980) 5502 h_vg = ( ( (swc_res_h(l)%val(kw,m) - swc_sat_h(l)%val(kw,m)) &5503 / ( swc_res_h(l)%val(kw,m) - &5504 MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) ) )** 5505 ( surf%n_vg_green(m) / (surf%n_vg_green(m) - 1.0_wp ) ) &5506 - 1.0_wp 5501 h_vg = ( ( (swc_res_h(l)%val(kw,m) - swc_sat_h(l)%val(kw,m)) & 5502 / ( swc_res_h(l)%val(kw,m) - & 5503 MAX( swc_h(l)%val(kw,m), wilt_h(l)%val(kw,m) ) ) )** & 5504 ( surf%n_vg_green(m) / (surf%n_vg_green(m) - 1.0_wp ) ) & 5505 - 1.0_wp & 5507 5506 )** ( 1.0_wp / surf%n_vg_green(m) ) / surf%alpha_vg_green(m) 5508 5507 5509 5508 5510 gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) &5511 * ( ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )** &5512 surf%n_vg_green(m) )** &5513 ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) ) &5514 - ( surf%alpha_vg_green(m) * h_vg )** &5515 ( surf%n_vg_green(m) - 1.0_wp) )**2 &5516 ) / ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )** &5517 surf%n_vg_green(m) )** &5518 ( ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) ) &5519 *( surf%l_vg_green(m) + 2.0_wp) ) &5509 gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) & 5510 * ( ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )** & 5511 surf%n_vg_green(m) )** & 5512 ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) ) & 5513 - ( surf%alpha_vg_green(m) * h_vg )** & 5514 ( surf%n_vg_green(m) - 1.0_wp) )**2 & 5515 ) / ( ( 1.0_wp + ( surf%alpha_vg_green(m) * h_vg )** & 5516 surf%n_vg_green(m) )** & 5517 ( ( 1.0_wp - 1.0_wp / surf%n_vg_green(m) ) & 5518 *( surf%l_vg_green(m) + 2.0_wp) ) & 5520 5519 ) 5521 5520 … … 5523 5522 !-- Parametrization of Clapp & Hornberger 5524 5523 ELSE 5525 gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) * ( swc_h(l)%val(kw,m) &5524 gamma_green_temp(kw) = surf%gamma_w_green_sat(kw,m) * ( swc_h(l)%val(kw,m) & 5526 5525 / swc_sat_h(l)%val(kw,m) )**( 2.0_wp * b_ch + 3.0_wp ) 5527 5526 ENDIF … … 5538 5537 DO kw = nzb_wall, nzt_wall-1 5539 5538 5540 surf%lambda_w_green(kw,m) = ( lambda_green_temp(kw+1) &5541 + lambda_green_temp(kw) ) 5539 surf%lambda_w_green(kw,m) = ( lambda_green_temp(kw+1) & 5540 + lambda_green_temp(kw) ) & 5542 5541 * 0.5_wp 5543 surf%gamma_w_green(kw,m) = ( gamma_green_temp(kw+1) &5544 + gamma_green_temp(kw) ) 5542 surf%gamma_w_green(kw,m) = ( gamma_green_temp(kw+1) & 5543 + gamma_green_temp(kw) ) & 5545 5544 * 0.5_wp 5546 5545 … … 5588 5587 tend(:) = 0.0_wp 5589 5588 5590 tend(nzb_wall) = ( surf_usm_h(l)%lambda_w_green(nzb_wall,m) 5591 * ( swc_h(l)%val(nzb_wall+1,m) - swc_h(l)%val(nzb_wall,m) ) 5592 * surf_usm_h(l)%ddz_green(nzb_wall+1,m) 5593 - surf_usm_h(l)%gamma_w_green(nzb_wall,m) 5594 - ( root_extr_green(nzb_wall) * surf_usm_h(l)%qsws_veg(m) 5595 ! + surf_usm_h(l)%qsws_soil_green(m) 5596 ) * drho_l_lv ) 5589 tend(nzb_wall) = ( surf_usm_h(l)%lambda_w_green(nzb_wall,m) & 5590 * ( swc_h(l)%val(nzb_wall+1,m) - swc_h(l)%val(nzb_wall,m) ) & 5591 * surf_usm_h(l)%ddz_green(nzb_wall+1,m) & 5592 - surf_usm_h(l)%gamma_w_green(nzb_wall,m) & 5593 - ( root_extr_green(nzb_wall) * surf_usm_h(l)%qsws_veg(m) & 5594 ! + surf_usm_h(l)%qsws_soil_green(m) & 5595 ) * drho_l_lv ) & 5597 5596 * surf_usm_h(l)%ddz_green_stag(nzb_wall,m) 5598 5597 5599 5598 DO kw = nzb_wall+1, nzt_wall-1 5600 tend(kw) = ( surf_usm_h(l)%lambda_w_green(kw,m) 5601 * ( swc_h(l)%val(kw+1,m) - swc_h(l)%val(kw,m) ) 5602 * surf_usm_h(l)%ddz_green(kw+1,m) 5603 - surf_usm_h(l)%gamma_w_green(kw,m) 5604 - surf_usm_h(l)%lambda_w_green(kw-1,m) 5605 * ( swc_h(l)%val(kw,m) - swc_h(l)%val(kw-1,m) ) 5606 * surf_usm_h(l)%ddz_green(kw,m) 5607 + surf_usm_h(l)%gamma_w_green(kw-1,m) 5608 - (root_extr_green(kw) 5609 * surf_usm_h(l)%qsws_veg(m) 5610 * drho_l_lv) 5599 tend(kw) = ( surf_usm_h(l)%lambda_w_green(kw,m) & 5600 * ( swc_h(l)%val(kw+1,m) - swc_h(l)%val(kw,m) ) & 5601 * surf_usm_h(l)%ddz_green(kw+1,m) & 5602 - surf_usm_h(l)%gamma_w_green(kw,m) & 5603 - surf_usm_h(l)%lambda_w_green(kw-1,m) & 5604 * ( swc_h(l)%val(kw,m) - swc_h(l)%val(kw-1,m) ) & 5605 * surf_usm_h(l)%ddz_green(kw,m) & 5606 + surf_usm_h(l)%gamma_w_green(kw-1,m) & 5607 - (root_extr_green(kw) & 5608 * surf_usm_h(l)%qsws_veg(m) & 5609 * drho_l_lv) & 5611 5610 ) * surf_usm_h(l)%ddz_green_stag(kw,m) 5612 5611 5613 5612 ENDDO 5614 tend(nzt_wall) = ( - surf_usm_h(l)%gamma_w_green(nzt_wall,m) 5615 - surf_usm_h(l)%lambda_w_green(nzt_wall-1,m) 5616 * (swc_h(l)%val(nzt_wall,m) 5617 - swc_h(l)%val(nzt_wall-1,m)) 5618 * surf_usm_h(l)%ddz_green(nzt_wall,m) 5619 + surf_usm_h(l)%gamma_w_green(nzt_wall-1,m) 5620 - ( root_extr_green(nzt_wall) 5621 * surf_usm_h(l)%qsws_veg(m) 5622 * drho_l_lv ) 5613 tend(nzt_wall) = ( - surf_usm_h(l)%gamma_w_green(nzt_wall,m) & 5614 - surf_usm_h(l)%lambda_w_green(nzt_wall-1,m) & 5615 * (swc_h(l)%val(nzt_wall,m) & 5616 - swc_h(l)%val(nzt_wall-1,m)) & 5617 * surf_usm_h(l)%ddz_green(nzt_wall,m) & 5618 + surf_usm_h(l)%gamma_w_green(nzt_wall-1,m) & 5619 - ( root_extr_green(nzt_wall) & 5620 * surf_usm_h(l)%qsws_veg(m) & 5621 * drho_l_lv ) & 5623 5622 ) * surf_usm_h(l)%ddz_green_stag(nzt_wall,m) 5624 5623 5625 swc_h_p(l)%val(nzb_wall:nzt_wall,m) = swc_h(l)%val(nzb_wall:nzt_wall,m) + dt_3d 5626 * ( tsc(2) * tend(:) + tsc(3) 5627 * surf_usm_h(l)%tswc_h_m(:,m) 5624 swc_h_p(l)%val(nzb_wall:nzt_wall,m) = swc_h(l)%val(nzb_wall:nzt_wall,m) + dt_3d & 5625 * ( tsc(2) * tend(:) + tsc(3) & 5626 * surf_usm_h(l)%tswc_h_m(:,m) & 5628 5627 ) 5629 5628 … … 5643 5642 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 5644 5643 DO kw = nzb_wall, nzt_wall 5645 surf_usm_h(l)%tswc_h_m(kw,m) = -9.5625_wp * tend(kw) + 5.3125_wp 5644 surf_usm_h(l)%tswc_h_m(kw,m) = -9.5625_wp * tend(kw) + 5.3125_wp & 5646 5645 * surf_usm_h(l)%tswc_h_m(kw,m) 5647 5646 ENDDO … … 5652 5651 ENDIF 5653 5652 ENDDO 5654 !$OMP END PARALLEL5655 5653 ELSE 5656 5654 IF ( horizontal) THEN … … 5665 5663 t_green => t_green_v(l) 5666 5664 ENDIF 5667 !$OMP DOSCHEDULE (STATIC)5665 !$OMP PARALLEL DO PRIVATE (m, i, j, k, kw) SCHEDULE (STATIC) 5668 5666 DO m = 1, surf%ns 5669 5667 IF (surf%frac(m,ind_pav_green) > 0.0_wp) THEN … … 5684 5682 ! !-- Prognostic equation for green temperature t_green_v 5685 5683 ! gtend(:) = 0.0_wp 5686 ! gtend(nzb_wall) = (1.0_wp / surf%rho_c_green(nzb_wall,m)) * &5687 ! ( surf%lambda_h_green(nzb_wall,m) * &5688 ! ( t_green%val(nzb_wall+1,m) &5689 ! - t_green%val(nzb_wall,m) ) * &5690 ! surf%ddz_green(nzb_wall+1,m) &5691 ! + surf%wghf_eb(m) ) * &5684 ! gtend(nzb_wall) = (1.0_wp / surf%rho_c_green(nzb_wall,m)) * & 5685 ! ( surf%lambda_h_green(nzb_wall,m) * & 5686 ! ( t_green%val(nzb_wall+1,m) & 5687 ! - t_green%val(nzb_wall,m) ) * & 5688 ! surf%ddz_green(nzb_wall+1,m) & 5689 ! + surf%wghf_eb(m) ) * & 5692 5690 ! surf%ddz_green_stag(nzb_wall,m) 5693 5691 ! 5694 5692 ! DO kw = nzb_wall+1, nzt_wall 5695 ! gtend(kw) = (1.0_wp / surf%rho_c_green(kw,m)) &5696 ! * ( surf%lambda_h_green(kw,m) &5697 ! * ( t_green%val(kw+1,m) - t_green%val(kw,m) ) &5698 ! * surf%ddz_green(kw+1,m) &5699 ! - surf%lambda_h(kw-1,m) &5700 ! * ( t_green%val(kw,m) - t_green%val(kw-1,m) ) &5701 ! * surf%ddz_green(kw,m) ) &5693 ! gtend(kw) = (1.0_wp / surf%rho_c_green(kw,m)) & 5694 ! * ( surf%lambda_h_green(kw,m) & 5695 ! * ( t_green%val(kw+1,m) - t_green%val(kw,m) ) & 5696 ! * surf%ddz_green(kw+1,m) & 5697 ! - surf%lambda_h(kw-1,m) & 5698 ! * ( t_green%val(kw,m) - t_green%val(kw-1,m) ) & 5699 ! * surf%ddz_green(kw,m) ) & 5702 5700 ! * surf%ddz_green_stag(kw,m) 5703 5701 ! ENDDO 5704 5702 ! 5705 ! t_green_v_p(l)%val(nzb_wall:nzt_wall,m) = &5706 ! t_green%val(nzb_wall:nzt_wall,m) &5707 ! + dt_3d * ( tsc(2) &5708 ! * gtend(nzb_wall:nzt_wall) + tsc(3) &5703 ! t_green_v_p(l)%val(nzb_wall:nzt_wall,m) = & 5704 ! t_green%val(nzb_wall:nzt_wall,m) & 5705 ! + dt_3d * ( tsc(2) & 5706 ! * gtend(nzb_wall:nzt_wall) + tsc(3) & 5709 5707 ! * surf%tt_green_m(nzb_wall:nzt_wall,m) ) 5710 5708 ! … … 5716 5714 ! surf%tt_green_m(kw,m) = gtend(kw) 5717 5715 ! ENDDO 5718 ! ELSEIF ( intermediate_timestep_count < &5716 ! ELSEIF ( intermediate_timestep_count < & 5719 5717 ! intermediate_timestep_count_max ) THEN 5720 5718 ! DO kw = nzb_wall, nzt_wall 5721 ! surf%tt_green_m(kw,m) = &5722 ! - 9.5625_wp * gtend(kw) + &5719 ! surf%tt_green_m(kw,m) = & 5720 ! - 9.5625_wp * gtend(kw) + & 5723 5721 ! 5.3125_wp * surf%tt_green_m(kw,m) 5724 5722 ! ENDDO … … 6768 6766 ! 6769 6767 !-- First, treat horizontal surface elements 6770 !$OMP PARALLEL PRIVATE (m, i, j, k, frac_win, frac_wall, frac_green, lambda_surface, ueff, & 6771 !$OMP& lambda_surface_window, lambda_surface_green, qv1, rho_cp, rho_lv, & 6772 !$OMP& drho_l_lv, f_shf, f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, & 6773 !$OMP& f3, f_qsws_veg, q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1, & 6774 !$OMP& coef_window_1, coef_green_1, coef_2, coef_window_2, coef_green_2, & 6775 !$OMP& stend_wall, stend_window, stend_green, tend, m_liq_max) 6776 !$OMP DO SCHEDULE (STATIC) 6768 !$OMP PARALLEL DO PRIVATE (m, i, j, k, frac_win, frac_wall, frac_green, lambda_surface, & 6769 !$OMP& lambda_surface_window, lambda_surface_green, ueff, qv1, rho_cp, rho_lv, & 6770 !$OMP& drho_l_lv, f_shf, f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, & 6771 !$OMP& f3, f_qsws_veg, q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1, & 6772 !$OMP& coef_window_1, coef_green_1, coef_2, coef_window_2, coef_green_2, & 6773 !$OMP& stend_wall, stend_window, stend_green, tend, m_liq_max) SCHEDULE (STATIC) 6777 6774 DO m = 1, surf%ns 6778 6775 ! … … 7200 7197 7201 7198 ENDDO 7202 !$OMP END PARALLEL7203 7199 7204 7200 !
Note: See TracChangeset
for help on using the changeset viewer.
