source: palm/trunk/SOURCE/plant_canopy_model_mod.f90 @ 1981

Last change on this file since 1981 was 1961, checked in by suehring, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 54.1 KB
Line 
1!> @file plant_canopy_model_mod.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-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: plant_canopy_model_mod.f90 1961 2016-07-12 16:37:58Z suehring $
26!
27! 1960 2016-07-12 16:34:24Z suehring
28! Separate humidity and passive scalar
29!
30! 1953 2016-06-21 09:28:42Z suehring
31! Bugfix, lad_s and lad must be public
32!
33! 1826 2016-04-07 12:01:39Z maronga
34! Further modularization
35!
36! 1721 2015-11-16 12:56:48Z raasch
37! bugfixes: shf is reduced in areas covered with canopy only,
38!           canopy is set on top of topography
39!
40! 1682 2015-10-07 23:56:08Z knoop
41! Code annotations made doxygen readable
42!
43! 1484 2014-10-21 10:53:05Z kanani
44! Changes due to new module structure of the plant canopy model:
45!   module plant_canopy_model_mod now contains a subroutine for the
46!   initialization of the canopy model (pcm_init),
47!   limitation of the canopy drag (previously accounted for by calculation of
48!   a limiting canopy timestep for the determination of the maximum LES timestep
49!   in subroutine timestep) is now realized by the calculation of pre-tendencies
50!   and preliminary velocities in subroutine pcm_tendency,
51!   some redundant MPI communication removed in subroutine pcm_init
52!   (was previously in init_3d_model),
53!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
54!   respective grid is now provided only by lad_s (e.g. in the calculation of
55!   the tendency terms or of cum_lai_hf),
56!   drag_coefficient, lai, leaf_surface_concentration,
57!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
58!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
59!   respectively,
60!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
61!   USE-statements and ONLY-lists modified accordingly
62!
63! 1340 2014-03-25 19:45:13Z kanani
64! REAL constants defined as wp-kind
65!
66! 1320 2014-03-20 08:40:49Z raasch
67! ONLY-attribute added to USE-statements,
68! kind-parameters added to all INTEGER and REAL declaration statements,
69! kinds are defined in new module kinds,
70! old module precision_kind is removed,
71! revision history before 2012 removed,
72! comment fields (!:) to be used for variable explanations added to
73! all variable declaration statements
74!
75! 1036 2012-10-22 13:43:42Z raasch
76! code put under GPL (PALM 3.9)
77!
78! 138 2007-11-28 10:03:58Z letzel
79! Initial revision
80!
81! Description:
82! ------------
83!> 1) Initialization of the canopy model, e.g. construction of leaf area density
84!> profile (subroutine pcm_init).
85!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
86!> due to canopy elements (subroutine pcm_tendency).
87!------------------------------------------------------------------------------!
88 MODULE plant_canopy_model_mod
89 
90    USE arrays_3d,                                                             &
91        ONLY:  dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw 
92
93    USE indices,                                                               &
94        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
95               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
96
97    USE kinds
98
99
100    IMPLICIT NONE
101
102
103    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !< canopy coverage
104
105    INTEGER(iwp) ::  pch_index = 0                 !< plant canopy height/top index
106    INTEGER(iwp) ::                                                            &
107       lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
108
109    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
110    LOGICAL ::  plant_canopy = .FALSE.          !< switch for use of canopy model
111
112    REAL(wp) ::  alpha_lad = 9999999.9_wp   !< coefficient for lad calculation
113    REAL(wp) ::  beta_lad = 9999999.9_wp    !< coefficient for lad calculation
114    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
115    REAL(wp) ::  cdc = 0.0_wp               !< canopy drag coeff. (abbreviation used in equations)
116    REAL(wp) ::  cthf = 0.0_wp              !< canopy top heat flux
117    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !< timestep account. for canopy drag
118    REAL(wp) ::  lad_surface = 0.0_wp       !< lad surface value
119    REAL(wp) ::  lai_beta = 0.0_wp          !< leaf area index (lai) for lad calc.
120    REAL(wp) ::                                                                &
121       leaf_scalar_exch_coeff = 0.0_wp      !< canopy scalar exchange coeff.
122    REAL(wp) ::                                                                &
123       leaf_surface_conc = 0.0_wp           !< leaf surface concentration
124    REAL(wp) ::  lsec = 0.0_wp              !< leaf scalar exchange coeff.
125    REAL(wp) ::  lsc = 0.0_wp               !< leaf surface concentration
126
127    REAL(wp) ::                                                                &
128       lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
129    REAL(wp) ::                                                                &
130       lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
131
132    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
133    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
134   
135    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
136       canopy_heat_flux                                    !< canopy heat flux
137    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !< cumulative lai for heatflux calc.
138    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !< lad on scalar-grid
139
140
141    SAVE
142
143
144    PRIVATE
145 
146!
147!-- Public functions
148    PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
149
150!
151!-- Public variables and constants
152    PUBLIC canopy_mode, cthf, dt_plant_canopy, lad, lad_s, plant_canopy
153
154
155    INTERFACE pcm_check_parameters
156       MODULE PROCEDURE pcm_check_parameters
157    END INTERFACE pcm_check_parameters     
158   
159     INTERFACE pcm_header
160       MODULE PROCEDURE pcm_header
161    END INTERFACE pcm_header       
162   
163    INTERFACE pcm_init
164       MODULE PROCEDURE pcm_init
165    END INTERFACE pcm_init
166
167    INTERFACE pcm_parin
168       MODULE PROCEDURE pcm_parin
169    END INTERFACE pcm_parin 
170   
171    INTERFACE pcm_tendency
172       MODULE PROCEDURE pcm_tendency
173       MODULE PROCEDURE pcm_tendency_ij
174    END INTERFACE pcm_tendency
175
176
177 CONTAINS
178
179 
180!------------------------------------------------------------------------------!
181! Description:
182! ------------
183!> Check parameters routine for plant canopy model
184!------------------------------------------------------------------------------!
185    SUBROUTINE pcm_check_parameters
186
187       USE control_parameters,                                                 &
188           ONLY: cloud_physics, message_string, microphysics_seifert
189                 
190   
191       IMPLICIT NONE
192
193   
194       IF ( canopy_drag_coeff == 0.0_wp )  THEN
195          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
196                           'coefficient & given value is canopy_drag_coeff = 0.0'
197          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
198       ENDIF
199   
200       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
201              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
202          message_string = 'using the beta function for the construction ' //  &
203                           'of the leaf area density profile requires '    //  &
204                           'both alpha_lad and beta_lad to be /= 9999999.9'
205          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
206       ENDIF
207   
208       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
209          message_string = 'using the beta function for the construction ' //  &
210                           'of the leaf area density profile requires '    //  &
211                           'a non-zero lai_beta, but given value is '      //  &
212                           'lai_beta = 0.0'
213          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
214       ENDIF
215
216       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
217          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
218                           'and lad_surface /= 0.0 is not possible, '       // &
219                           'use either vertical gradients or the beta '     // &
220                           'function for the construction of the leaf area '// &
221                           'density profile'
222          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
223       ENDIF
224
225       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
226          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
227                          ' seifert_beheng'
228          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
229       ENDIF
230
231 
232    END SUBROUTINE pcm_check_parameters 
233 
234
235!------------------------------------------------------------------------------!
236! Description:
237! ------------
238!> Header output for plant canopy model
239!------------------------------------------------------------------------------!
240    SUBROUTINE pcm_header ( io )
241
242       USE control_parameters,                                                 &
243           ONLY: dz, passive_scalar
244
245
246       IMPLICIT NONE
247 
248       CHARACTER (LEN=10) ::  coor_chr            !<
249
250       CHARACTER (LEN=86) ::  coordinates         !<
251       CHARACTER (LEN=86) ::  gradients           !<
252       CHARACTER (LEN=86) ::  leaf_area_density   !<
253       CHARACTER (LEN=86) ::  slices              !<
254 
255       INTEGER(iwp) :: i                !<
256       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
257       INTEGER(iwp) :: k                !<       
258   
259       REAL(wp) ::  canopy_height       !< canopy height (in m)
260       
261       canopy_height = pch_index * dz
262
263       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
264                          canopy_drag_coeff
265       IF ( passive_scalar )  THEN
266          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
267                             leaf_surface_conc
268       ENDIF
269
270!
271!--    Heat flux at the top of vegetation
272       WRITE ( io, 3 )  cthf
273
274!
275!--    Leaf area density profile, calculated either from given vertical
276!--    gradients or from beta probability density function.
277       IF (  .NOT.  calc_beta_lad_profile )  THEN
278
279!--       Building output strings, starting with surface value
280          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
281          gradients = '------'
282          slices = '     0'
283          coordinates = '   0.0'
284          i = 1
285          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i)        &
286                      /= -9999 )
287
288             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
289             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
290                                 TRIM( coor_chr )
291 
292             WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
293             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
294
295             WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
296             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
297
298             WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
299             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
300
301             i = i + 1
302          ENDDO
303
304          WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
305                             TRIM( gradients ), TRIM( slices )
306
307       ELSE
308       
309          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
310          coordinates = '   0.0'
311         
312          DO  k = 1, pch_index
313
314             WRITE (coor_chr,'(F7.2)')  lad(k)
315             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
316                                 TRIM( coor_chr )
317 
318             WRITE (coor_chr,'(F7.1)')  zu(k)
319             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
320
321          ENDDO       
322
323          WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
324                          alpha_lad, beta_lad, lai_beta
325
326       ENDIF 
327
3281 FORMAT (//' Vegetation canopy (drag) model:'/                                &
329              ' ------------------------------'//                              &
330              ' Canopy mode: ', A /                                            &
331              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
332              ' Leaf drag coefficient: ',F6.2 /)
3332 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
334              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
3353 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
336          ' K m/s')
3374 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
338              ' Height:              ',A,'  m'/                                &
339              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
340              ' Gradient:            ',A,'  m**2/m**4'/                        &
341              ' Gridpoint:           ',A)
3425 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
343          //  ' Height:              ',A,'  m'/                                &
344              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
345              ' Coefficient alpha: ',F6.2 /                                    &
346              ' Coefficient beta: ',F6.2 /                                     &
347              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
348       
349    END SUBROUTINE pcm_header
350 
351 
352!------------------------------------------------------------------------------!
353! Description:
354! ------------
355!> Initialization of the plant canopy model
356!------------------------------------------------------------------------------!
357    SUBROUTINE pcm_init
358   
359
360       USE control_parameters,                                                 &
361           ONLY: dz, ocean, passive_scalar
362
363
364       IMPLICIT NONE
365
366       INTEGER(iwp) ::  i !< running index
367       INTEGER(iwp) ::  j !< running index
368       INTEGER(iwp) ::  k !< running index
369
370       REAL(wp) ::  int_bpdf      !< vertical integral for lad-profile construction
371       REAL(wp) ::  dzh           !< vertical grid spacing in units of canopy height
372       REAL(wp) ::  gradient      !< gradient for lad-profile construction
373       REAL(wp) ::  canopy_height !< canopy height for lad-profile construction
374
375!
376!--    Allocate one-dimensional arrays for the computation of the
377!--    leaf area density (lad) profile
378       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
379       lad = 0.0_wp
380       pre_lad = 0.0_wp
381
382!
383!--    Set flag that indicates that the lad-profile shall be calculated by using
384!--    a beta probability density function
385       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
386          calc_beta_lad_profile = .TRUE.
387       ENDIF
388       
389       
390!
391!--    Compute the profile of leaf area density used in the plant
392!--    canopy model. The profile can either be constructed from
393!--    prescribed vertical gradients of the leaf area density or by
394!--    using a beta probability density function (see e.g. Markkanen et al.,
395!--    2003: Boundary-Layer Meteorology, 106, 437-459)
396       IF (  .NOT.  calc_beta_lad_profile )  THEN   
397
398!
399!--       Use vertical gradients for lad-profile construction   
400          i = 1
401          gradient = 0.0_wp
402
403          IF (  .NOT.  ocean )  THEN
404
405             lad(0) = lad_surface
406             lad_vertical_gradient_level_ind(1) = 0
407 
408             DO k = 1, pch_index
409                IF ( i < 11 )  THEN
410                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
411                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
412                      gradient = lad_vertical_gradient(i)
413                      lad_vertical_gradient_level_ind(i) = k - 1
414                      i = i + 1
415                   ENDIF
416                ENDIF
417                IF ( gradient /= 0.0_wp )  THEN
418                   IF ( k /= 1 )  THEN
419                      lad(k) = lad(k-1) + dzu(k) * gradient
420                   ELSE
421                      lad(k) = lad_surface + dzu(k) * gradient
422                   ENDIF
423                ELSE
424                   lad(k) = lad(k-1)
425                ENDIF
426             ENDDO
427
428          ENDIF
429
430!
431!--       In case of no given leaf area density gradients, choose a vanishing
432!--       gradient. This information is used for the HEADER and the RUN_CONTROL
433!--       file.
434          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
435             lad_vertical_gradient_level(1) = 0.0_wp
436          ENDIF
437
438       ELSE
439
440!
441!--       Use beta function for lad-profile construction
442          int_bpdf = 0.0_wp
443          canopy_height = pch_index * dz
444
445          DO k = nzb, pch_index
446             int_bpdf = int_bpdf +                                             &
447                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
448                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
449                          beta_lad-1.0_wp ) )                                  &
450                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
451          ENDDO
452
453!
454!--       Preliminary lad profile (defined on w-grid)
455          DO k = nzb, pch_index
456             pre_lad(k) =  lai_beta *                                          &
457                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
458                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
459                              beta_lad-1.0_wp ) ) / int_bpdf                   &
460                        ) / canopy_height
461          ENDDO
462
463!
464!--       Final lad profile (defined on scalar-grid level, since most prognostic
465!--       quantities are defined there, hence, less interpolation is required
466!--       when calculating the canopy tendencies)
467          lad(0) = pre_lad(0)
468          DO k = nzb+1, pch_index
469             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
470          ENDDO         
471
472       ENDIF
473
474!
475!--    Allocate 3D-array for the leaf area density (lad_s). In case of a
476!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
477!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
478       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
479
480       IF ( cthf /= 0.0_wp )  THEN
481          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
482                    canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
483       ENDIF
484
485!
486!--    Initialize canopy parameters cdc (canopy drag coefficient),
487!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
488!--    with the prescribed values
489       cdc = canopy_drag_coeff
490       lsec = leaf_scalar_exch_coeff
491       lsc = leaf_surface_conc
492
493!
494!--    Initialization of the canopy coverage in the model domain:
495!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
496!--    fully covers the domain surface
497       SELECT CASE ( TRIM( canopy_mode ) )
498
499          CASE( 'block' )
500
501             DO  i = nxlg, nxrg
502                DO  j = nysg, nyng
503                   lad_s(:,j,i) = lad(:)
504                ENDDO
505             ENDDO
506
507          CASE DEFAULT
508
509!
510!--       The DEFAULT case is reached either if the parameter
511!--       canopy mode contains a wrong character string or if the
512!--       user has coded a special case in the user interface.
513!--       There, the subroutine user_init_plant_canopy checks
514!--       which of these two conditions applies.
515          CALL user_init_plant_canopy
516 
517       END SELECT
518
519!
520!--    Initialization of the canopy heat source distribution
521       IF ( cthf /= 0.0_wp )  THEN
522!
523!--       Piecewise calculation of the leaf area index by vertical
524!--       integration of the leaf area density
525          cum_lai_hf(:,:,:) = 0.0_wp
526          DO  i = nxlg, nxrg
527             DO  j = nysg, nyng
528                DO  k = pch_index-1, 0, -1
529                   IF ( k == pch_index-1 )  THEN
530                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
531                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
532                           ( zw(k+1) - zu(k+1) ) )  +                          &
533                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
534                                                 lad_s(k,j,i) ) +              &
535                                      lad_s(k+1,j,i) ) *                       &
536                           ( zu(k+1) - zw(k) ) ) 
537                   ELSE
538                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
539                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
540                                                 lad_s(k+1,j,i) ) +            &
541                                      lad_s(k+1,j,i) ) *                       &
542                           ( zw(k+1) - zu(k+1) ) )  +                          &
543                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
544                                                 lad_s(k,j,i) ) +              &
545                                      lad_s(k+1,j,i) ) *                       &
546                           ( zu(k+1) - zw(k) ) )
547                   ENDIF
548                ENDDO
549             ENDDO
550          ENDDO
551
552!
553!--       Calculation of the upward kinematic vertical heat flux within the
554!--       canopy
555          DO  i = nxlg, nxrg
556             DO  j = nysg, nyng
557                DO  k = 0, pch_index
558                   canopy_heat_flux(k,j,i) = cthf *                            &
559                                             exp( -0.6_wp * cum_lai_hf(k,j,i) )
560                ENDDO
561             ENDDO
562          ENDDO
563
564!
565!--       In areas covered with canopy, the surface heat flux is set to
566!--       the surface value of the above calculated in-canopy heat flux
567!--       distribution
568          DO  i = nxlg,nxrg
569             DO  j = nysg, nyng
570                IF ( canopy_heat_flux(0,j,i) /= cthf )  THEN
571                   shf(j,i) = canopy_heat_flux(0,j,i)
572                ENDIF
573             ENDDO
574          ENDDO
575
576       ENDIF
577
578
579
580    END SUBROUTINE pcm_init
581
582
583
584!------------------------------------------------------------------------------!
585! Description:
586! ------------
587!> Calculation of the tendency terms, accounting for the effect of the plant
588!> canopy on momentum and scalar quantities.
589!>
590!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
591!> (defined on scalar grid), as initialized in subroutine pcm_init.
592!> The lad on the w-grid is vertically interpolated from the surrounding
593!> lad_s. The upper boundary of the canopy is defined on the w-grid at
594!> k = pch_index. Here, the lad is zero.
595!>
596!> The canopy drag must be limited (previously accounted for by calculation of
597!> a limiting canopy timestep for the determination of the maximum LES timestep
598!> in subroutine timestep), since it is physically impossible that the canopy
599!> drag alone can locally change the sign of a velocity component. This
600!> limitation is realized by calculating preliminary tendencies and velocities.
601!> It is subsequently checked if the preliminary new velocity has a different
602!> sign than the current velocity. If so, the tendency is limited in a way that
603!> the velocity can at maximum be reduced to zero by the canopy drag.
604!>
605!>
606!> Call for all grid points
607!------------------------------------------------------------------------------!
608    SUBROUTINE pcm_tendency( component )
609
610
611       USE control_parameters,                                                 &
612           ONLY:  dt_3d, message_string
613
614       USE kinds
615
616       IMPLICIT NONE
617
618       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
619       INTEGER(iwp) ::  i         !< running index
620       INTEGER(iwp) ::  j         !< running index
621       INTEGER(iwp) ::  k         !< running index
622       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
623
624       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
625       REAL(wp) ::  lad_local !< local lad value
626       REAL(wp) ::  pre_tend  !< preliminary tendency
627       REAL(wp) ::  pre_u     !< preliminary u-value
628       REAL(wp) ::  pre_v     !< preliminary v-value
629       REAL(wp) ::  pre_w     !< preliminary w-value
630
631
632       ddt_3d = 1.0_wp / dt_3d
633 
634!
635!--    Compute drag for the three velocity components and the SGS-TKE:
636       SELECT CASE ( component )
637
638!
639!--       u-component
640          CASE ( 1 )
641             DO  i = nxlu, nxr
642                DO  j = nys, nyn
643                   DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
644
645                      kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
646!
647!--                   In order to create sharp boundaries of the plant canopy,
648!--                   the lad on the u-grid at index (k,j,i) is equal to
649!--                   lad_s(k,j,i), rather than being interpolated from the
650!--                   surrounding lad_s, because this would yield smaller lad
651!--                   at the canopy boundaries than inside of the canopy.
652!--                   For the same reason, the lad at the rightmost(i+1)canopy
653!--                   boundary on the u-grid equals lad_s(k,j,i).
654                      lad_local = lad_s(kk,j,i)
655                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
656                      THEN
657                         lad_local = lad_s(kk,j,i-1)
658                      ENDIF
659
660                      pre_tend = 0.0_wp
661                      pre_u = 0.0_wp
662!
663!--                   Calculate preliminary value (pre_tend) of the tendency
664                      pre_tend = - cdc *                                       &
665                                   lad_local *                                 &
666                                   SQRT( u(k,j,i)**2 +                         &
667                                         ( 0.25_wp * ( v(k,j,i-1) +            &
668                                                       v(k,j,i)   +            &
669                                                       v(k,j+1,i) +            &
670                                                       v(k,j+1,i-1) )          &
671                                         )**2 +                                &
672                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
673                                                       w(k-1,j,i)   +          &
674                                                       w(k,j,i-1)   +          &
675                                                       w(k,j,i) )              &
676                                         )**2                                  &
677                                       ) *                                     &
678                                   u(k,j,i)
679
680!
681!--                   Calculate preliminary new velocity, based on pre_tend
682                      pre_u = u(k,j,i) + dt_3d * pre_tend
683!
684!--                   Compare sign of old velocity and new preliminary velocity,
685!--                   and in case the signs are different, limit the tendency
686                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
687                         pre_tend = - u(k,j,i) * ddt_3d
688                      ELSE
689                         pre_tend = pre_tend
690                      ENDIF
691!
692!--                   Calculate final tendency
693                      tend(k,j,i) = tend(k,j,i) + pre_tend
694
695                   ENDDO
696                ENDDO
697             ENDDO
698
699!
700!--       v-component
701          CASE ( 2 )
702             DO  i = nxl, nxr
703                DO  j = nysv, nyn
704                   DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
705
706                      kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
707!
708!--                   In order to create sharp boundaries of the plant canopy,
709!--                   the lad on the v-grid at index (k,j,i) is equal to
710!--                   lad_s(k,j,i), rather than being interpolated from the
711!--                   surrounding lad_s, because this would yield smaller lad
712!--                   at the canopy boundaries than inside of the canopy.
713!--                   For the same reason, the lad at the northmost(j+1) canopy
714!--                   boundary on the v-grid equals lad_s(k,j,i).
715                      lad_local = lad_s(kk,j,i)
716                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
717                      THEN
718                         lad_local = lad_s(kk,j-1,i)
719                      ENDIF
720
721                      pre_tend = 0.0_wp
722                      pre_v = 0.0_wp
723!
724!--                   Calculate preliminary value (pre_tend) of the tendency
725                      pre_tend = - cdc *                                       &
726                                   lad_local *                                 &
727                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
728                                                       u(k,j-1,i+1) +          &
729                                                       u(k,j,i)     +          &
730                                                       u(k,j,i+1) )            &
731                                         )**2 +                                &
732                                         v(k,j,i)**2 +                         &
733                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
734                                                       w(k-1,j,i)   +          &
735                                                       w(k,j-1,i)   +          &
736                                                       w(k,j,i) )              &
737                                         )**2                                  &
738                                       ) *                                     &
739                                   v(k,j,i)
740
741!
742!--                   Calculate preliminary new velocity, based on pre_tend
743                      pre_v = v(k,j,i) + dt_3d * pre_tend
744!
745!--                   Compare sign of old velocity and new preliminary velocity,
746!--                   and in case the signs are different, limit the tendency
747                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
748                         pre_tend = - v(k,j,i) * ddt_3d
749                      ELSE
750                         pre_tend = pre_tend
751                      ENDIF
752!
753!--                   Calculate final tendency
754                      tend(k,j,i) = tend(k,j,i) + pre_tend
755
756                   ENDDO
757                ENDDO
758             ENDDO
759
760!
761!--       w-component
762          CASE ( 3 )
763             DO  i = nxl, nxr
764                DO  j = nys, nyn
765                   DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
766
767                      kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
768
769                      pre_tend = 0.0_wp
770                      pre_w = 0.0_wp
771!
772!--                   Calculate preliminary value (pre_tend) of the tendency
773                      pre_tend = - cdc *                                       &
774                                   (0.5_wp *                                   &
775                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
776                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
777                                                       u(k,j,i+1) +            &
778                                                       u(k+1,j,i) +            &
779                                                       u(k+1,j,i+1) )          &
780                                         )**2 +                                &
781                                         ( 0.25_wp * ( v(k,j,i)   +            &
782                                                       v(k,j+1,i) +            &
783                                                       v(k+1,j,i) +            &
784                                                       v(k+1,j+1,i) )          &
785                                         )**2 +                                &
786                                         w(k,j,i)**2                           &
787                                       ) *                                     &
788                                   w(k,j,i)
789!
790!--                   Calculate preliminary new velocity, based on pre_tend
791                      pre_w = w(k,j,i) + dt_3d * pre_tend
792!
793!--                   Compare sign of old velocity and new preliminary velocity,
794!--                   and in case the signs are different, limit the tendency
795                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
796                         pre_tend = - w(k,j,i) * ddt_3d
797                      ELSE
798                         pre_tend = pre_tend
799                      ENDIF
800!
801!--                   Calculate final tendency
802                      tend(k,j,i) = tend(k,j,i) + pre_tend
803
804                   ENDDO
805                ENDDO
806             ENDDO
807
808!
809!--       potential temperature
810          CASE ( 4 )
811             DO  i = nxl, nxr
812                DO  j = nys, nyn
813                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
814                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
815                      tend(k,j,i) = tend(k,j,i) +                              &
816                                       ( canopy_heat_flux(kk,j,i) -            &
817                                         canopy_heat_flux(kk-1,j,i) ) / dzw(k)
818                   ENDDO
819                ENDDO
820             ENDDO
821
822!
823!--       humidity
824          CASE ( 5 )
825             DO  i = nxl, nxr
826                DO  j = nys, nyn
827                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
828                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
829                      tend(k,j,i) = tend(k,j,i) -                              &
830                                       lsec *                                  &
831                                       lad_s(kk,j,i) *                         &
832                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
833                                                          u(k,j,i+1) )         &
834                                             )**2 +                            &
835                                             ( 0.5_wp * ( v(k,j,i) +           &
836                                                          v(k,j+1,i) )         &
837                                             )**2 +                            &
838                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
839                                                          w(k,j,i) )           &
840                                             )**2                              &
841                                           ) *                                 &
842                                       ( q(k,j,i) - lsc )
843                   ENDDO
844                ENDDO
845             ENDDO
846
847!
848!--       sgs-tke
849          CASE ( 6 )
850             DO  i = nxl, nxr
851                DO  j = nys, nyn
852                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
853                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
854                      tend(k,j,i) = tend(k,j,i) -                              &
855                                       2.0_wp * cdc *                          &
856                                       lad_s(kk,j,i) *                         &
857                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
858                                                          u(k,j,i+1) )         &
859                                             )**2 +                            &
860                                             ( 0.5_wp * ( v(k,j,i) +           &
861                                                          v(k,j+1,i) )         &
862                                             )**2 +                            &
863                                             ( 0.5_wp * ( w(k,j,i) +           &
864                                                          w(k+1,j,i) )         &
865                                             )**2                              &
866                                           ) *                                 &
867                                       e(k,j,i)
868                   ENDDO
869                ENDDO
870             ENDDO 
871!
872!--       scalar concentration
873          CASE ( 7 )
874             DO  i = nxl, nxr
875                DO  j = nys, nyn
876                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
877                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
878                      tend(k,j,i) = tend(k,j,i) -                              &
879                                       lsec *                                  &
880                                       lad_s(kk,j,i) *                         &
881                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
882                                                          u(k,j,i+1) )         &
883                                             )**2 +                            &
884                                             ( 0.5_wp * ( v(k,j,i) +           &
885                                                          v(k,j+1,i) )         &
886                                             )**2 +                            &
887                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
888                                                          w(k,j,i) )           &
889                                             )**2                              &
890                                           ) *                                 &
891                                       ( s(k,j,i) - lsc )
892                   ENDDO
893                ENDDO
894             ENDDO   
895
896
897
898          CASE DEFAULT
899
900             WRITE( message_string, * ) 'wrong component: ', component
901             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
902
903       END SELECT
904
905    END SUBROUTINE pcm_tendency
906
907
908!------------------------------------------------------------------------------!
909! Description:
910! ------------
911!> Parin for &canopy_par for plant canopy model
912!------------------------------------------------------------------------------!
913    SUBROUTINE pcm_parin
914
915
916       IMPLICIT NONE
917
918       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
919       
920        NAMELIST /canopy_par/     alpha_lad, beta_lad, canopy_drag_coeff,      &
921                                  canopy_mode, cthf,                           &
922                                  lad_surface,                                 & 
923                                  lad_vertical_gradient,                       &
924                                  lad_vertical_gradient_level,                 &
925                                  lai_beta,                                    &
926                                  leaf_scalar_exch_coeff,                      &
927                                  leaf_surface_conc, pch_index
928       
929       line = ' '
930       
931!
932!--    Try to find radiation model package
933       REWIND ( 11 )
934       line = ' '
935       DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
936          READ ( 11, '(A)', END=10 )  line
937       ENDDO
938       BACKSPACE ( 11 )
939
940!
941!--    Read user-defined namelist
942       READ ( 11, canopy_par )
943
944!
945!--    Set flag that indicates that the radiation model is switched on
946       plant_canopy = .TRUE.
947
948 10    CONTINUE
949       
950
951    END SUBROUTINE pcm_parin   
952   
953
954
955!------------------------------------------------------------------------------!
956! Description:
957! ------------
958!> Calculation of the tendency terms, accounting for the effect of the plant
959!> canopy on momentum and scalar quantities.
960!>
961!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
962!> (defined on scalar grid), as initialized in subroutine pcm_init.
963!> The lad on the w-grid is vertically interpolated from the surrounding
964!> lad_s. The upper boundary of the canopy is defined on the w-grid at
965!> k = pch_index. Here, the lad is zero.
966!>
967!> The canopy drag must be limited (previously accounted for by calculation of
968!> a limiting canopy timestep for the determination of the maximum LES timestep
969!> in subroutine timestep), since it is physically impossible that the canopy
970!> drag alone can locally change the sign of a velocity component. This
971!> limitation is realized by calculating preliminary tendencies and velocities.
972!> It is subsequently checked if the preliminary new velocity has a different
973!> sign than the current velocity. If so, the tendency is limited in a way that
974!> the velocity can at maximum be reduced to zero by the canopy drag.
975!>
976!>
977!> Call for grid point i,j
978!------------------------------------------------------------------------------!
979    SUBROUTINE pcm_tendency_ij( i, j, component )
980
981
982       USE control_parameters,                                                 &
983           ONLY:  dt_3d, message_string
984
985       USE kinds
986
987       IMPLICIT NONE
988
989       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
990       INTEGER(iwp) ::  i         !< running index
991       INTEGER(iwp) ::  j         !< running index
992       INTEGER(iwp) ::  k         !< running index
993       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
994
995       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
996       REAL(wp) ::  lad_local !< local lad value
997       REAL(wp) ::  pre_tend  !< preliminary tendency
998       REAL(wp) ::  pre_u     !< preliminary u-value
999       REAL(wp) ::  pre_v     !< preliminary v-value
1000       REAL(wp) ::  pre_w     !< preliminary w-value
1001
1002
1003       ddt_3d = 1.0_wp / dt_3d
1004
1005!
1006!--    Compute drag for the three velocity components and the SGS-TKE
1007       SELECT CASE ( component )
1008
1009!
1010!--       u-component
1011          CASE ( 1 )
1012             DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
1013
1014                kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
1015!
1016!--             In order to create sharp boundaries of the plant canopy,
1017!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1018!--             rather than being interpolated from the surrounding lad_s,
1019!--             because this would yield smaller lad at the canopy boundaries
1020!--             than inside of the canopy.
1021!--             For the same reason, the lad at the rightmost(i+1)canopy
1022!--             boundary on the u-grid equals lad_s(k,j,i).
1023                lad_local = lad_s(kk,j,i)
1024                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
1025                   lad_local = lad_s(kk,j,i-1)
1026                ENDIF
1027
1028                pre_tend = 0.0_wp
1029                pre_u = 0.0_wp
1030!
1031!--             Calculate preliminary value (pre_tend) of the tendency
1032                pre_tend = - cdc *                                             &
1033                             lad_local *                                       &   
1034                             SQRT( u(k,j,i)**2 +                               &
1035                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
1036                                                 v(k,j,i)    +                 &
1037                                                 v(k,j+1,i)  +                 &
1038                                                 v(k,j+1,i-1) )                &
1039                                   )**2 +                                      &
1040                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
1041                                                 w(k-1,j,i)   +                &
1042                                                 w(k,j,i-1)   +                &
1043                                                 w(k,j,i) )                    &
1044                                   )**2                                        &
1045                                 ) *                                           &
1046                             u(k,j,i)
1047
1048!
1049!--             Calculate preliminary new velocity, based on pre_tend
1050                pre_u = u(k,j,i) + dt_3d * pre_tend
1051!
1052!--             Compare sign of old velocity and new preliminary velocity,
1053!--             and in case the signs are different, limit the tendency
1054                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1055                   pre_tend = - u(k,j,i) * ddt_3d
1056                ELSE
1057                   pre_tend = pre_tend
1058                ENDIF
1059!
1060!--             Calculate final tendency
1061                tend(k,j,i) = tend(k,j,i) + pre_tend
1062             ENDDO
1063
1064
1065!
1066!--       v-component
1067          CASE ( 2 )
1068             DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
1069
1070                kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
1071!
1072!--             In order to create sharp boundaries of the plant canopy,
1073!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1074!--             rather than being interpolated from the surrounding lad_s,
1075!--             because this would yield smaller lad at the canopy boundaries
1076!--             than inside of the canopy.
1077!--             For the same reason, the lad at the northmost(j+1)canopy
1078!--             boundary on the v-grid equals lad_s(k,j,i).
1079                lad_local = lad_s(kk,j,i)
1080                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
1081                   lad_local = lad_s(kk,j-1,i)
1082                ENDIF
1083
1084                pre_tend = 0.0_wp
1085                pre_v = 0.0_wp
1086!
1087!--             Calculate preliminary value (pre_tend) of the tendency
1088                pre_tend = - cdc *                                             &
1089                             lad_local *                                       &
1090                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
1091                                                 u(k,j-1,i+1) +                &
1092                                                 u(k,j,i)     +                &
1093                                                 u(k,j,i+1) )                  &
1094                                   )**2 +                                      &
1095                                   v(k,j,i)**2 +                               &
1096                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
1097                                                 w(k-1,j,i)   +                &
1098                                                 w(k,j-1,i)   +                &
1099                                                 w(k,j,i) )                    &
1100                                   )**2                                        &
1101                                 ) *                                           &
1102                             v(k,j,i)
1103
1104!
1105!--             Calculate preliminary new velocity, based on pre_tend
1106                pre_v = v(k,j,i) + dt_3d * pre_tend
1107!
1108!--             Compare sign of old velocity and new preliminary velocity,
1109!--             and in case the signs are different, limit the tendency
1110                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1111                   pre_tend = - v(k,j,i) * ddt_3d
1112                ELSE
1113                   pre_tend = pre_tend
1114                ENDIF
1115!
1116!--             Calculate final tendency
1117                tend(k,j,i) = tend(k,j,i) + pre_tend
1118             ENDDO
1119
1120
1121!
1122!--       w-component
1123          CASE ( 3 )
1124             DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
1125
1126                kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
1127
1128                pre_tend = 0.0_wp
1129                pre_w = 0.0_wp
1130!
1131!--             Calculate preliminary value (pre_tend) of the tendency
1132                pre_tend = - cdc *                                             &
1133                             (0.5_wp *                                         &
1134                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
1135                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
1136                                                 u(k,j,i+1)  +                 &
1137                                                 u(k+1,j,i)  +                 &
1138                                                 u(k+1,j,i+1) )                &
1139                                   )**2 +                                      &
1140                                   ( 0.25_wp * ( v(k,j,i)    +                 &
1141                                                 v(k,j+1,i)  +                 &
1142                                                 v(k+1,j,i)  +                 &
1143                                                 v(k+1,j+1,i) )                &
1144                                   )**2 +                                      &
1145                                   w(k,j,i)**2                                 &
1146                                 ) *                                           &
1147                             w(k,j,i)
1148!
1149!--             Calculate preliminary new velocity, based on pre_tend
1150                pre_w = w(k,j,i) + dt_3d * pre_tend
1151!
1152!--             Compare sign of old velocity and new preliminary velocity,
1153!--             and in case the signs are different, limit the tendency
1154                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1155                   pre_tend = - w(k,j,i) * ddt_3d
1156                ELSE
1157                   pre_tend = pre_tend
1158                ENDIF
1159!
1160!--             Calculate final tendency
1161                tend(k,j,i) = tend(k,j,i) + pre_tend
1162             ENDDO
1163
1164!
1165!--       potential temperature
1166          CASE ( 4 )
1167             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1168                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1169                tend(k,j,i) = tend(k,j,i) +                                    &
1170                                 ( canopy_heat_flux(kk,j,i) -                  &
1171                                   canopy_heat_flux(kk-1,j,i) ) / dzw(k)
1172             ENDDO
1173
1174
1175!
1176!--       humidity
1177          CASE ( 5 )
1178             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1179                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1180                tend(k,j,i) = tend(k,j,i) -                                    &
1181                                 lsec *                                        &
1182                                 lad_s(kk,j,i) *                               &
1183                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1184                                                    u(k,j,i+1) )               &
1185                                       )**2  +                                 &
1186                                       ( 0.5_wp * ( v(k,j,i) +                 &
1187                                                    v(k,j+1,i) )               &
1188                                       )**2 +                                  &
1189                                       ( 0.5_wp * ( w(k-1,j,i) +               &
1190                                                    w(k,j,i) )                 &
1191                                       )**2                                    &
1192                                     ) *                                       &
1193                                 ( q(k,j,i) - lsc )
1194             ENDDO   
1195
1196!
1197!--       sgs-tke
1198          CASE ( 6 )
1199             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1200                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1201                tend(k,j,i) = tend(k,j,i) -                                    &
1202                                 2.0_wp * cdc *                                &
1203                                 lad_s(kk,j,i) *                               &
1204                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1205                                                    u(k,j,i+1) )               &
1206                                       )**2 +                                  & 
1207                                       ( 0.5_wp * ( v(k,j,i) +                 &
1208                                                    v(k,j+1,i) )               &
1209                                       )**2 +                                  &
1210                                       ( 0.5_wp * ( w(k,j,i) +                 &
1211                                                    w(k+1,j,i) )               &
1212                                       )**2                                    &
1213                                     ) *                                       &
1214                                 e(k,j,i)
1215             ENDDO
1216             
1217!
1218!--       scalar concentration
1219          CASE ( 7 )
1220             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1221                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1222                tend(k,j,i) = tend(k,j,i) -                                    &
1223                                 lsec *                                        &
1224                                 lad_s(kk,j,i) *                               &
1225                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1226                                                    u(k,j,i+1) )               &
1227                                       )**2  +                                 &
1228                                       ( 0.5_wp * ( v(k,j,i) +                 &
1229                                                    v(k,j+1,i) )               &
1230                                       )**2 +                                  &
1231                                       ( 0.5_wp * ( w(k-1,j,i) +               &
1232                                                    w(k,j,i) )                 &
1233                                       )**2                                    &
1234                                     ) *                                       &
1235                                 ( s(k,j,i) - lsc )
1236             ENDDO               
1237
1238       CASE DEFAULT
1239
1240          WRITE( message_string, * ) 'wrong component: ', component
1241          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1242
1243       END SELECT
1244
1245    END SUBROUTINE pcm_tendency_ij
1246
1247 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.