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

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

last commit documented

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