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

Last change on this file since 2932 was 2932, checked in by maronga, 3 years ago

renamed all Fortran NAMELISTS

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