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

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

last commit documented

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