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

Last change on this file since 2258 was 2233, checked in by suehring, 7 years ago

last commit documented

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