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

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

further modularization of radiation model and plant canopy model

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