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

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

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