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

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

last commit documented

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