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

Last change on this file since 2669 was 2669, checked in by raasch, 4 years ago

file attributes and activation strings in .palm.iofiles revised, file extensions for nesting, masked output, wind turbine data, etc. changed

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