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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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