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

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

Fixes for radiative transfer model

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