Changeset 1484 for palm/trunk/SOURCE/plant_canopy_model.f90
- Timestamp:
- Oct 21, 2014 10:53:05 AM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/plant_canopy_model.f90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/plant_canopy_model.f90
r1341 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! module plant_canopy_model_mod now contains a subroutine for the 24 ! initialization of the canopy model (init_plant_canopy), 25 ! limitation of the canopy drag (previously accounted for by calculation of 26 ! a limiting canopy timestep for the determination of the maximum LES timestep 27 ! in subroutine timestep) is now realized by the calculation of pre-tendencies 28 ! and preliminary velocities in subroutine plant_canopy_model, 29 ! some redundant MPI communication removed in subroutine init_plant_canopy 30 ! (was previously in init_3d_model), 31 ! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the 32 ! respective grid is now provided only by lad_s (e.g. in the calculation of 33 ! the tendency terms or of cum_lai_hf), 34 ! drag_coefficient, lai, leaf_surface_concentration, 35 ! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff, 36 ! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc, 37 ! respectively, 38 ! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants, 39 ! USE-statements and ONLY-lists modified accordingly 23 40 ! 24 41 ! Former revisions: … … 46 63 ! Description: 47 64 ! ------------ 48 ! Evaluation of sinks and sources of momentum, heat and scalar concentration 49 ! due to canopy elements 65 ! 1) Initialization of the canopy model, e.g. construction of leaf area density 66 ! profile (subroutine init_plant_canopy). 67 ! 2) Calculation of sinks and sources of momentum, heat and scalar concentration 68 ! due to canopy elements (subroutine plant_canopy_model). 50 69 !------------------------------------------------------------------------------! 70 USE arrays_3d, & 71 ONLY: dzu, dzw, e, q, shf, tend, u, v, w, zu, zw 72 73 USE indices, & 74 ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, & 75 nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt 76 77 USE kinds 78 79 80 IMPLICIT NONE 81 82 83 CHARACTER (LEN=20) :: canopy_mode = 'block' !: canopy coverage 84 85 INTEGER(iwp) :: pch_index = 0 !: plant canopy height/top index 86 INTEGER(iwp) :: & 87 lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index) 88 89 LOGICAL :: calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func. 90 LOGICAL :: plant_canopy = .FALSE. !: switch for use of canopy model 91 92 REAL(wp) :: alpha_lad = 9999999.9_wp !: coefficient for lad calculation 93 REAL(wp) :: beta_lad = 9999999.9_wp !: coefficient for lad calculation 94 REAL(wp) :: canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter) 95 REAL(wp) :: cdc = 0.0_wp !: canopy drag coeff. (abbreviation used in equations) 96 REAL(wp) :: cthf = 0.0_wp !: canopy top heat flux 97 REAL(wp) :: dt_plant_canopy = 0.0_wp !: timestep account. for canopy drag 98 REAL(wp) :: lad_surface = 0.0_wp !: lad surface value 99 REAL(wp) :: lai_beta = 0.0_wp !: leaf area index (lai) for lad calc. 100 REAL(wp) :: & 101 leaf_scalar_exch_coeff = 0.0_wp !: canopy scalar exchange coeff. 102 REAL(wp) :: & 103 leaf_surface_conc = 0.0_wp !: leaf surface concentration 104 REAL(wp) :: lsec = 0.0_wp !: leaf scalar exchange coeff. 105 REAL(wp) :: lsc = 0.0_wp !: leaf surface concentration 106 107 REAL(wp) :: & 108 lad_vertical_gradient(10) = 0.0_wp !: lad gradient 109 REAL(wp) :: & 110 lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m) 111 112 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !: leaf area density 113 REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !: preliminary lad 114 115 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 116 canopy_heat_flux !: canopy heat flux 117 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !: cumulative lai for heatflux calc. 118 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !: lad on scalar-grid 119 120 121 SAVE 122 51 123 52 124 PRIVATE 53 PUBLIC plant_canopy_model 125 PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, & 126 canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad, & 127 lad_s, lad_surface, lad_vertical_gradient, & 128 lad_vertical_gradient_level, lad_vertical_gradient_level_ind, & 129 lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec, & 130 pch_index, plant_canopy, plant_canopy_model 131 132 133 INTERFACE init_plant_canopy 134 MODULE PROCEDURE init_plant_canopy 135 END INTERFACE init_plant_canopy 54 136 55 137 INTERFACE plant_canopy_model … … 58 140 END INTERFACE plant_canopy_model 59 141 142 143 144 60 145 CONTAINS 61 146 62 147 63 148 !------------------------------------------------------------------------------! 64 ! Call for all grid points 149 ! Description: 150 ! ------------ 151 !-- Initialization of the plant canopy model 152 !------------------------------------------------------------------------------! 153 SUBROUTINE init_plant_canopy 154 155 156 USE control_parameters, & 157 ONLY: dz, ocean, passive_scalar 158 159 160 IMPLICIT NONE 161 162 INTEGER(iwp) :: i !: running index 163 INTEGER(iwp) :: j !: running index 164 INTEGER(iwp) :: k !: running index 165 166 REAL(wp) :: int_bpdf !: vertical integral for lad-profile construction 167 REAL(wp) :: dzh !: vertical grid spacing in units of canopy height 168 REAL(wp) :: gradient !: gradient for lad-profile construction 169 REAL(wp) :: canopy_height !: canopy height for lad-profile construction 170 171 ! 172 !-- Allocate one-dimensional arrays for the computation of the 173 !-- leaf area density (lad) profile 174 ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) ) 175 lad = 0.0_wp 176 pre_lad = 0.0_wp 177 178 ! 179 !-- Compute the profile of leaf area density used in the plant 180 !-- canopy model. The profile can either be constructed from 181 !-- prescribed vertical gradients of the leaf area density or by 182 !-- using a beta probability density function (see e.g. Markkanen et al., 183 !-- 2003: Boundary-Layer Meteorology, 106, 437-459) 184 IF ( .NOT. calc_beta_lad_profile ) THEN 185 186 ! 187 !-- Use vertical gradients for lad-profile construction 188 i = 1 189 gradient = 0.0_wp 190 191 IF ( .NOT. ocean ) THEN 192 193 lad(0) = lad_surface 194 lad_vertical_gradient_level_ind(1) = 0 195 196 DO k = 1, pch_index 197 IF ( i < 11 ) THEN 198 IF ( lad_vertical_gradient_level(i) < zu(k) .AND. & 199 lad_vertical_gradient_level(i) >= 0.0_wp ) THEN 200 gradient = lad_vertical_gradient(i) 201 lad_vertical_gradient_level_ind(i) = k - 1 202 i = i + 1 203 ENDIF 204 ENDIF 205 IF ( gradient /= 0.0_wp ) THEN 206 IF ( k /= 1 ) THEN 207 lad(k) = lad(k-1) + dzu(k) * gradient 208 ELSE 209 lad(k) = lad_surface + dzu(k) * gradient 210 ENDIF 211 ELSE 212 lad(k) = lad(k-1) 213 ENDIF 214 ENDDO 215 216 ENDIF 217 218 ! 219 !-- In case of no given leaf area density gradients, choose a vanishing 220 !-- gradient. This information is used for the HEADER and the RUN_CONTROL 221 !-- file. 222 IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN 223 lad_vertical_gradient_level(1) = 0.0_wp 224 ENDIF 225 226 ELSE 227 228 ! 229 !-- Use beta function for lad-profile construction 230 int_bpdf = 0.0_wp 231 canopy_height = pch_index * dz 232 233 DO k = nzb, pch_index 234 int_bpdf = int_bpdf + & 235 ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * & 236 ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) * & 237 ( ( zw(k+1)-zw(k) ) / canopy_height ) ) 238 ENDDO 239 240 ! 241 !-- Preliminary lad profile (defined on w-grid) 242 DO k = nzb, pch_index 243 pre_lad(k) = lai_beta * & 244 ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * & 245 ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) / & 246 int_bpdf & 247 ) / canopy_height 248 ENDDO 249 250 ! 251 !-- Final lad profile (defined on scalar-grid level, since most prognostic 252 !-- quantities are defined there, hence, less interpolation is required 253 !-- when calculating the canopy tendencies) 254 lad(0) = pre_lad(0) 255 DO k = nzb+1, pch_index 256 lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) ) 257 ENDDO 258 259 ENDIF 260 261 ! 262 !-- Allocate 3D-array for the leaf area density (lad_s). In case of a 263 !-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for 264 !-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux. 265 ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 266 267 IF ( cthf /= 0.0_wp ) THEN 268 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 269 canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 270 ENDIF 271 272 ! 273 !-- Initialize canopy parameters cdc (canopy drag coefficient), 274 !-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration) 275 !-- with the prescribed values 276 cdc = canopy_drag_coeff 277 lsec = leaf_scalar_exch_coeff 278 lsc = leaf_surface_conc 279 280 ! 281 !-- Initialization of the canopy coverage in the model domain: 282 !-- Setting the parameter canopy_mode = 'block' initializes a canopy, which 283 !-- fully covers the domain surface 284 SELECT CASE ( TRIM( canopy_mode ) ) 285 286 CASE( 'block' ) 287 288 DO i = nxlg, nxrg 289 DO j = nysg, nyng 290 lad_s(:,j,i) = lad(:) 291 ENDDO 292 ENDDO 293 294 CASE DEFAULT 295 296 ! 297 !-- The DEFAULT case is reached either if the parameter 298 !-- canopy mode contains a wrong character string or if the 299 !-- user has coded a special case in the user interface. 300 !-- There, the subroutine user_init_plant_canopy checks 301 !-- which of these two conditions applies. 302 CALL user_init_plant_canopy 303 304 END SELECT 305 306 ! 307 !-- Initialization of the canopy heat source distribution 308 IF ( cthf /= 0.0_wp ) THEN 309 ! 310 !-- Piecewise calculation of the leaf area index by vertical 311 !-- integration of the leaf area density 312 cum_lai_hf(:,:,:) = 0.0_wp 313 DO i = nxlg, nxrg 314 DO j = nysg, nyng 315 DO k = pch_index-1, 0, -1 316 IF ( k == pch_index-1 ) THEN 317 cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + & 318 ( 0.5_wp * lad_s(k+1,j,i) * & 319 ( zw(k+1) - zu(k+1) ) ) + & 320 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + & 321 lad_s(k,j,i) ) + & 322 lad_s(k+1,j,i) ) * & 323 ( zu(k+1) - zw(k) ) ) 324 ELSE 325 cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + & 326 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + & 327 lad_s(k+1,j,i) ) + & 328 lad_s(k+1,j,i) ) * & 329 ( zw(k+1) - zu(k+1) ) ) + & 330 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + & 331 lad_s(k,j,i) ) + & 332 lad_s(k+1,j,i) ) * & 333 ( zu(k+1) - zw(k) ) ) 334 ENDIF 335 ENDDO 336 ENDDO 337 ENDDO 338 339 ! 340 !-- Calculation of the upward kinematic vertical heat flux within the 341 !-- canopy 342 DO i = nxlg, nxrg 343 DO j = nysg, nyng 344 DO k = 0, pch_index 345 canopy_heat_flux(k,j,i) = cthf * & 346 exp( -0.6_wp * cum_lai_hf(k,j,i) ) 347 ENDDO 348 ENDDO 349 ENDDO 350 351 ! 352 !-- The surface heat flux is set to the surface value of the calculated 353 !-- in-canopy heat flux distribution 354 shf(:,:) = canopy_heat_flux(0,:,:) 355 356 ENDIF 357 358 359 360 END SUBROUTINE init_plant_canopy 361 362 363 364 !------------------------------------------------------------------------------! 365 ! Description: 366 ! ------------ 367 !-- Calculation of the tendency terms, accounting for the effect of the plant 368 !-- canopy on momentum and scalar quantities. 369 !-- 370 !-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 371 !-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 372 !-- The lad on the w-grid is vertically interpolated from the surrounding 373 !-- lad_s. The upper boundary of the canopy is defined on the w-grid at 374 !-- k = pch_index. Here, the lad is zero. 375 !-- 376 !-- The canopy drag must be limited (previously accounted for by calculation of 377 !-- a limiting canopy timestep for the determination of the maximum LES timestep 378 !-- in subroutine timestep), since it is physically impossible that the canopy 379 !-- drag alone can locally change the sign of a velocity component. This 380 !-- limitation is realized by calculating preliminary tendencies and velocities. 381 !-- It is subsequently checked if the preliminary new velocity has a different 382 !-- sign than the current velocity. If so, the tendency is limited in a way that 383 !-- the velocity can at maximum be reduced to zero by the canopy drag. 384 !-- 385 !-- 386 !-- Call for all grid points 65 387 !------------------------------------------------------------------------------! 66 388 SUBROUTINE plant_canopy_model( component ) 67 389 68 USE arrays_3d, &69 ONLY: canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w, &70 q, sec, sls, tend, u, v, w71 390 72 391 USE control_parameters, & 73 ONLY: pch_index, message_string 74 75 USE indices, & 76 ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner, & 77 nzb_v_inner, nzb_w_inner 392 ONLY: dt_3d, message_string 78 393 79 394 USE kinds … … 81 396 IMPLICIT NONE 82 397 83 INTEGER(iwp) :: component !: 84 INTEGER(iwp) :: i !: 85 INTEGER(iwp) :: j !: 86 INTEGER(iwp) :: k !: 398 INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e) 399 INTEGER(iwp) :: i !: running index 400 INTEGER(iwp) :: j !: running index 401 INTEGER(iwp) :: k !: running index 402 403 REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d) 404 REAL(wp) :: lad_local !: local lad value 405 REAL(wp) :: pre_tend !: preliminary tendency 406 REAL(wp) :: pre_u !: preliminary u-value 407 REAL(wp) :: pre_v !: preliminary v-value 408 REAL(wp) :: pre_w !: preliminary w-value 409 410 411 ddt_3d = 1.0_wp / dt_3d 87 412 88 413 ! 89 !-- Compute drag for the three velocity components and the SGS-TKE 414 !-- Compute drag for the three velocity components and the SGS-TKE: 90 415 SELECT CASE ( component ) 91 416 … … 96 421 DO j = nys, nyn 97 422 DO k = nzb_u_inner(j,i)+1, pch_index 98 tend(k,j,i) = tend(k,j,i) - & 99 cdc(k,j,i) * lad_u(k,j,i) * & 100 SQRT( u(k,j,i)**2 + & 101 ( ( v(k,j,i-1) + & 102 v(k,j,i) + & 103 v(k,j+1,i) + & 104 v(k,j+1,i-1) ) & 105 / 4.0_wp )**2 + & 106 ( ( w(k-1,j,i-1) + & 107 w(k-1,j,i) + & 108 w(k,j,i-1) + & 109 w(k,j,i) ) & 110 / 4.0_wp )**2 ) * & 111 u(k,j,i) 423 424 ! 425 !-- In order to create sharp boundaries of the plant canopy, 426 !-- the lad on the u-grid at index (k,j,i) is equal to 427 !-- lad_s(k,j,i), rather than being interpolated from the 428 !-- surrounding lad_s, because this would yield smaller lad 429 !-- at the canopy boundaries than inside of the canopy. 430 !-- For the same reason, the lad at the rightmost(i+1)canopy 431 !-- boundary on the u-grid equals lad_s(k,j,i). 432 lad_local = lad_s(k,j,i) 433 IF ( lad_local == 0.0_wp .AND. & 434 lad_s(k,j,i-1) > 0.0_wp ) THEN 435 lad_local = lad_s(k,j,i-1) 436 ENDIF 437 438 pre_tend = 0.0_wp 439 pre_u = 0.0_wp 440 ! 441 !-- Calculate preliminary value (pre_tend) of the tendency 442 pre_tend = - cdc * & 443 lad_local * & 444 SQRT( u(k,j,i)**2 + & 445 ( 0.25_wp * ( v(k,j,i-1) + & 446 v(k,j,i) + & 447 v(k,j+1,i) + & 448 v(k,j+1,i-1) ) & 449 )**2 + & 450 ( 0.25_wp * ( w(k-1,j,i-1) + & 451 w(k-1,j,i) + & 452 w(k,j,i-1) + & 453 w(k,j,i) ) & 454 )**2 & 455 ) * & 456 u(k,j,i) 457 458 ! 459 !-- Calculate preliminary new velocity, based on pre_tend 460 pre_u = u(k,j,i) + dt_3d * pre_tend 461 ! 462 !-- Compare sign of old velocity and new preliminary velocity, 463 !-- and in case the signs are different, limit the tendency 464 IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN 465 pre_tend = - u(k,j,i) * ddt_3d 466 ELSE 467 pre_tend = pre_tend 468 ENDIF 469 ! 470 !-- Calculate final tendency 471 tend(k,j,i) = tend(k,j,i) + pre_tend 472 112 473 ENDDO 113 474 ENDDO … … 120 481 DO j = nysv, nyn 121 482 DO k = nzb_v_inner(j,i)+1, pch_index 122 tend(k,j,i) = tend(k,j,i) - & 123 cdc(k,j,i) * lad_v(k,j,i) * & 124 SQRT( ( ( u(k,j-1,i) + & 125 u(k,j-1,i+1) + & 126 u(k,j,i) + & 127 u(k,j,i+1) ) & 128 / 4.0_wp )**2 + & 129 v(k,j,i)**2 + & 130 ( ( w(k-1,j-1,i) + & 131 w(k-1,j,i) + & 132 w(k,j-1,i) + & 133 w(k,j,i) ) & 134 / 4.0_wp )**2 ) * & 135 v(k,j,i) 483 484 ! 485 !-- In order to create sharp boundaries of the plant canopy, 486 !-- the lad on the v-grid at index (k,j,i) is equal to 487 !-- lad_s(k,j,i), rather than being interpolated from the 488 !-- surrounding lad_s, because this would yield smaller lad 489 !-- at the canopy boundaries than inside of the canopy. 490 !-- For the same reason, the lad at the northmost(j+1) canopy 491 !-- boundary on the v-grid equals lad_s(k,j,i). 492 lad_local = lad_s(k,j,i) 493 IF ( lad_local == 0.0_wp .AND. & 494 lad_s(k,j-1,i) > 0.0_wp ) THEN 495 lad_local = lad_s(k,j-1,i) 496 ENDIF 497 498 pre_tend = 0.0_wp 499 pre_v = 0.0_wp 500 ! 501 !-- Calculate preliminary value (pre_tend) of the tendency 502 pre_tend = - cdc * & 503 lad_local * & 504 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 505 u(k,j-1,i+1) + & 506 u(k,j,i) + & 507 u(k,j,i+1) ) & 508 )**2 + & 509 v(k,j,i)**2 + & 510 ( 0.25_wp * ( w(k-1,j-1,i) + & 511 w(k-1,j,i) + & 512 w(k,j-1,i) + & 513 w(k,j,i) ) & 514 )**2 & 515 ) * & 516 v(k,j,i) 517 518 ! 519 !-- Calculate preliminary new velocity, based on pre_tend 520 pre_v = v(k,j,i) + dt_3d * pre_tend 521 ! 522 !-- Compare sign of old velocity and new preliminary velocity, 523 !-- and in case the signs are different, limit the tendency 524 IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN 525 pre_tend = - v(k,j,i) * ddt_3d 526 ELSE 527 pre_tend = pre_tend 528 ENDIF 529 ! 530 !-- Calculate final tendency 531 tend(k,j,i) = tend(k,j,i) + pre_tend 532 136 533 ENDDO 137 534 ENDDO … … 143 540 DO i = nxl, nxr 144 541 DO j = nys, nyn 145 DO k = nzb_w_inner(j,i)+1, pch_index 146 tend(k,j,i) = tend(k,j,i) - & 147 cdc(k,j,i) * lad_w(k,j,i) * & 148 SQRT( ( ( u(k,j,i) + & 149 u(k,j,i+1) + & 150 u(k+1,j,i) + & 151 u(k+1,j,i+1) ) & 152 / 4.0_wp )**2 + & 153 ( ( v(k,j,i) + & 154 v(k,j+1,i) + & 155 v(k+1,j,i) + & 156 v(k+1,j+1,i) ) & 157 / 4.0_wp )**2 + & 158 w(k,j,i)**2 ) * & 159 w(k,j,i) 542 DO k = nzb_w_inner(j,i)+1, pch_index-1 543 544 pre_tend = 0.0_wp 545 pre_w = 0.0_wp 546 ! 547 !-- Calculate preliminary value (pre_tend) of the tendency 548 pre_tend = - cdc * & 549 (0.5_wp * & 550 ( lad_s(k+1,j,i) + lad_s(k,j,i) )) * & 551 SQRT( ( 0.25_wp * ( u(k,j,i) + & 552 u(k,j,i+1) + & 553 u(k+1,j,i) + & 554 u(k+1,j,i+1) ) & 555 )**2 + & 556 ( 0.25_wp * ( v(k,j,i) + & 557 v(k,j+1,i) + & 558 v(k+1,j,i) + & 559 v(k+1,j+1,i) ) & 560 )**2 + & 561 w(k,j,i)**2 & 562 ) * & 563 w(k,j,i) 564 ! 565 !-- Calculate preliminary new velocity, based on pre_tend 566 pre_w = w(k,j,i) + dt_3d * pre_tend 567 ! 568 !-- Compare sign of old velocity and new preliminary velocity, 569 !-- and in case the signs are different, limit the tendency 570 IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN 571 pre_tend = - w(k,j,i) * ddt_3d 572 ELSE 573 pre_tend = pre_tend 574 ENDIF 575 ! 576 !-- Calculate final tendency 577 tend(k,j,i) = tend(k,j,i) + pre_tend 578 160 579 ENDDO 161 580 ENDDO … … 168 587 DO j = nys, nyn 169 588 DO k = nzb_s_inner(j,i)+1, pch_index 170 tend(k,j,i) = tend(k,j,i) + & 171 ( canopy_heat_flux(k,j,i) - & 172 canopy_heat_flux(k-1,j,i) ) / & 173 dzw(k) 589 tend(k,j,i) = tend(k,j,i) + & 590 ( canopy_heat_flux(k,j,i) - & 591 canopy_heat_flux(k-1,j,i) ) / dzw(k) 174 592 ENDDO 175 593 ENDDO … … 182 600 DO j = nys, nyn 183 601 DO k = nzb_s_inner(j,i)+1, pch_index 184 tend(k,j,i) = tend(k,j,i) - & 185 sec(k,j,i) * lad_s(k,j,i) * & 186 SQRT( ( ( u(k,j,i) + & 187 u(k,j,i+1) ) & 188 / 2.0_wp )**2 + & 189 ( ( v(k,j,i) + & 190 v(k,j+1,i) ) & 191 / 2.0_wp )**2 + & 192 ( ( w(k-1,j,i) + & 193 w(k,j,i) ) & 194 / 2.0_wp )**2 ) * & 195 ( q(k,j,i) - sls(k,j,i) ) 602 tend(k,j,i) = tend(k,j,i) - & 603 lsec * & 604 lad_s(k,j,i) * & 605 SQRT( ( 0.5_wp * ( u(k,j,i) + & 606 u(k,j,i+1) ) & 607 )**2 + & 608 ( 0.5_wp * ( v(k,j,i) + & 609 v(k,j+1,i) ) & 610 )**2 + & 611 ( 0.5_wp * ( w(k-1,j,i) + & 612 w(k,j,i) ) & 613 )**2 & 614 ) * & 615 ( q(k,j,i) - lsc ) 196 616 ENDDO 197 617 ENDDO … … 204 624 DO j = nys, nyn 205 625 DO k = nzb_s_inner(j,i)+1, pch_index 206 tend(k,j,i) = tend(k,j,i) - & 207 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 208 SQRT( ( ( u(k,j,i) + & 209 u(k,j,i+1) ) & 210 / 2.0_wp )**2 + & 211 ( ( v(k,j,i) + & 212 v(k,j+1,i) ) & 213 / 2.0_wp )**2 + & 214 ( ( w(k,j,i) + & 215 w(k+1,j,i) ) & 216 / 2.0_wp )**2 ) * & 217 e(k,j,i) 626 tend(k,j,i) = tend(k,j,i) - & 627 2.0_wp * cdc * & 628 lad_s(k,j,i) * & 629 SQRT( ( 0.5_wp * ( u(k,j,i) + & 630 u(k,j,i+1) ) & 631 )**2 + & 632 ( 0.5_wp * ( v(k,j,i) + & 633 v(k,j+1,i) ) & 634 )**2 + & 635 ( 0.5_wp * ( w(k,j,i) + & 636 w(k+1,j,i) ) & 637 )**2 & 638 ) * & 639 e(k,j,i) 218 640 ENDDO 219 641 ENDDO 220 642 ENDDO 221 643 644 222 645 CASE DEFAULT 223 646 … … 231 654 232 655 !------------------------------------------------------------------------------! 233 ! Call for grid point i,j 656 ! Description: 657 ! ------------ 658 !-- Calculation of the tendency terms, accounting for the effect of the plant 659 !-- canopy on momentum and scalar quantities. 660 !-- 661 !-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 662 !-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 663 !-- The lad on the w-grid is vertically interpolated from the surrounding 664 !-- lad_s. The upper boundary of the canopy is defined on the w-grid at 665 !-- k = pch_index. Here, the lad is zero. 666 !-- 667 !-- The canopy drag must be limited (previously accounted for by calculation of 668 !-- a limiting canopy timestep for the determination of the maximum LES timestep 669 !-- in subroutine timestep), since it is physically impossible that the canopy 670 !-- drag alone can locally change the sign of a velocity component. This 671 !-- limitation is realized by calculating preliminary tendencies and velocities. 672 !-- It is subsequently checked if the preliminary new velocity has a different 673 !-- sign than the current velocity. If so, the tendency is limited in a way that 674 !-- the velocity can at maximum be reduced to zero by the canopy drag. 675 !-- 676 !-- 677 !-- Call for grid point i,j 234 678 !------------------------------------------------------------------------------! 235 679 SUBROUTINE plant_canopy_model_ij( i, j, component ) 236 680 237 USE arrays_3d, &238 ONLY: canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w, &239 q, sec, sls, tend, u, v, w240 681 241 682 USE control_parameters, & 242 ONLY: pch_index, message_string 243 244 USE indices, & 245 ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner, & 246 nzb_v_inner, nzb_w_inner 683 ONLY: dt_3d, message_string 247 684 248 685 USE kinds … … 250 687 IMPLICIT NONE 251 688 252 INTEGER(iwp) :: component !: 253 INTEGER(iwp) :: i !: 254 INTEGER(iwp) :: j !: 255 INTEGER(iwp) :: k !: 256 257 ! 258 !-- Compute drag for the three velocity components 689 INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e) 690 INTEGER(iwp) :: i !: running index 691 INTEGER(iwp) :: j !: running index 692 INTEGER(iwp) :: k !: running index 693 694 REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d) 695 REAL(wp) :: lad_local !: local lad value 696 REAL(wp) :: pre_tend !: preliminary tendency 697 REAL(wp) :: pre_u !: preliminary u-value 698 REAL(wp) :: pre_v !: preliminary v-value 699 REAL(wp) :: pre_w !: preliminary w-value 700 701 702 ddt_3d = 1.0_wp / dt_3d 703 704 ! 705 !-- Compute drag for the three velocity components and the SGS-TKE 259 706 SELECT CASE ( component ) 260 707 261 708 ! 262 709 !-- u-component 263 CASE ( 1 ) 264 DO k = nzb_u_inner(j,i)+1, pch_index 265 tend(k,j,i) = tend(k,j,i) - & 266 cdc(k,j,i) * lad_u(k,j,i) * & 267 SQRT( u(k,j,i)**2 + & 268 ( ( v(k,j,i-1) + & 269 v(k,j,i) + & 270 v(k,j+1,i) + & 271 v(k,j+1,i-1) ) & 272 / 4.0_wp )**2 + & 273 ( ( w(k-1,j,i-1) + & 274 w(k-1,j,i) + & 275 w(k,j,i-1) + & 276 w(k,j,i) ) & 277 / 4.0_wp )**2 ) * & 278 u(k,j,i) 279 ENDDO 710 CASE ( 1 ) 711 DO k = nzb_u_inner(j,i)+1, pch_index 712 713 ! 714 !-- In order to create sharp boundaries of the plant canopy, 715 !-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i), 716 !-- rather than being interpolated from the surrounding lad_s, 717 !-- because this would yield smaller lad at the canopy boundaries 718 !-- than inside of the canopy. 719 !-- For the same reason, the lad at the rightmost(i+1)canopy 720 !-- boundary on the u-grid equals lad_s(k,j,i). 721 lad_local = lad_s(k,j,i) 722 IF ( lad_local == 0.0_wp .AND. & 723 lad_s(k,j,i-1) > 0.0_wp ) THEN 724 lad_local = lad_s(k,j,i-1) 725 ENDIF 726 727 pre_tend = 0.0_wp 728 pre_u = 0.0_wp 729 ! 730 !-- Calculate preliminary value (pre_tend) of the tendency 731 pre_tend = - cdc * & 732 lad_local * & 733 SQRT( u(k,j,i)**2 + & 734 ( 0.25_wp * ( v(k,j,i-1) + & 735 v(k,j,i) + & 736 v(k,j+1,i) + & 737 v(k,j+1,i-1) ) & 738 )**2 + & 739 ( 0.25_wp * ( w(k-1,j,i-1) + & 740 w(k-1,j,i) + & 741 w(k,j,i-1) + & 742 w(k,j,i) ) & 743 )**2 & 744 ) * & 745 u(k,j,i) 746 747 ! 748 !-- Calculate preliminary new velocity, based on pre_tend 749 pre_u = u(k,j,i) + dt_3d * pre_tend 750 ! 751 !-- Compare sign of old velocity and new preliminary velocity, 752 !-- and in case the signs are different, limit the tendency 753 IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN 754 pre_tend = - u(k,j,i) * ddt_3d 755 ELSE 756 pre_tend = pre_tend 757 ENDIF 758 ! 759 !-- Calculate final tendency 760 tend(k,j,i) = tend(k,j,i) + pre_tend 761 ENDDO 762 280 763 281 764 ! 282 765 !-- v-component 283 CASE ( 2 ) 284 DO k = nzb_v_inner(j,i)+1, pch_index 285 tend(k,j,i) = tend(k,j,i) - & 286 cdc(k,j,i) * lad_v(k,j,i) * & 287 SQRT( ( ( u(k,j-1,i) + & 288 u(k,j-1,i+1) + & 289 u(k,j,i) + & 290 u(k,j,i+1) ) & 291 / 4.0_wp )**2 + & 292 v(k,j,i)**2 + & 293 ( ( w(k-1,j-1,i) + & 294 w(k-1,j,i) + & 295 w(k,j-1,i) + & 296 w(k,j,i) ) & 297 / 4.0_wp )**2 ) * & 298 v(k,j,i) 299 ENDDO 766 CASE ( 2 ) 767 DO k = nzb_v_inner(j,i)+1, pch_index 768 769 ! 770 !-- In order to create sharp boundaries of the plant canopy, 771 !-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i), 772 !-- rather than being interpolated from the surrounding lad_s, 773 !-- because this would yield smaller lad at the canopy boundaries 774 !-- than inside of the canopy. 775 !-- For the same reason, the lad at the northmost(j+1)canopy 776 !-- boundary on the v-grid equals lad_s(k,j,i). 777 lad_local = lad_s(k,j,i) 778 IF ( lad_local == 0.0_wp .AND. & 779 lad_s(k,j-1,i) > 0.0_wp ) THEN 780 lad_local = lad_s(k,j-1,i) 781 ENDIF 782 783 pre_tend = 0.0_wp 784 pre_v = 0.0_wp 785 ! 786 !-- Calculate preliminary value (pre_tend) of the tendency 787 pre_tend = - cdc * & 788 lad_local * & 789 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 790 u(k,j-1,i+1) + & 791 u(k,j,i) + & 792 u(k,j,i+1) ) & 793 )**2 + & 794 v(k,j,i)**2 + & 795 ( 0.25_wp * ( w(k-1,j-1,i) + & 796 w(k-1,j,i) + & 797 w(k,j-1,i) + & 798 w(k,j,i) ) & 799 )**2 & 800 ) * & 801 v(k,j,i) 802 803 ! 804 !-- Calculate preliminary new velocity, based on pre_tend 805 pre_v = v(k,j,i) + dt_3d * pre_tend 806 ! 807 !-- Compare sign of old velocity and new preliminary velocity, 808 !-- and in case the signs are different, limit the tendency 809 IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN 810 pre_tend = - v(k,j,i) * ddt_3d 811 ELSE 812 pre_tend = pre_tend 813 ENDIF 814 ! 815 !-- Calculate final tendency 816 tend(k,j,i) = tend(k,j,i) + pre_tend 817 ENDDO 818 300 819 301 820 ! 302 821 !-- w-component 303 CASE ( 3 ) 304 DO k = nzb_w_inner(j,i)+1, pch_index 305 tend(k,j,i) = tend(k,j,i) - & 306 cdc(k,j,i) * lad_w(k,j,i) * & 307 SQRT( ( ( u(k,j,i) + & 308 u(k,j,i+1) + & 309 u(k+1,j,i) + & 310 u(k+1,j,i+1) ) & 311 / 4.0_wp )**2 + & 312 ( ( v(k,j,i) + & 313 v(k,j+1,i) + & 314 v(k+1,j,i) + & 315 v(k+1,j+1,i) ) & 316 / 4.0_wp )**2 + & 317 w(k,j,i)**2 ) * & 318 w(k,j,i) 319 320 ENDDO 822 CASE ( 3 ) 823 DO k = nzb_w_inner(j,i)+1, pch_index-1 824 825 pre_tend = 0.0_wp 826 pre_w = 0.0_wp 827 ! 828 !-- Calculate preliminary value (pre_tend) of the tendency 829 pre_tend = - cdc * & 830 (0.5_wp * & 831 ( lad_s(k+1,j,i) + lad_s(k,j,i) )) * & 832 SQRT( ( 0.25_wp * ( u(k,j,i) + & 833 u(k,j,i+1) + & 834 u(k+1,j,i) + & 835 u(k+1,j,i+1) ) & 836 )**2 + & 837 ( 0.25_wp * ( v(k,j,i) + & 838 v(k,j+1,i) + & 839 v(k+1,j,i) + & 840 v(k+1,j+1,i) ) & 841 )**2 + & 842 w(k,j,i)**2 & 843 ) * & 844 w(k,j,i) 845 ! 846 !-- Calculate preliminary new velocity, based on pre_tend 847 pre_w = w(k,j,i) + dt_3d * pre_tend 848 ! 849 !-- Compare sign of old velocity and new preliminary velocity, 850 !-- and in case the signs are different, limit the tendency 851 IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN 852 pre_tend = - w(k,j,i) * ddt_3d 853 ELSE 854 pre_tend = pre_tend 855 ENDIF 856 ! 857 !-- Calculate final tendency 858 tend(k,j,i) = tend(k,j,i) + pre_tend 859 ENDDO 321 860 322 861 ! … … 324 863 CASE ( 4 ) 325 864 DO k = nzb_s_inner(j,i)+1, pch_index 326 tend(k,j,i) = tend(k,j,i) + & 327 ( canopy_heat_flux(k,j,i) - & 328 canopy_heat_flux(k-1,j,i) ) / & 329 dzw(k) 865 tend(k,j,i) = tend(k,j,i) + & 866 ( canopy_heat_flux(k,j,i) - & 867 canopy_heat_flux(k-1,j,i) ) / dzw(k) 330 868 ENDDO 331 869 … … 335 873 CASE ( 5 ) 336 874 DO k = nzb_s_inner(j,i)+1, pch_index 337 tend(k,j,i) = tend(k,j,i) - & 338 sec(k,j,i) * lad_s(k,j,i) * & 339 SQRT( ( ( u(k,j,i) + & 340 u(k,j,i+1) ) & 341 / 2.0_wp )**2 + & 342 ( ( v(k,j,i) + & 343 v(k,j+1,i) ) & 344 / 2.0_wp )**2 + & 345 ( ( w(k-1,j,i) + & 346 w(k,j,i) ) & 347 / 2.0_wp )**2 ) * & 348 ( q(k,j,i) - sls(k,j,i) ) 875 tend(k,j,i) = tend(k,j,i) - & 876 lsec * & 877 lad_s(k,j,i) * & 878 SQRT( ( 0.5_wp * ( u(k,j,i) + & 879 u(k,j,i+1) ) & 880 )**2 + & 881 ( 0.5_wp * ( v(k,j,i) + & 882 v(k,j+1,i) ) & 883 )**2 + & 884 ( 0.5_wp * ( w(k-1,j,i) + & 885 w(k,j,i) ) & 886 )**2 & 887 ) * & 888 ( q(k,j,i) - lsc ) 349 889 ENDDO 350 890 351 891 ! 352 892 !-- sgs-tke 353 CASE ( 6 ) 354 DO k = nzb_s_inner(j,i)+1, pch_index 355 tend(k,j,i) = tend(k,j,i) - & 356 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 357 SQRT( ( ( u(k,j,i) + & 358 u(k,j,i+1) ) & 359 / 2.0_wp )**2 + & 360 ( ( v(k,j,i) + & 361 v(k,j+1,i) ) & 362 / 2.0_wp )**2 + & 363 ( ( w(k,j,i) + & 364 w(k+1,j,i) ) & 365 / 2.0_wp )**2 ) * & 366 e(k,j,i) 367 ENDDO 893 CASE ( 6 ) 894 DO k = nzb_s_inner(j,i)+1, pch_index 895 tend(k,j,i) = tend(k,j,i) - & 896 2.0_wp * cdc * & 897 lad_s(k,j,i) * & 898 SQRT( ( 0.5_wp * ( u(k,j,i) + & 899 u(k,j,i+1) ) & 900 )**2 + & 901 ( 0.5_wp * ( v(k,j,i) + & 902 v(k,j+1,i) ) & 903 )**2 + & 904 ( 0.5_wp * ( w(k,j,i) + & 905 w(k+1,j,i) ) & 906 )**2 & 907 ) * & 908 e(k,j,i) 909 ENDDO 368 910 369 911 CASE DEFAULT
Note: See TracChangeset
for help on using the changeset viewer.
