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

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

upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny; no output if redundant ghost layer data to NetCDF files

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