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

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

last commit documented

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