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

Last change on this file since 2318 was 2318, checked in by suehring, 4 years ago

get topograpyhy top index via function call

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