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

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

Correction of "Former revisions" section

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