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

Last change on this file since 1834 was 1827, checked in by maronga, 5 years ago

last commit documented

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