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

Last change on this file since 2024 was 2024, checked in by kanani, 5 years ago

changes related to urban surface model and output of ssws

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