source: palm/trunk/SOURCE/plant_canopy_model.f90 @ 1805

Last change on this file since 1805 was 1722, checked in by raasch, 8 years ago

last commit documented

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