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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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