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

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

changes related to steering and formating of urban surface model

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