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

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

Particle reflections at downward-facing walls; revision of particle speed interpolations at walls; bugfixes in get_topography_index and in date constants

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