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

Last change on this file since 1484 was 1484, checked in by kanani, 7 years ago

New:
---
Subroutine init_plant_canopy added to module plant_canopy_model_mod. (plant_canopy_model)
Alternative method for lad-profile construction added, also, new parameters added.
(header, package_parin, plant_canopy_model, read_var_list, write_var_list)
plant_canopy_model-dependency added to several subroutines. (Makefile)
New package/namelist canopy_par for canopy-related parameters added. (package_parin)

Changed:
---
Code structure of the plant canopy model changed, all canopy-model related code
combined to module plant_canopy_model_mod. (check_parameters, init_3d_model,
modules, timestep)
Module plant_canopy_model_mod added in USE-lists of some subroutines. (check_parameters,
header, init_3d_model, package_parin, read_var_list, user_init_plant_canopy, write_var_list)
Canopy initialization moved to new subroutine init_plant_canopy. (check_parameters,
init_3d_model, plant_canopy_model)
Calculation of canopy timestep-criterion removed, instead, the canopy
drag is now directly limited in the calculation of the canopy tendency terms.
(plant_canopy_model, timestep)
Some parameters renamed. (check_parameters, header, init_plant_canopy,
plant_canopy_model, read_var_list, write_var_list)
Unnecessary 3d-arrays removed. (init_plant_canopy, plant_canopy_model, user_init_plant_canopy)
Parameter checks regarding canopy initialization added. (check_parameters)
All canopy steering parameters moved from namelist inipar to canopy_par. (package_parin, parin)
Some redundant MPI communication removed. (init_plant_canopy)

Bugfix:
---
Missing KIND-attribute for REAL constant added. (check_parameters)
DO-WHILE-loop for lad-profile output restricted. (header)
Removed double-listing of use_upstream_for_tke in ONLY-list of module
control_parameters. (prognostic_equations)

  • Property svn:keywords set to Id
File size: 40.7 KB
Line 
1 MODULE plant_canopy_model_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
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
40!
41! Former revisions:
42! -----------------
43! $Id: plant_canopy_model.f90 1484 2014-10-21 10:53:05Z kanani $
44!
45! 1340 2014-03-25 19:45:13Z kanani
46! REAL constants defined as wp-kind
47!
48! 1320 2014-03-20 08:40:49Z raasch
49! ONLY-attribute added to USE-statements,
50! kind-parameters added to all INTEGER and REAL declaration statements,
51! kinds are defined in new module kinds,
52! old module precision_kind is removed,
53! revision history before 2012 removed,
54! comment fields (!:) to be used for variable explanations added to
55! all variable declaration statements
56!
57! 1036 2012-10-22 13:43:42Z raasch
58! code put under GPL (PALM 3.9)
59!
60! 138 2007-11-28 10:03:58Z letzel
61! Initial revision
62!
63! Description:
64! ------------
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).
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
123
124    PRIVATE
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
136
137    INTERFACE plant_canopy_model
138       MODULE PROCEDURE plant_canopy_model
139       MODULE PROCEDURE plant_canopy_model_ij
140    END INTERFACE plant_canopy_model
141
142
143
144
145 CONTAINS
146
147
148!------------------------------------------------------------------------------!
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
387!------------------------------------------------------------------------------!
388    SUBROUTINE plant_canopy_model( component )
389
390
391       USE control_parameters,                                                 &
392           ONLY:  dt_3d, message_string
393
394       USE kinds
395
396       IMPLICIT NONE
397
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
412 
413!
414!--    Compute drag for the three velocity components and the SGS-TKE:
415       SELECT CASE ( component )
416
417!
418!--       u-component
419          CASE ( 1 )
420             DO  i = nxlu, nxr
421                DO  j = nys, nyn
422                   DO  k = nzb_u_inner(j,i)+1, pch_index
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
473                   ENDDO
474                ENDDO
475             ENDDO
476
477!
478!--       v-component
479          CASE ( 2 )
480             DO  i = nxl, nxr
481                DO  j = nysv, nyn
482                   DO  k = nzb_v_inner(j,i)+1, pch_index
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
533                   ENDDO
534                ENDDO
535             ENDDO
536
537!
538!--       w-component
539          CASE ( 3 )
540             DO  i = nxl, nxr
541                DO  j = nys, nyn
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
579                   ENDDO
580                ENDDO
581             ENDDO
582
583!
584!--       potential temperature
585          CASE ( 4 )
586             DO  i = nxl, nxr
587                DO  j = nys, nyn
588                   DO  k = nzb_s_inner(j,i)+1, pch_index
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)
592                   ENDDO
593                ENDDO
594             ENDDO
595
596!
597!--       scalar concentration
598          CASE ( 5 )
599             DO  i = nxl, nxr
600                DO  j = nys, nyn
601                   DO  k = nzb_s_inner(j,i)+1, pch_index
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 )
616                   ENDDO
617                ENDDO
618             ENDDO
619
620!
621!--       sgs-tke
622          CASE ( 6 )
623             DO  i = nxl, nxr
624                DO  j = nys, nyn
625                   DO  k = nzb_s_inner(j,i)+1, pch_index
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)
640                   ENDDO
641                ENDDO
642             ENDDO
643
644
645          CASE DEFAULT
646
647             WRITE( message_string, * ) 'wrong component: ', component
648             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
649
650       END SELECT
651
652    END SUBROUTINE plant_canopy_model
653
654
655!------------------------------------------------------------------------------!
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
678!------------------------------------------------------------------------------!
679    SUBROUTINE plant_canopy_model_ij( i, j, component )
680
681
682       USE control_parameters,                                                 &
683           ONLY:  dt_3d, message_string
684
685       USE kinds
686
687       IMPLICIT NONE
688
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
706       SELECT CASE ( component )
707
708!
709!--       u-component
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
763
764!
765!--       v-component
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
819
820!
821!--       w-component
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
860
861!
862!--       potential temperature
863          CASE ( 4 )
864             DO  k = nzb_s_inner(j,i)+1, pch_index
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)
868             ENDDO
869
870
871!
872!--       scalar concentration
873          CASE ( 5 )
874             DO  k = nzb_s_inner(j,i)+1, pch_index
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 )
889             ENDDO   
890
891!
892!--       sgs-tke
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
910
911       CASE DEFAULT
912
913          WRITE( message_string, * ) 'wrong component: ', component
914          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
915
916       END SELECT
917
918    END SUBROUTINE plant_canopy_model_ij
919
920 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.