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

Last change on this file since 3163 was 3065, checked in by Giersch, 6 years ago

New vertical stretching procedure has been introduced

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