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

Last change on this file since 2274 was 2274, checked in by Giersch, 4 years ago

Changed error messages

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