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

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

changes in the course of urban surface model implementation

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