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

Last change on this file since 2213 was 2213, checked in by kanani, 4 years ago

bugfix in PCM output, minor formatting and clean-up

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