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

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

Bugfix in message calls for local checks; error messages in init_grid slightly revised; bugfix in time_integration (uninitialized emission time index); read_restart_data (change from J.Resler): use dynamic array allocation rather than automatic arrays to avoid problems with stack memory

  • Property svn:keywords set to Id
File size: 97.3 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-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Institute of Computer Science of the
19!                     Czech Academy of Sciences, Prague
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: plant_canopy_model_mod.f90 4356 2019-12-20 17:09:33Z suehring $
29! Correct single message call, local check must be given by the respective
30! mpi rank.
31!
32! 4346 2019-12-18 11:55:56Z motisi
33! Introduction of wall_flags_total_0, which currently sets bits based on static
34! topography information used in wall_flags_static_0
35!
36! 4342 2019-12-16 13:49:14Z Giersch
37! Use statements moved to module level, ocean dependency removed, redundant
38! variables removed
39!
40! 4341 2019-12-16 10:43:49Z motisi
41! - Unification of variable names: pc_-variables now pcm_-variables
42!   (pc_latent_rate, pc_heating_rate, pc_transpiration_rate)
43! - Removal of pcm_bowenratio output
44! - Renamed canopy-mode 'block' to 'homogeneous'
45! - Renamed value 'read_from_file_3d' to 'read_from_file'
46! - Removal of confusing comment lines
47! - Replacement of k_wall by topo_top_ind
48! - Removal of Else-Statement in tendency-calculation
49!
50! 4335 2019-12-12 16:39:05Z suehring
51! Fix for LAD at building edges also implemented in vector branch.
52!
53! 4331 2019-12-10 18:25:02Z suehring
54! Typo corrected
55!
56! 4329 2019-12-10 15:46:36Z motisi
57! Renamed wall_flags_0 to wall_flags_static_0
58!
59! 4314 2019-11-29 10:29:20Z suehring
60! - Bugfix, plant canopy was still considered at building edges on for the u-
61!   and v-component.
62! - Relax restriction of LAD on building tops. LAD is only omitted at
63!   locations where building grid points emerged artificially by the
64!   topography filtering.
65!
66! 4309 2019-11-26 18:49:59Z suehring
67! Typo
68!
69! 4302 2019-11-22 13:15:56Z suehring
70! Omit tall canopy mapped on top of buildings
71!
72! 4279 2019-10-29 08:48:17Z scharf
73! unused variables removed
74!
75! 4258 2019-10-07 13:29:08Z scharf
76! changed check for static driver and fixed bugs in initialization and header
77!
78! 4258 2019-10-07 13:29:08Z suehring
79! Check if any LAD is prescribed when plant-canopy model is applied.
80!
81! 4226 2019-09-10 17:03:24Z suehring
82! Bugfix, missing initialization of heating rate
83!
84! 4221 2019-09-09 08:50:35Z suehring
85! Further bugfix in 3d data output for plant canopy
86!
87! 4216 2019-09-04 09:09:03Z suehring
88! Bugfixes in 3d data output
89!
90! 4205 2019-08-30 13:25:00Z suehring
91! Missing working precision + bugfix in calculation of wind speed
92!
93! 4188 2019-08-26 14:15:47Z suehring
94! Minor adjustment in error number
95!
96! 4187 2019-08-26 12:43:15Z suehring
97! Give specific error numbers instead of PA0999
98!
99! 4182 2019-08-22 15:20:23Z scharf
100! Corrected "Former revisions" section
101!
102! 4168 2019-08-16 13:50:17Z suehring
103! Replace function get_topography_top_index by topo_top_ind
104!
105! 4127 2019-07-30 14:47:10Z suehring
106! Output of 3D plant canopy variables changed. It is now relative to the local
107! terrain rather than located at the acutal vertical level in the model. This
108! way, the vertical dimension of the output can be significantly reduced.
109! (merge from branch resler)
110!
111! 3885 2019-04-11 11:29:34Z kanani
112! Changes related to global restructuring of location messages and introduction
113! of additional debug messages
114!
115! 3864 2019-04-05 09:01:56Z monakurppa
116! unsed variables removed
117!
118! 3745 2019-02-15 18:57:56Z suehring
119! Bugfix in transpiration, floating invalid when temperature
120! becomes > 40 degrees
121!
122! 3744 2019-02-15 18:38:58Z suehring
123! Some interface calls moved to module_interface + cleanup
124!
125! 3655 2019-01-07 16:51:22Z knoop
126! unused variables removed
127!
128! 138 2007-11-28 10:03:58Z letzel
129! Initial revision
130!
131! Description:
132! ------------
133!> 1) Initialization of the canopy model, e.g. construction of leaf area density
134!> profile (subroutine pcm_init).
135!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
136!> due to canopy elements (subroutine pcm_tendency).
137!
138! @todo - precalculate constant terms in pcm_calc_transpiration_rate
139! @todo - unify variable names (pcm_, pc_, ...)
140!------------------------------------------------------------------------------!
141 MODULE plant_canopy_model_mod
142 
143    USE arrays_3d,                                                             &
144        ONLY:  dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
145
146    USE basic_constants_and_equations_mod,                                     &
147        ONLY:  c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
148       
149    USE bulk_cloud_model_mod,                                                  &
150        ONLY: bulk_cloud_model, microphysics_seifert
151
152    USE control_parameters,                                                    &
153        ONLY: average_count_3d, coupling_char, debug_output, dt_3d, dz,        &
154              humidity, message_string, ocean_mode, passive_scalar,            &
155              plant_canopy, urban_surface
156             
157    USE grid_variables,                                                        &
158        ONLY:  dx, dy
159
160    USE indices,                                                               &
161        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
162               nz, nzb, nzt, topo_top_ind, wall_flags_total_0
163
164    USE kinds
165   
166    USE netcdf_data_input_mod,                                                 &
167        ONLY:  input_pids_static, leaf_area_density_f
168
169    USE pegrid
170       
171    USE surface_mod,                                                           &
172        ONLY: surf_def_h, surf_lsm_h, surf_usm_h
173
174
175    IMPLICIT NONE
176
177    CHARACTER (LEN=30)   ::  canopy_mode = 'homogeneous'           !< canopy coverage
178    LOGICAL              ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
179                                                                   !< for resolved plant canopy inside radiation model
180                                                                   !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
181
182    INTEGER(iwp) ::  pch_index = 0                                 !< plant canopy height/top index
183    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999   !< lad-profile levels (index)
184
185    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji     !< local plant canopy top
186
187    LOGICAL ::  calc_beta_lad_profile = .FALSE.   !< switch for calc. of lad from beta func.
188
189    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
190    REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
191    REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
192    REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
193    REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
194    REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
195    REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
196    REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
197    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
198    REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
199
200    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
201    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
202
203    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
204                                                !< of plant canopy (e.g. deciduous tree during winter)
205
206    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
207    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
208   
209    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf               !< cumulative lai for heatflux calc.
210    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
211    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heating_rate         !< plant canopy heating rate
212    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpiration_rate   !< plant canopy transpiration rate
213    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latent_rate          !< plant canopy latent heating rate
214
215    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heatrate_av          !< array for averaging plant canopy sensible heating rate
216    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latentrate_av        !< array for averaging plant canopy latent heating rate
217    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
218
219    SAVE
220
221    PRIVATE
222   
223!
224!-- Public functions
225    PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output,                &
226           pcm_check_parameters, pcm_3d_data_averaging,                       &
227           pcm_data_output_3d, pcm_define_netcdf_grid,                        &
228           pcm_header, pcm_init, pcm_parin, pcm_tendency
229
230!
231!-- Public variables and constants
232    PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate,       &
233           pcm_latent_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
234           pch_index, plant_canopy_transpiration
235
236    INTERFACE pcm_calc_transpiration_rate
237       MODULE PROCEDURE pcm_calc_transpiration_rate
238    END INTERFACE pcm_calc_transpiration_rate
239
240    INTERFACE pcm_check_data_output
241       MODULE PROCEDURE pcm_check_data_output
242    END INTERFACE pcm_check_data_output
243   
244    INTERFACE pcm_check_parameters
245       MODULE PROCEDURE pcm_check_parameters
246    END INTERFACE pcm_check_parameters
247
248    INTERFACE pcm_3d_data_averaging
249       MODULE PROCEDURE pcm_3d_data_averaging
250    END INTERFACE pcm_3d_data_averaging
251
252    INTERFACE pcm_data_output_3d
253       MODULE PROCEDURE pcm_data_output_3d
254    END INTERFACE pcm_data_output_3d
255
256    INTERFACE pcm_define_netcdf_grid
257       MODULE PROCEDURE pcm_define_netcdf_grid
258    END INTERFACE pcm_define_netcdf_grid
259   
260     INTERFACE pcm_header
261       MODULE PROCEDURE pcm_header
262    END INTERFACE pcm_header       
263   
264    INTERFACE pcm_init
265       MODULE PROCEDURE pcm_init
266    END INTERFACE pcm_init
267
268    INTERFACE pcm_parin
269       MODULE PROCEDURE pcm_parin
270    END INTERFACE pcm_parin
271
272    INTERFACE pcm_read_plant_canopy_3d
273       MODULE PROCEDURE pcm_read_plant_canopy_3d
274    END INTERFACE pcm_read_plant_canopy_3d
275   
276    INTERFACE pcm_tendency
277       MODULE PROCEDURE pcm_tendency
278       MODULE PROCEDURE pcm_tendency_ij
279    END INTERFACE pcm_tendency
280
281
282 CONTAINS
283 
284 
285!------------------------------------------------------------------------------!
286! Description:
287! ------------
288!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
289!> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
290!> Meteorol. 97) and Ngao, Adam and Saudreau (2017;  Agricult. and Forest Meteorol
291!> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
292!> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
293!> of values. Funcion f4 used in form present in van Wijk et al. (1998;
294!> Tree Physiology 20).
295!>
296!> This subroutine is called from subroutine radiation_interaction
297!> after the calculation of radiation in plant canopy boxes.
298!> (arrays pcbinsw and pcbinlw).
299!>
300!------------------------------------------------------------------------------!
301 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
302
303!
304!--  input parameters
305     INTEGER(iwp), INTENT(IN) ::  i, j, k, kk        !< indices of the pc gridbox
306     REAL(wp), INTENT(IN)     ::  pcbsw              !< sw radiation in gridbox (W)
307     REAL(wp), INTENT(IN)     ::  pcblw              !< lw radiation in gridbox (W)
308     REAL(wp), INTENT(OUT)    ::  pcbtr              !< transpiration rate dq/dt (kg/kg/s)
309     REAL(wp), INTENT(OUT)    ::  pcblh              !< latent heat from transpiration dT/dt (K/s)
310
311!--  variables and parameters for calculation of transpiration rate
312     REAL(wp)             ::  sat_press, sat_press_d, temp, v_lad
313     REAL(wp)             ::  d_fact, g_b, g_s, wind_speed, evapor_rate
314     REAL(wp)             ::  f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
315     REAL(wp), PARAMETER  ::  gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
316     REAL(wp), PARAMETER  ::  g_s_max = 0.01        !< maximum stomatal conductivity (m/s)
317     REAL(wp), PARAMETER  ::  m_soil = 0.4_wp       !< soil water content (needs to adjust or take from LSM)
318     REAL(wp), PARAMETER  ::  m_wilt = 0.01_wp      !< wilting point soil water content (needs to adjust or take from LSM)
319     REAL(wp), PARAMETER  ::  m_sat = 0.51_wp       !< saturation soil water content (needs to adjust or take from LSM)
320     REAL(wp), PARAMETER  ::  t2_min = 0.0_wp       !< minimal temperature for calculation of f2
321     REAL(wp), PARAMETER  ::  t2_max = 40.0_wp      !< maximal temperature for calculation of f2
322
323
324!--  Temperature (deg C)
325     temp = pt(k,j,i) * exner(k) - degc_to_k
326!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
327     v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
328!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
329!--  There are updated formulas available, kept consistent with the rest of the parametrization
330     sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
331!--  Saturation pressure derivative (derivative of the above)
332     sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
333!--  Wind speed
334     wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +            &
335                        ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 +            &
336                        ( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
337!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
338     g_b = 0.01_wp * wind_speed + 0.0071_wp
339!--  Radiation flux per leaf surface unit
340     rad = pcbsw * v_lad
341!--  First function for calculation of stomatal conductivity (radiation dependency)
342!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
343     f1 = rad * (1000.0_wp+42.1_wp) / 1000.0_wp / (rad+42.1_wp)
344!--  Second function for calculation of stomatal conductivity (temperature dependency)
345!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
346     f2 = MAX(t2_min, (temp-t2_min) * MAX(0.0_wp,t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
347              ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
348!--  Water pressure deficit
349!--  Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
350     vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
351!--  Third function for calculation of stomatal conductivity (water pressure deficit dependency)
352!--  Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
353!--  The coefficients of the linear dependence should better correspond to broad-leaved trees
354!--  than the coefficients from Stewart (1988) which correspond to conifer trees.
355     vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
356     f3 = -2E-4_wp * vpd + 1.154_wp
357!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
358!--  Residual soil water content
359!--  van Wijk et al. (1998; Tree Physiology 20) eq. 7
360!--  TODO - over LSM surface might be calculated from LSM parameters
361     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
362!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
363     f4 = MAX(0.0_wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
364!--  Stomatal conductivity
365!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
366!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
367     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
368!--  Decoupling factor
369!--  Daudet et al. (1999) eq. 6
370     d_fact = (sat_press_d / gama_psychr + 2.0_wp ) /                        &
371              (sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
372!--  Equilibrium evaporation rate
373!--  Daudet et al. (1999) eq. 4
374     e_eq = (pcbsw + pcblw) * v_lad * sat_press_d /                         &
375                 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
376!--  Imposed evaporation rate
377!--  Daudet et al. (1999) eq. 5
378     e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
379!--  Evaporation rate
380!--  Daudet et al. (1999) eq. 3
381!--  (evaporation rate is limited to non-negative values)
382     evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
383!--  Conversion of evaporation rate to q tendency in gridbox
384!--  dq/dt = E * LAD * V_g / (rho_air * V_g)
385     pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k)  !-- = dq/dt
386!--  latent heat from evaporation
387     pcblh = pcbtr * lv_d_cp  !-- = - dT/dt
388
389 END SUBROUTINE pcm_calc_transpiration_rate
390
391
392!------------------------------------------------------------------------------!
393! Description:
394! ------------
395!> Check data output for plant canopy model
396!------------------------------------------------------------------------------!
397 SUBROUTINE pcm_check_data_output( var, unit )
398
399    CHARACTER (LEN=*) ::  unit  !<
400    CHARACTER (LEN=*) ::  var   !<
401
402
403    SELECT CASE ( TRIM( var ) )
404
405       CASE ( 'pcm_heatrate' )
406          IF ( cthf == 0.0_wp  .AND. .NOT.  urban_surface )  THEN
407             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
408                              'res setting of parameter cthf /= 0.0'
409             CALL message( 'pcm_check_data_output', 'PA1000', 1, 2, 0, 6, 0 )
410          ENDIF
411          unit = 'K s-1'
412   
413       CASE ( 'pcm_transpirationrate' )
414          unit = 'kg kg-1 s-1'
415
416       CASE ( 'pcm_latentrate' )
417          unit = 'K s-1'
418
419       CASE ( 'pcm_lad' )
420          unit = 'm2 m-3'
421
422
423       CASE DEFAULT
424          unit = 'illegal'
425
426    END SELECT
427
428
429 END SUBROUTINE pcm_check_data_output
430 
431 
432!------------------------------------------------------------------------------!
433! Description:
434! ------------
435!> Check parameters routine for plant canopy model
436!------------------------------------------------------------------------------!
437    SUBROUTINE pcm_check_parameters
438           
439       IF ( ocean_mode )  THEN
440          message_string = 'plant_canopy = .TRUE. is not allowed in the '//    &
441                           'ocean'
442          CALL message( 'pcm_check_parameters', 'PA0696', 1, 2, 0, 6, 0 )
443       ENDIF
444           
445       IF ( canopy_drag_coeff == 0.0_wp )  THEN
446          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
447                           'coefficient & given value is canopy_drag_coeff = 0.0'
448          CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
449       ENDIF
450
451       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp ) .OR.&
452              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
453          message_string = 'using the beta function for the construction ' //  &
454                           'of the leaf area density profile requires '    //  &
455                           'both alpha_lad and beta_lad to be /= 9999999.9'
456          CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
457       ENDIF
458
459       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
460          message_string = 'using the beta function for the construction ' //  &
461                           'of the leaf area density profile requires '    //  &
462                           'a non-zero lai_beta, but given value is '      //  &
463                           'lai_beta = 0.0'
464          CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
465       ENDIF
466
467       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
468          message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
469                           'combined with beta_lad /= 9999999.9 '           // &
470                           'and lad_surface /= 0.0 is not possible, '       // &
471                           'use either vertical gradients or the beta '     // &
472                           'function for the construction of the leaf area '// &
473                           'density profile'
474          CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
475       ENDIF
476
477       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
478          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
479                          ' seifert_beheng'
480          CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
481       ENDIF
482!
483!--    If canopy shall be read from file, static input file must be present
484       IF ( TRIM( canopy_mode ) == 'read_from_file' .AND.                   &
485            .NOT. input_pids_static )  THEN
486          message_string = 'canopy_mode = read_from_file requires ' //      &
487                           'static input file'
488          CALL message( 'pcm_check_parameters', 'PA0672', 1, 2, 0, 6, 0 )
489       ENDIF
490
491    END SUBROUTINE pcm_check_parameters 
492 
493
494!------------------------------------------------------------------------------!
495!
496! Description:
497! ------------
498!> Subroutine for averaging 3D data
499!------------------------------------------------------------------------------!
500 SUBROUTINE pcm_3d_data_averaging( mode, variable )
501
502    CHARACTER (LEN=*) ::  mode    !<
503    CHARACTER (LEN=*) :: variable !<
504
505    INTEGER(iwp) ::  i            !<
506    INTEGER(iwp) ::  j            !<
507    INTEGER(iwp) ::  k            !<
508
509
510    IF ( mode == 'allocate' )  THEN
511
512       SELECT CASE ( TRIM( variable ) )
513
514          CASE ( 'pcm_heatrate' )
515             IF ( .NOT. ALLOCATED( pcm_heatrate_av ) )  THEN
516                ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
517             ENDIF
518             pcm_heatrate_av = 0.0_wp
519
520
521          CASE ( 'pcm_latentrate' )
522             IF ( .NOT. ALLOCATED( pcm_latentrate_av ) )  THEN
523                ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
524             ENDIF
525             pcm_latentrate_av = 0.0_wp
526
527
528          CASE ( 'pcm_transpirationrate' )
529             IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) )  THEN
530                ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
531             ENDIF
532             pcm_transpirationrate_av = 0.0_wp
533
534          CASE DEFAULT
535             CONTINUE
536
537       END SELECT
538
539    ELSEIF ( mode == 'sum' )  THEN
540
541       SELECT CASE ( TRIM( variable ) )
542
543          CASE ( 'pcm_heatrate' )
544             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
545                DO  i = nxl, nxr
546                   DO  j = nys, nyn
547                      IF ( pch_index_ji(j,i) /= 0 )  THEN
548                         DO  k = 0, pch_index_ji(j,i)
549                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + pcm_heating_rate(k,j,i)
550                         ENDDO
551                      ENDIF
552                   ENDDO
553                ENDDO
554             ENDIF
555
556
557          CASE ( 'pcm_latentrate' )
558             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
559                DO  i = nxl, nxr
560                   DO  j = nys, nyn
561                      IF ( pch_index_ji(j,i) /= 0 )  THEN
562                         DO  k = 0, pch_index_ji(j,i)
563                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + pcm_latent_rate(k,j,i)
564                         ENDDO
565                      ENDIF
566                   ENDDO
567                ENDDO
568             ENDIF
569
570
571          CASE ( 'pcm_transpirationrate' )
572             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
573                DO  i = nxl, nxr
574                   DO  j = nys, nyn
575                      IF ( pch_index_ji(j,i) /= 0 )  THEN
576                         DO  k = 0, pch_index_ji(j,i)
577                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + pcm_transpiration_rate(k,j,i)
578                         ENDDO
579                      ENDIF
580                   ENDDO
581                ENDDO
582             ENDIF
583
584          CASE DEFAULT
585             CONTINUE
586
587       END SELECT
588
589    ELSEIF ( mode == 'average' )  THEN
590
591       SELECT CASE ( TRIM( variable ) )
592
593          CASE ( 'pcm_heatrate' )
594             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
595                DO  i = nxlg, nxrg
596                   DO  j = nysg, nyng
597                      IF ( pch_index_ji(j,i) /= 0 )  THEN
598                         DO  k = 0, pch_index_ji(j,i)
599                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i)                 &
600                                                     / REAL( average_count_3d, KIND=wp )
601                         ENDDO
602                      ENDIF
603                   ENDDO
604                ENDDO
605             ENDIF
606
607
608          CASE ( 'pcm_latentrate' )
609             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
610                DO  i = nxlg, nxrg
611                   DO  j = nysg, nyng
612                      IF ( pch_index_ji(j,i) /= 0 )  THEN
613                         DO  k = 0, pch_index_ji(j,i)
614                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i)              &
615                                                       / REAL( average_count_3d, KIND=wp )
616                         ENDDO
617                      ENDIF
618                   ENDDO
619                ENDDO
620             ENDIF
621
622
623          CASE ( 'pcm_transpirationrate' )
624             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
625                DO  i = nxlg, nxrg
626                   DO  j = nysg, nyng
627                      IF ( pch_index_ji(j,i) /= 0 )  THEN
628                         DO  k = 0, pch_index_ji(j,i)
629                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i)  &
630                                                              / REAL( average_count_3d, KIND=wp )
631                         ENDDO
632                      ENDIF
633                   ENDDO
634                ENDDO
635             ENDIF
636
637       END SELECT
638
639    ENDIF
640
641 END SUBROUTINE pcm_3d_data_averaging
642
643!------------------------------------------------------------------------------!
644!
645! Description:
646! ------------
647!> Subroutine defining 3D output variables.
648!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning
649!> that 3D output is relative to the model surface now rather than at the actual
650!> grid point where the plant canopy is located.
651!------------------------------------------------------------------------------!
652 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
653                                nzb_do, nzt_do )
654
655    CHARACTER (LEN=*) ::  variable !< treated variable
656
657    INTEGER(iwp) ::  av     !< flag indicating instantaneous or averaged data output
658    INTEGER(iwp) ::  i      !< grid index x-direction
659    INTEGER(iwp) ::  j      !< grid index y-direction
660    INTEGER(iwp) ::  k      !< grid index z-direction
661    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
662    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
663
664    LOGICAL      ::  found  !< flag indicating if variable is found
665
666    REAL(wp)     ::  fill_value !< fill value
667    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< data output array
668
669
670    found = .TRUE.
671
672    local_pf = REAL( fill_value, KIND = 4 )
673
674    SELECT CASE ( TRIM( variable ) )
675!
676!--    Note, to save memory arrays for heating are allocated from 0:pch_index.
677!--    Thus, output must be relative to these array indices. Further, check
678!--    whether the output is within the vertical output range,
679!--    i.e. nzb_do:nzt_do, which is necessary as local_pf is only allocated
680!--    for this index space. Note, plant-canopy output has a separate
681!--    vertical output coordinate zlad, so that output is mapped down to the
682!--    surface.
683       CASE ( 'pcm_heatrate' )
684          IF ( av == 0 )  THEN
685             DO  i = nxl, nxr
686                DO  j = nys, nyn
687                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
688                      local_pf(i,j,k) = pcm_heating_rate(k,j,i)
689                   ENDDO
690                ENDDO
691             ENDDO
692          ELSE
693             DO  i = nxl, nxr
694                DO  j = nys, nyn
695                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
696                      local_pf(i,j,k) = pcm_heatrate_av(k,j,i)
697                   ENDDO
698                ENDDO
699             ENDDO
700          ENDIF
701
702       CASE ( 'pcm_latentrate' )
703          IF ( av == 0 )  THEN
704             DO  i = nxl, nxr
705                DO  j = nys, nyn
706                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
707                      local_pf(i,j,k) = pcm_latent_rate(k,j,i)
708                   ENDDO
709                ENDDO
710             ENDDO
711          ELSE
712             DO  i = nxl, nxr
713                DO  j = nys, nyn
714                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
715                      local_pf(i,j,k) = pcm_latentrate_av(k,j,i)
716                   ENDDO
717                ENDDO
718             ENDDO
719          ENDIF
720
721       CASE ( 'pcm_transpirationrate' )
722          IF ( av == 0 )  THEN
723             DO  i = nxl, nxr
724                DO  j = nys, nyn
725                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
726                      local_pf(i,j,k) = pcm_transpiration_rate(k,j,i)
727                   ENDDO
728                ENDDO
729             ENDDO
730          ELSE
731             DO  i = nxl, nxr
732                DO  j = nys, nyn
733                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
734                      local_pf(i,j,k) = pcm_transpirationrate_av(k,j,i)
735                   ENDDO
736                ENDDO
737             ENDDO
738          ENDIF
739
740       CASE ( 'pcm_lad' )
741          IF ( av == 0 )  THEN
742             DO  i = nxl, nxr
743                DO  j = nys, nyn
744                   DO  k = MAX( 1, nzb_do ), MIN( pch_index, nzt_do )
745                      local_pf(i,j,k) = lad_s(k,j,i)
746                   ENDDO
747                ENDDO
748             ENDDO
749          ENDIF
750
751       CASE DEFAULT
752          found = .FALSE.
753
754    END SELECT
755
756 END SUBROUTINE pcm_data_output_3d 
757         
758!------------------------------------------------------------------------------!
759!
760! Description:
761! ------------
762!> Subroutine defining appropriate grid for netcdf variables.
763!> It is called from subroutine netcdf.
764!------------------------------------------------------------------------------!
765 SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
766
767     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
768     LOGICAL, INTENT(OUT)           ::  found       !<
769     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
770     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
771     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
772
773     found  = .TRUE.
774
775!
776!--  Check for the grid. zpc is zu(nzb:nzb+pch_index)
777     SELECT CASE ( TRIM( var ) )
778
779        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
780           grid_x = 'x'
781           grid_y = 'y'
782           grid_z = 'zpc'
783
784        CASE DEFAULT
785           found  = .FALSE.
786           grid_x = 'none'
787           grid_y = 'none'
788           grid_z = 'none'
789     END SELECT
790
791 END SUBROUTINE pcm_define_netcdf_grid
792 
793 
794!------------------------------------------------------------------------------!
795! Description:
796! ------------
797!> Header output for plant canopy model
798!------------------------------------------------------------------------------!
799    SUBROUTINE pcm_header ( io )
800 
801       CHARACTER (LEN=10) ::  coor_chr            !<
802
803       CHARACTER (LEN=86) ::  coordinates         !<
804       CHARACTER (LEN=86) ::  gradients           !<
805       CHARACTER (LEN=86) ::  leaf_area_density   !<
806       CHARACTER (LEN=86) ::  slices              !<
807 
808       INTEGER(iwp) :: i                !<
809       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
810       INTEGER(iwp) :: k                !<       
811   
812       REAL(wp) ::  canopy_height       !< canopy height (in m)
813       
814       canopy_height = zw(pch_index)
815
816       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
817                          canopy_drag_coeff
818       IF ( passive_scalar )  THEN
819          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
820                             leaf_surface_conc
821       ENDIF
822
823!
824!--    Heat flux at the top of vegetation
825       WRITE ( io, 3 )  cthf
826
827!
828!--    Leaf area density profile, calculated either from given vertical
829!--    gradients or from beta probability density function.
830       IF (  .NOT.  calc_beta_lad_profile )  THEN
831
832!--       Building output strings, starting with surface value
833          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
834          gradients = '------'
835          slices = '     0'
836          coordinates = '   0.0'
837          DO i = 1, UBOUND(lad_vertical_gradient_level_ind, DIM=1)
838             IF  ( lad_vertical_gradient_level_ind(i) /= -9999 ) THEN
839
840                WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
841                leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
842
843                WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
844                gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
845
846                WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
847                slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
848
849                WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
850                coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
851             ELSE
852                EXIT
853             ENDIF
854          ENDDO
855
856          WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
857                             TRIM( gradients ), TRIM( slices )
858
859       ELSE
860       
861          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
862          coordinates = '   0.0'
863         
864          DO  k = 1, pch_index
865
866             WRITE (coor_chr,'(F7.2)')  lad(k)
867             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
868                                 TRIM( coor_chr )
869 
870             WRITE (coor_chr,'(F7.1)')  zu(k)
871             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
872
873          ENDDO       
874
875          WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
876                          alpha_lad, beta_lad, lai_beta
877
878       ENDIF 
879
8801 FORMAT (//' Vegetation canopy (drag) model:'/                                &
881              ' ------------------------------'//                              &
882              ' Canopy mode: ', A /                                            &
883              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
884              ' Leaf drag coefficient: ',F6.2 /)
8852 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
886              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
8873 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
888          ' K m/s')
8894 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
890              ' Height:              ',A,'  m'/                                &
891              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
892              ' Gradient:            ',A,'  m**2/m**4'/                        &
893              ' Gridpoint:           ',A)
8945 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
895          //  ' Height:              ',A,'  m'/                                &
896              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
897              ' Coefficient alpha: ',F6.2 /                                    &
898              ' Coefficient beta: ',F6.2 /                                     &
899              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
900       
901    END SUBROUTINE pcm_header
902 
903 
904!------------------------------------------------------------------------------!
905! Description:
906! ------------
907!> Initialization of the plant canopy model
908!------------------------------------------------------------------------------!
909    SUBROUTINE pcm_init
910
911       INTEGER(iwp) ::  i   !< running index
912       INTEGER(iwp) ::  j   !< running index
913       INTEGER(iwp) ::  k   !< running index
914       INTEGER(iwp) ::  m   !< running index
915
916       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
917       REAL(wp) ::  gradient        !< gradient for lad-profile construction
918       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction     
919       REAL(wp) ::  lad_max         !< maximum LAD value in the model domain, used to perform a check
920
921       IF ( debug_output )  CALL debug_message( 'pcm_init', 'start' )
922!
923!--    Allocate one-dimensional arrays for the computation of the
924!--    leaf area density (lad) profile
925       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
926       lad = 0.0_wp
927       pre_lad = 0.0_wp
928
929!
930!--    Set flag that indicates that the lad-profile shall be calculated by using
931!--    a beta probability density function
932       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
933          calc_beta_lad_profile = .TRUE.
934       ENDIF
935       
936       
937!
938!--    Compute the profile of leaf area density used in the plant
939!--    canopy model. The profile can either be constructed from
940!--    prescribed vertical gradients of the leaf area density or by
941!--    using a beta probability density function (see e.g. Markkanen et al.,
942!--    2003: Boundary-Layer Meteorology, 106, 437-459)
943       IF (  .NOT.  calc_beta_lad_profile )  THEN   
944
945!
946!--       Use vertical gradients for lad-profile construction   
947          i = 1
948          gradient = 0.0_wp
949
950          lad(0) = lad_surface
951          lad_vertical_gradient_level_ind(1) = 0
952
953          DO k = 1, pch_index
954             IF ( i < 11 )  THEN
955                IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
956                     lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
957                   gradient = lad_vertical_gradient(i)
958                   lad_vertical_gradient_level_ind(i) = k - 1
959                   i = i + 1
960                ENDIF
961             ENDIF
962             IF ( gradient /= 0.0_wp )  THEN
963                IF ( k /= 1 )  THEN
964                   lad(k) = lad(k-1) + dzu(k) * gradient
965                ELSE
966                   lad(k) = lad_surface + dzu(k) * gradient
967                ENDIF
968             ELSE
969                lad(k) = lad(k-1)
970             ENDIF
971          ENDDO
972
973!
974!--       In case of no given leaf area density gradients, choose a vanishing
975!--       gradient. This information is used for the HEADER and the RUN_CONTROL
976!--       file.
977          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
978             lad_vertical_gradient_level(1) = 0.0_wp
979          ENDIF
980
981       ELSE
982
983!
984!--       Use beta function for lad-profile construction
985          int_bpdf = 0.0_wp
986          canopy_height = zw(pch_index)
987
988          DO k = 0, pch_index
989             int_bpdf = int_bpdf +                                             &
990                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
991                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
992                          beta_lad-1.0_wp ) )                                  &
993                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
994          ENDDO
995
996!
997!--       Preliminary lad profile (defined on w-grid)
998          DO k = 0, pch_index
999             pre_lad(k) =  lai_beta *                                          &
1000                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
1001                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
1002                              beta_lad-1.0_wp ) ) / int_bpdf                   &
1003                        ) / canopy_height
1004          ENDDO
1005
1006!
1007!--       Final lad profile (defined on scalar-grid level, since most prognostic
1008!--       quantities are defined there, hence, less interpolation is required
1009!--       when calculating the canopy tendencies)
1010          lad(0) = pre_lad(0)
1011          DO k = 1, pch_index
1012             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
1013          ENDDO
1014
1015       ENDIF
1016
1017!
1018!--    Allocate 3D-array for the leaf area density (lad_s).
1019       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1020
1021!
1022!--    Initialization of the canopy coverage in the model domain:
1023!--    Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which
1024!--    fully covers the domain surface
1025       SELECT CASE ( TRIM( canopy_mode ) )
1026
1027          CASE( 'homogeneous' )
1028
1029             DO  i = nxlg, nxrg
1030                DO  j = nysg, nyng
1031                   lad_s(:,j,i) = lad(:)
1032                ENDDO
1033             ENDDO
1034
1035          CASE ( 'read_from_file' )
1036!
1037!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
1038!--          use these values, else take LAD profiles from ASCII file.
1039!--          Please note, in NetCDF file LAD is only given up to the maximum
1040!--          canopy top, indicated by leaf_area_density_f%nz. 
1041             lad_s = 0.0_wp
1042             IF ( leaf_area_density_f%from_file )  THEN
1043!
1044!--             Set also pch_index, used to be the upper bound of the vertical
1045!--             loops. Therefore, use the global top of the canopy layer.
1046                pch_index = leaf_area_density_f%nz - 1
1047
1048                DO  i = nxl, nxr
1049                   DO  j = nys, nyn
1050                      DO  k = 0, leaf_area_density_f%nz - 1
1051                         IF ( leaf_area_density_f%var(k,j,i) /=                &
1052                              leaf_area_density_f%fill )                       &
1053                         lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
1054                      ENDDO
1055!
1056!--                   Check if resolved vegetation is mapped onto buildings.
1057!--                   In general, this is allowed and also meaningful, e.g.
1058!--                   when trees carry across roofs. However, due to the
1059!--                   topography filtering, new building grid points can emerge
1060!--                   at location where also plant canopy is defined. As a
1061!--                   result, plant canopy is mapped on top of roofs, with
1062!--                   siginficant impact on the downstream flow field and the
1063!--                   nearby surface radiation. In order to avoid that
1064!--                   plant canopy is mistakenly mapped onto building roofs,
1065!--                   check for building grid points (bit 6) that emerge from
1066!--                   the filtering (bit 4) and set LAD to zero at these
1067!--                   artificially  created building grid points. This case,
1068!--                   an informative message is given.
1069                      IF ( ANY( lad_s(:,j,i) /= 0.0_wp )  .AND.                &
1070                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
1071                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
1072                         lad_s(:,j,i) = 0.0_wp
1073                         WRITE( message_string, * )                            &
1074                                        'Resolved plant-canopy is ' //         &
1075                                        'defined on top of an artificially '// &
1076                                        'created building grid point ' //      &
1077                                        '(emerged from the filtering) - ' //   &
1078                                        'LAD profile is omitted at this ' //   &
1079                                        'grid point: (i,j) = ', i, j
1080                         CALL message( 'pcm_init', 'PA0313', 0, 0, myid, 6, 0 )
1081                      ENDIF
1082                   ENDDO
1083                ENDDO
1084                CALL exchange_horiz( lad_s, nbgp )
1085!
1086!            ASCII file
1087!--          Initialize canopy parameters canopy_drag_coeff,
1088!--          leaf_scalar_exch_coeff, leaf_surface_conc
1089!--          from file which contains complete 3D data (separate vertical profiles for
1090!--          each location).
1091             ELSE
1092                CALL pcm_read_plant_canopy_3d
1093             ENDIF
1094
1095          CASE DEFAULT
1096!
1097!--          The DEFAULT case is reached either if the parameter
1098!--          canopy mode contains a wrong character string or if the
1099!--          user has coded a special case in the user interface.
1100!--          There, the subroutine user_init_plant_canopy checks
1101!--          which of these two conditions applies.
1102             CALL user_init_plant_canopy
1103 
1104       END SELECT
1105!
1106!--    Check that at least one grid point has an LAD /= 0, else this may
1107!--    cause errors in the radiation model.
1108       lad_max = MAXVAL( lad_s )
1109#if defined( __parallel )
1110       CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX,        &
1111                           comm2d, ierr)
1112#endif
1113       IF ( lad_max <= 0.0_wp )  THEN
1114          message_string = 'Plant-canopy model is switched-on but no ' //      &
1115                           'plant canopy is present in the model domain.'
1116          CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
1117       ENDIF
1118   
1119!
1120!--    Initialize 2D index array indicating canopy top index.
1121       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
1122       pch_index_ji = 0
1123       
1124       DO  i = nxl, nxr
1125          DO  j = nys, nyn
1126             DO  k = 0, pch_index
1127                IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
1128             ENDDO
1129!
1130!--          Check whether topography and local vegetation on top exceed
1131!--          height of the model domain.
1132             IF ( topo_top_ind(j,i,0) + pch_index_ji(j,i) >= nzt + 1 )  THEN
1133                message_string =  'Local vegetation height on top of ' //      &
1134                                  'topography exceeds height of model domain.'
1135                CALL message( 'pcm_init', 'PA0674', 2, 2, myid, 6, 0 )
1136             ENDIF
1137
1138          ENDDO
1139       ENDDO
1140
1141       CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp )
1142!
1143!--    Calculate global pch_index value (index of top of plant canopy from ground)
1144       pch_index = MAXVAL( pch_index_ji )
1145       
1146       
1147!
1148!--    Exchange pch_index from all processors
1149#if defined( __parallel )
1150       CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER,            &
1151                           MPI_MAX, comm2d, ierr)
1152#endif
1153
1154!--    Allocation of arrays pcm_heating_rate, pcm_transpiration_rate and pcm_latent_rate
1155       ALLOCATE( pcm_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1156       pcm_heating_rate = 0.0_wp
1157       
1158       IF ( humidity )  THEN
1159          ALLOCATE( pcm_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1160          pcm_transpiration_rate = 0.0_wp
1161          ALLOCATE( pcm_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1162          pcm_latent_rate = 0.0_wp
1163       ENDIF
1164!
1165!--    Initialization of the canopy heat source distribution due to heating
1166!--    of the canopy layers by incoming solar radiation, in case that a non-zero
1167!--    value is set for the canopy top heat flux (cthf), which equals the
1168!--    available net radiation at canopy top.
1169!--    The heat source distribution is calculated by a decaying exponential
1170!--    function of the downward cumulative leaf area index (cum_lai_hf),
1171!--    assuming that the foliage inside the plant canopy is heated by solar
1172!--    radiation penetrating the canopy layers according to the distribution
1173!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
1174!--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
1175!--    Bound.-Layer Meteorol. 61, 47–64).
1176!--    When using the radiation_interactions, canopy heating (pcm_heating_rate)
1177!--    and plant canopy transpiration (pcm_transpiration_rate, pcm_latent_rate)
1178!--    are calculated in the RTM after the calculation of radiation.
1179       IF ( cthf /= 0.0_wp )  THEN
1180
1181          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1182!
1183!--       Piecewise calculation of the cumulative leaf area index by vertical
1184!--       integration of the leaf area density
1185          cum_lai_hf(:,:,:) = 0.0_wp
1186          DO  i = nxlg, nxrg
1187             DO  j = nysg, nyng
1188                DO  k = pch_index_ji(j,i)-1, 0, -1
1189                   IF ( k == pch_index_ji(j,i)-1 )  THEN
1190                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1191                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
1192                           ( zw(k+1) - zu(k+1) ) )  +                          &
1193                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1194                                                 lad_s(k,j,i) ) +              &
1195                                      lad_s(k+1,j,i) ) *                       &
1196                           ( zu(k+1) - zw(k) ) ) 
1197                   ELSE
1198                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1199                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
1200                                                 lad_s(k+1,j,i) ) +            &
1201                                      lad_s(k+1,j,i) ) *                       &
1202                           ( zw(k+1) - zu(k+1) ) )  +                          &
1203                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1204                                                 lad_s(k,j,i) ) +              &
1205                                      lad_s(k+1,j,i) ) *                       &
1206                           ( zu(k+1) - zw(k) ) )
1207                   ENDIF
1208                ENDDO
1209             ENDDO
1210          ENDDO
1211
1212!
1213!--       In areas with canopy the surface value of the canopy heat
1214!--       flux distribution overrides the surface heat flux (shf)
1215!--       Start with default surface type
1216          DO  m = 1, surf_def_h(0)%ns
1217             i = surf_def_h(0)%i(m)
1218             j = surf_def_h(0)%j(m)
1219             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1220                surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1221          ENDDO
1222!
1223!--       Natural surfaces
1224          DO  m = 1, surf_lsm_h%ns
1225             i = surf_lsm_h%i(m)
1226             j = surf_lsm_h%j(m)
1227             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1228                surf_lsm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1229          ENDDO
1230!
1231!--       Urban surfaces
1232          DO  m = 1, surf_usm_h%ns
1233             i = surf_usm_h%i(m)
1234             j = surf_usm_h%j(m)
1235             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1236                surf_usm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1237          ENDDO
1238!
1239!
1240!--       Calculation of the heating rate (K/s) within the different layers of
1241!--       the plant canopy. Calculation is only necessary in areas covered with
1242!--       canopy.
1243!--       Within the different canopy layers the plant-canopy heating
1244!--       rate (pcm_heating_rate) is calculated as the vertical
1245!--       divergence of the canopy heat fluxes at the top and bottom
1246!--       of the respective layer
1247          DO  i = nxlg, nxrg
1248             DO  j = nysg, nyng
1249                DO  k = 1, pch_index_ji(j,i)
1250                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
1251                      pcm_heating_rate(k,j,i) = cthf *                          &
1252                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -           &
1253                                  exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1254                   ENDIF
1255                ENDDO
1256             ENDDO
1257          ENDDO
1258
1259       ENDIF
1260
1261       IF ( debug_output )  CALL debug_message( 'pcm_init', 'end' )
1262
1263    END SUBROUTINE pcm_init
1264
1265
1266!------------------------------------------------------------------------------!
1267! Description:
1268! ------------
1269!> Parin for &plant_canopy_parameters for plant canopy model
1270!------------------------------------------------------------------------------!
1271    SUBROUTINE pcm_parin
1272
1273       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1274       
1275       NAMELIST /plant_canopy_parameters/                                      &
1276                                  alpha_lad, beta_lad, canopy_drag_coeff,      &
1277                                  canopy_mode, cthf,                           &
1278                                  lad_surface, lad_type_coef,                  & 
1279                                  lad_vertical_gradient,                       &
1280                                  lad_vertical_gradient_level,                 &
1281                                  lai_beta,                                    &
1282                                  leaf_scalar_exch_coeff,                      &
1283                                  leaf_surface_conc, pch_index,                &
1284                                  plant_canopy_transpiration
1285
1286       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
1287                                  canopy_mode, cthf,                           &
1288                                  lad_surface, lad_type_coef,                  & 
1289                                  lad_vertical_gradient,                       &
1290                                  lad_vertical_gradient_level,                 &
1291                                  lai_beta,                                    &
1292                                  leaf_scalar_exch_coeff,                      &
1293                                  leaf_surface_conc, pch_index,                &
1294                                  plant_canopy_transpiration
1295
1296       line = ' '
1297
1298!
1299!--    Try to find plant-canopy model package
1300       REWIND ( 11 )
1301       line = ' '
1302       DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
1303          READ ( 11, '(A)', END=12 )  line
1304       ENDDO
1305       BACKSPACE ( 11 )
1306
1307!
1308!--    Read user-defined namelist
1309       READ ( 11, plant_canopy_parameters, ERR = 10 )
1310
1311!
1312!--    Set flag that indicates that the plant-canopy model is switched on
1313       plant_canopy = .TRUE.
1314
1315       GOTO 14
1316
1317 10    BACKSPACE( 11 )
1318       READ( 11 , '(A)') line
1319       CALL parin_fail_message( 'plant_canopy_parameters', line )
1320!
1321!--    Try to find old namelist
1322 12    REWIND ( 11 )
1323       line = ' '
1324       DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
1325          READ ( 11, '(A)', END=14 )  line
1326       ENDDO
1327       BACKSPACE ( 11 )
1328
1329!
1330!--    Read user-defined namelist
1331       READ ( 11, canopy_par, ERR = 13, END = 14 )
1332
1333       message_string = 'namelist canopy_par is deprecated and will be ' // &
1334                     'removed in near future. Please use namelist ' //      &
1335                     'plant_canopy_parameters instead'
1336       CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
1337
1338!
1339!--    Set flag that indicates that the plant-canopy model is switched on
1340       plant_canopy = .TRUE.
1341
1342       GOTO 14
1343
1344 13    BACKSPACE( 11 )
1345       READ( 11 , '(A)') line
1346       CALL parin_fail_message( 'canopy_par', line )
1347
1348 14    CONTINUE
1349
1350    END SUBROUTINE pcm_parin
1351
1352
1353!------------------------------------------------------------------------------!
1354! Description:
1355! ------------
1356!
1357!> Loads 3D plant canopy data from file. File format is as follows:
1358!>
1359!> num_levels
1360!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1361!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1362!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1363!> ...
1364!>
1365!> i.e. first line determines number of levels and further lines represent plant
1366!> canopy data, one line per column and variable. In each data line,
1367!> dtype represents variable to be set:
1368!>
1369!> dtype=1: leaf area density (lad_s)
1370!> dtype=2....n: some additional plant canopy input data quantity
1371!>
1372!> Zeros are added automatically above num_levels until top of domain.  Any
1373!> non-specified (x,y) columns have zero values as default.
1374!------------------------------------------------------------------------------!
1375    SUBROUTINE pcm_read_plant_canopy_3d
1376
1377       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
1378       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
1379       INTEGER(iwp)                        ::  i, j      !< running index
1380       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
1381       INTEGER(iwp)                        ::  nzpltop   !<
1382       INTEGER(iwp)                        ::  nzpl      !<
1383       INTEGER(iwp)                        ::  kk        !<
1384       
1385       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
1386
1387!
1388!--    Initialize lad_s array
1389       lad_s = 0.0_wp
1390       
1391!
1392!--    Open and read plant canopy input data
1393       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
1394                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
1395                 FORM='FORMATTED', ERR=515)
1396       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
1397       nzpltop = MIN(nzt+1, nzb+nzp-1)
1398       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
1399       ALLOCATE( col(0:nzp-1) )
1400
1401       DO
1402          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
1403          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
1404         
1405          SELECT CASE (dtype)
1406             CASE( 1 )   !< leaf area density
1407!
1408!--             This is just the pure canopy layer assumed to be grounded to
1409!--             a flat domain surface. At locations where plant canopy sits
1410!--             on top of any kind of topography, the vertical plant column
1411!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
1412                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
1413                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
1414                                              'Allowed values 0 <= pctype <= 10, ' //  &
1415                                              'but pctype is ', pctype
1416                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1417                ENDIF
1418                kk = topo_top_ind(j,i,0)
1419                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
1420             CASE DEFAULT
1421                WRITE(message_string, '(a,i2,a)')   &
1422                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1423                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1424          END SELECT
1425       ENDDO
1426
1427515    message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1428       CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
1429
1430516    message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1431       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1432
1433517    CLOSE(152)
1434       DEALLOCATE( col )
1435       
1436       CALL exchange_horiz( lad_s, nbgp )
1437       
1438    END SUBROUTINE pcm_read_plant_canopy_3d
1439   
1440   
1441!------------------------------------------------------------------------------!
1442! Description:
1443! ------------
1444!> Calculation of the tendency terms, accounting for the effect of the plant
1445!> canopy on momentum and scalar quantities.
1446!>
1447!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1448!> (defined on scalar grid), as initialized in subroutine pcm_init.
1449!> The lad on the w-grid is vertically interpolated from the surrounding
1450!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1451!> k = pch_index. Here, the lad is zero.
1452!>
1453!> The canopy drag must be limited (previously accounted for by calculation of
1454!> a limiting canopy timestep for the determination of the maximum LES timestep
1455!> in subroutine timestep), since it is physically impossible that the canopy
1456!> drag alone can locally change the sign of a velocity component. This
1457!> limitation is realized by calculating preliminary tendencies and velocities.
1458!> It is subsequently checked if the preliminary new velocity has a different
1459!> sign than the current velocity. If so, the tendency is limited in a way that
1460!> the velocity can at maximum be reduced to zero by the canopy drag.
1461!>
1462!>
1463!> Call for all grid points
1464!------------------------------------------------------------------------------!
1465    SUBROUTINE pcm_tendency( component )
1466
1467       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1468       INTEGER(iwp) ::  i         !< running index
1469       INTEGER(iwp) ::  j         !< running index
1470       INTEGER(iwp) ::  k         !< running index
1471       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1472
1473       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
1474       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
1475       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
1476       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
1477
1478       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1479       REAL(wp) ::  lad_local !< local lad value
1480       REAL(wp) ::  pre_tend  !< preliminary tendency
1481       REAL(wp) ::  pre_u     !< preliminary u-value
1482       REAL(wp) ::  pre_v     !< preliminary v-value
1483       REAL(wp) ::  pre_w     !< preliminary w-value
1484
1485
1486       ddt_3d = 1.0_wp / dt_3d
1487 
1488!
1489!--    Compute drag for the three velocity components and the SGS-TKE:
1490       SELECT CASE ( component )
1491
1492!
1493!--       u-component
1494          CASE ( 1 )
1495             DO  i = nxlu, nxr
1496                DO  j = nys, nyn
1497!
1498!--                Set control flags indicating east- and westward-orientated
1499!--                building edges. Note, building_egde_w is set from the perspective
1500!--                of the potential rooftop grid point, while building_edge_e is
1501!--                is set from the perspective of the non-building grid point.
1502                   building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
1503                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
1504                   building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )&
1505                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
1506!
1507!--                Determine topography-top index on u-grid
1508                   DO  k = topo_top_ind(j,i,1)+1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
1509
1510                      kk = k - topo_top_ind(j,i,1)   !- lad arrays are defined flat
1511!
1512!--                   In order to create sharp boundaries of the plant canopy,
1513!--                   the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1514!--                   rather than being interpolated from the surrounding lad_s,
1515!--                   because this would yield smaller lad at the canopy boundaries
1516!--                   than inside of the canopy.
1517!--                   For the same reason, the lad at the rightmost(i+1)canopy
1518!--                   boundary on the u-grid equals lad_s(k,j,i), which is considered
1519!--                   in the next if-statement. Note, at left-sided building edges
1520!--                   this is not applied, here the LAD is equals the LAD at grid
1521!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
1522!--                   on top of a roof where (usually) is no LAD is defined.
1523                      lad_local = lad_s(kk,j,i)
1524                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  &
1525                           .AND.  .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
1526!
1527!--                   In order to avoid that LAD is mistakenly considered at right-
1528!--                   sided building edges (here the topography-top index for the
1529!--                   u-component at index j,i is still on the building while the
1530!--                   topography top for the scalar isn't), LAD is taken from grid
1531!--                   point (j,i-1).
1532                      IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  &
1533                           .AND.  building_edge_e )  lad_local = lad_s(kk,j,i-1)
1534
1535                      pre_tend = 0.0_wp
1536                      pre_u = 0.0_wp
1537!
1538!--                   Calculate preliminary value (pre_tend) of the tendency
1539                      pre_tend = - canopy_drag_coeff *                         &
1540                                   lad_local *                                 &
1541                                   SQRT( u(k,j,i)**2 +                         &
1542                                         ( 0.25_wp * ( v(k,j,i-1) +            &
1543                                                       v(k,j,i)   +            &
1544                                                       v(k,j+1,i) +            &
1545                                                       v(k,j+1,i-1) )          &
1546                                         )**2 +                                &
1547                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
1548                                                       w(k-1,j,i)   +          &
1549                                                       w(k,j,i-1)   +          &
1550                                                       w(k,j,i) )              &
1551                                         )**2                                  &
1552                                       ) *                                     &
1553                                   u(k,j,i)
1554
1555!
1556!--                   Calculate preliminary new velocity, based on pre_tend
1557                      pre_u = u(k,j,i) + dt_3d * pre_tend
1558!
1559!--                   Compare sign of old velocity and new preliminary velocity,
1560!--                   and in case the signs are different, limit the tendency
1561                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1562                         pre_tend = - u(k,j,i) * ddt_3d
1563                      ENDIF
1564!
1565!--                   Calculate final tendency
1566                      tend(k,j,i) = tend(k,j,i) + pre_tend
1567
1568                   ENDDO
1569                ENDDO
1570             ENDDO
1571
1572!
1573!--       v-component
1574          CASE ( 2 )
1575             DO  i = nxl, nxr
1576                DO  j = nysv, nyn
1577!
1578!--                Set control flags indicating north- and southward-orientated
1579!--                building edges. Note, building_egde_s is set from the perspective
1580!--                of the potential rooftop grid point, while building_edge_n is
1581!--                is set from the perspective of the non-building grid point.
1582                   building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
1583                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
1584                   building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )&
1585                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
1586!
1587!--                Determine topography-top index on v-grid
1588                   DO  k = topo_top_ind(j,i,2)+1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
1589
1590                      kk = k - topo_top_ind(j,i,2)   !- lad arrays are defined flat
1591!
1592!--                   In order to create sharp boundaries of the plant canopy,
1593!--                   the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1594!--                   rather than being interpolated from the surrounding lad_s,
1595!--                   because this would yield smaller lad at the canopy boundaries
1596!--                   than inside of the canopy.
1597!--                   For the same reason, the lad at the northmost(j+1)canopy
1598!--                   boundary on the v-grid equals lad_s(k,j,i), which is considered
1599!--                   in the next if-statement. Note, at left-sided building edges
1600!--                   this is not applied, here the LAD is equals the LAD at grid
1601!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
1602!--                   on top of a roof where (usually) is no LAD is defined.
1603                      lad_local = lad_s(kk,j,i)
1604                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
1605                       .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
1606!
1607!--                   In order to avoid that LAD is mistakenly considered at right-
1608!--                   sided building edges (here the topography-top index for the
1609!--                   u-component at index j,i is still on the building while the
1610!--                   topography top for the scalar isn't), LAD is taken from grid
1611!--                   point (j,i-1).
1612                      IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  &
1613                       .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
1614
1615                      pre_tend = 0.0_wp
1616                      pre_v = 0.0_wp
1617!
1618!--                   Calculate preliminary value (pre_tend) of the tendency
1619                      pre_tend = - canopy_drag_coeff *                         &
1620                                   lad_local *                                 &
1621                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
1622                                                       u(k,j-1,i+1) +          &
1623                                                       u(k,j,i)     +          &
1624                                                       u(k,j,i+1) )            &
1625                                         )**2 +                                &
1626                                         v(k,j,i)**2 +                         &
1627                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
1628                                                       w(k-1,j,i)   +          &
1629                                                       w(k,j-1,i)   +          &
1630                                                       w(k,j,i) )              &
1631                                         )**2                                  &
1632                                       ) *                                     &
1633                                   v(k,j,i)
1634
1635!
1636!--                   Calculate preliminary new velocity, based on pre_tend
1637                      pre_v = v(k,j,i) + dt_3d * pre_tend
1638!
1639!--                   Compare sign of old velocity and new preliminary velocity,
1640!--                   and in case the signs are different, limit the tendency
1641                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1642                         pre_tend = - v(k,j,i) * ddt_3d
1643                      ELSE
1644                         pre_tend = pre_tend
1645                      ENDIF
1646!
1647!--                   Calculate final tendency
1648                      tend(k,j,i) = tend(k,j,i) + pre_tend
1649
1650                   ENDDO
1651                ENDDO
1652             ENDDO
1653
1654!
1655!--       w-component
1656          CASE ( 3 )
1657             DO  i = nxl, nxr
1658                DO  j = nys, nyn
1659!
1660!--                Determine topography-top index on w-grid
1661                   DO  k = topo_top_ind(j,i,3)+1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
1662
1663                      kk = k - topo_top_ind(j,i,3)   !- lad arrays are defined flat
1664
1665                      pre_tend = 0.0_wp
1666                      pre_w = 0.0_wp
1667!
1668!--                   Calculate preliminary value (pre_tend) of the tendency
1669                      pre_tend = - canopy_drag_coeff *                         &
1670                                   (0.5_wp *                                   &
1671                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
1672                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
1673                                                       u(k,j,i+1) +            &
1674                                                       u(k+1,j,i) +            &
1675                                                       u(k+1,j,i+1) )          &
1676                                         )**2 +                                &
1677                                         ( 0.25_wp * ( v(k,j,i)   +            &
1678                                                       v(k,j+1,i) +            &
1679                                                       v(k+1,j,i) +            &
1680                                                       v(k+1,j+1,i) )          &
1681                                         )**2 +                                &
1682                                         w(k,j,i)**2                           &
1683                                       ) *                                     &
1684                                   w(k,j,i)
1685!
1686!--                   Calculate preliminary new velocity, based on pre_tend
1687                      pre_w = w(k,j,i) + dt_3d * pre_tend
1688!
1689!--                   Compare sign of old velocity and new preliminary velocity,
1690!--                   and in case the signs are different, limit the tendency
1691                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1692                         pre_tend = - w(k,j,i) * ddt_3d
1693                      ELSE
1694                         pre_tend = pre_tend
1695                      ENDIF
1696!
1697!--                   Calculate final tendency
1698                      tend(k,j,i) = tend(k,j,i) + pre_tend
1699
1700                   ENDDO
1701                ENDDO
1702             ENDDO
1703
1704!
1705!--       potential temperature
1706          CASE ( 4 )
1707             IF ( humidity ) THEN
1708                DO  i = nxl, nxr
1709                   DO  j = nys, nyn
1710!--                   Determine topography-top index on scalar-grid
1711                      DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
1712                         kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
1713                         tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
1714                      ENDDO
1715                   ENDDO
1716                ENDDO
1717             ELSE
1718                DO  i = nxl, nxr
1719                   DO  j = nys, nyn
1720!--                   Determine topography-top index on scalar-grid
1721                      DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
1722                         kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
1723                         tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
1724                      ENDDO
1725                   ENDDO
1726                ENDDO
1727             ENDIF
1728
1729!
1730!--       humidity
1731          CASE ( 5 )
1732             DO  i = nxl, nxr
1733                DO  j = nys, nyn
1734!
1735!--                Determine topography-top index on scalar-grid
1736                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
1737
1738                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
1739
1740                      IF ( .NOT. plant_canopy_transpiration ) THEN
1741                         ! pcm_transpiration_rate is calculated in radiation model
1742                         ! in case of plant_canopy_transpiration = .T.
1743                         ! to include also the dependecy to the radiation
1744                         ! in the plant canopy box
1745                         pcm_transpiration_rate(kk,j,i) =  - leaf_scalar_exch_coeff &
1746                                          * lad_s(kk,j,i) *                         &
1747                                          SQRT( ( 0.5_wp * ( u(k,j,i) +             &
1748                                                             u(k,j,i+1) )           &
1749                                                )**2 +                              &
1750                                                ( 0.5_wp * ( v(k,j,i) +             &
1751                                                             v(k,j+1,i) )           &
1752                                                )**2 +                              &
1753                                                ( 0.5_wp * ( w(k-1,j,i) +           &
1754                                                             w(k,j,i) )             &
1755                                                )**2                                &
1756                                              ) *                                   &
1757                                          ( q(k,j,i) - leaf_surface_conc )
1758                      ENDIF
1759
1760                      tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
1761                   ENDDO
1762                ENDDO
1763             ENDDO
1764
1765!
1766!--       sgs-tke
1767          CASE ( 6 )
1768             DO  i = nxl, nxr
1769                DO  j = nys, nyn
1770!
1771!--                Determine topography-top index on scalar-grid
1772                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
1773
1774                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
1775                      tend(k,j,i) = tend(k,j,i) -                              &
1776                                       2.0_wp * canopy_drag_coeff *            &
1777                                       lad_s(kk,j,i) *                         &
1778                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1779                                                          u(k,j,i+1) )         &
1780                                             )**2 +                            &
1781                                             ( 0.5_wp * ( v(k,j,i) +           &
1782                                                          v(k,j+1,i) )         &
1783                                             )**2 +                            &
1784                                             ( 0.5_wp * ( w(k,j,i) +           &
1785                                                          w(k+1,j,i) )         &
1786                                             )**2                              &
1787                                           ) *                                 &
1788                                       e(k,j,i)
1789                   ENDDO
1790                ENDDO
1791             ENDDO 
1792!
1793!--       scalar concentration
1794          CASE ( 7 )
1795             DO  i = nxl, nxr
1796                DO  j = nys, nyn
1797!
1798!--                Determine topography-top index on scalar-grid
1799                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
1800
1801                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
1802                      tend(k,j,i) = tend(k,j,i) -                              &
1803                                       leaf_scalar_exch_coeff *                &
1804                                       lad_s(kk,j,i) *                         &
1805                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1806                                                          u(k,j,i+1) )         &
1807                                             )**2 +                            &
1808                                             ( 0.5_wp * ( v(k,j,i) +           &
1809                                                          v(k,j+1,i) )         &
1810                                             )**2 +                            &
1811                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
1812                                                          w(k,j,i) )           &
1813                                             )**2                              &
1814                                           ) *                                 &
1815                                       ( s(k,j,i) - leaf_surface_conc )
1816                   ENDDO
1817                ENDDO
1818             ENDDO   
1819
1820
1821
1822          CASE DEFAULT
1823
1824             WRITE( message_string, * ) 'wrong component: ', component
1825             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1826
1827       END SELECT
1828
1829    END SUBROUTINE pcm_tendency
1830
1831
1832!------------------------------------------------------------------------------!
1833! Description:
1834! ------------
1835!> Calculation of the tendency terms, accounting for the effect of the plant
1836!> canopy on momentum and scalar quantities.
1837!>
1838!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1839!> (defined on scalar grid), as initialized in subroutine pcm_init.
1840!> The lad on the w-grid is vertically interpolated from the surrounding
1841!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1842!> k = pch_index. Here, the lad is zero.
1843!>
1844!> The canopy drag must be limited (previously accounted for by calculation of
1845!> a limiting canopy timestep for the determination of the maximum LES timestep
1846!> in subroutine timestep), since it is physically impossible that the canopy
1847!> drag alone can locally change the sign of a velocity component. This
1848!> limitation is realized by calculating preliminary tendencies and velocities.
1849!> It is subsequently checked if the preliminary new velocity has a different
1850!> sign than the current velocity. If so, the tendency is limited in a way that
1851!> the velocity can at maximum be reduced to zero by the canopy drag.
1852!>
1853!>
1854!> Call for grid point i,j
1855!------------------------------------------------------------------------------!
1856    SUBROUTINE pcm_tendency_ij( i, j, component )
1857
1858       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1859       INTEGER(iwp) ::  i         !< running index
1860       INTEGER(iwp) ::  j         !< running index
1861       INTEGER(iwp) ::  k         !< running index
1862       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1863
1864       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
1865       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
1866       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
1867       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
1868
1869       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1870       REAL(wp) ::  lad_local !< local lad value
1871       REAL(wp) ::  pre_tend  !< preliminary tendency
1872       REAL(wp) ::  pre_u     !< preliminary u-value
1873       REAL(wp) ::  pre_v     !< preliminary v-value
1874       REAL(wp) ::  pre_w     !< preliminary w-value
1875
1876
1877       ddt_3d = 1.0_wp / dt_3d
1878!
1879!--    Compute drag for the three velocity components and the SGS-TKE
1880       SELECT CASE ( component )
1881
1882!
1883!--       u-component
1884          CASE ( 1 )
1885!
1886!--          Set control flags indicating east- and westward-orientated
1887!--          building edges. Note, building_egde_w is set from the perspective
1888!--          of the potential rooftop grid point, while building_edge_e is
1889!--          is set from the perspective of the non-building grid point.
1890             building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
1891                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
1892             building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )  .AND. &
1893                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
1894!
1895!--          Determine topography-top index on u-grid
1896             DO  k = topo_top_ind(j,i,1) + 1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
1897
1898                kk = k - topo_top_ind(j,i,1)  !- lad arrays are defined flat
1899
1900!
1901!--             In order to create sharp boundaries of the plant canopy,
1902!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1903!--             rather than being interpolated from the surrounding lad_s,
1904!--             because this would yield smaller lad at the canopy boundaries
1905!--             than inside of the canopy.
1906!--             For the same reason, the lad at the rightmost(i+1)canopy
1907!--             boundary on the u-grid equals lad_s(k,j,i), which is considered
1908!--             in the next if-statement. Note, at left-sided building edges
1909!--             this is not applied, here the LAD is equals the LAD at grid
1910!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
1911!--             on top of a roof where (usually) is no LAD is defined.
1912                lad_local = lad_s(kk,j,i)
1913                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  .AND. &
1914                     .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
1915!
1916!--             In order to avoid that LAD is mistakenly considered at right-
1917!--             sided building edges (here the topography-top index for the
1918!--             u-component at index j,i is still on the building while the
1919!--             topography top for the scalar isn't), LAD is taken from grid
1920!--             point (j,i-1).
1921                IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  .AND. &
1922                     building_edge_e )  lad_local = lad_s(kk,j,i-1)
1923
1924                pre_tend = 0.0_wp
1925                pre_u = 0.0_wp
1926!
1927!--             Calculate preliminary value (pre_tend) of the tendency
1928                pre_tend = - canopy_drag_coeff *                               &
1929                             lad_local *                                       &   
1930                             SQRT( u(k,j,i)**2 +                               &
1931                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
1932                                                 v(k,j,i)    +                 &
1933                                                 v(k,j+1,i)  +                 &
1934                                                 v(k,j+1,i-1) )                &
1935                                   )**2 +                                      &
1936                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
1937                                                 w(k-1,j,i)   +                &
1938                                                 w(k,j,i-1)   +                &
1939                                                 w(k,j,i) )                    &
1940                                   )**2                                        &
1941                                 ) *                                           &
1942                             u(k,j,i)
1943
1944!
1945!--             Calculate preliminary new velocity, based on pre_tend
1946                pre_u = u(k,j,i) + dt_3d * pre_tend
1947!
1948!--             Compare sign of old velocity and new preliminary velocity,
1949!--             and in case the signs are different, limit the tendency
1950                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1951                   pre_tend = - u(k,j,i) * ddt_3d
1952                ELSE
1953                   pre_tend = pre_tend
1954                ENDIF
1955!
1956!--             Calculate final tendency
1957                tend(k,j,i) = tend(k,j,i) + pre_tend
1958             ENDDO
1959
1960
1961!
1962!--       v-component
1963          CASE ( 2 )
1964!
1965!--          Set control flags indicating north- and southward-orientated
1966!--          building edges. Note, building_egde_s is set from the perspective
1967!--          of the potential rooftop grid point, while building_edge_n is
1968!--          is set from the perspective of the non-building grid point.
1969             building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
1970                         .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
1971             building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )  .AND. &
1972                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
1973!
1974!--          Determine topography-top index on v-grid
1975             DO  k = topo_top_ind(j,i,2) + 1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
1976
1977                kk = k - topo_top_ind(j,i,2)  !- lad arrays are defined flat
1978!
1979!--             In order to create sharp boundaries of the plant canopy,
1980!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1981!--             rather than being interpolated from the surrounding lad_s,
1982!--             because this would yield smaller lad at the canopy boundaries
1983!--             than inside of the canopy.
1984!--             For the same reason, the lad at the northmost(j+1)canopy
1985!--             boundary on the v-grid equals lad_s(k,j,i), which is considered
1986!--             in the next if-statement. Note, at left-sided building edges
1987!--             this is not applied, here the LAD is equals the LAD at grid
1988!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
1989!--             on top of a roof where (usually) is no LAD is defined.
1990                lad_local = lad_s(kk,j,i)
1991                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp  .AND. &
1992                     .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
1993!
1994!--             In order to avoid that LAD is mistakenly considered at right-
1995!--             sided building edges (here the topography-top index for the
1996!--             u-component at index j,i is still on the building while the
1997!--             topography top for the scalar isn't), LAD is taken from grid
1998!--             point (j,i-1).
1999                IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  .AND. &
2000                     building_edge_n )  lad_local = lad_s(kk,j-1,i)
2001
2002                pre_tend = 0.0_wp
2003                pre_v = 0.0_wp
2004!
2005!--             Calculate preliminary value (pre_tend) of the tendency
2006                pre_tend = - canopy_drag_coeff *                               &
2007                             lad_local *                                       &
2008                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
2009                                                 u(k,j-1,i+1) +                &
2010                                                 u(k,j,i)     +                &
2011                                                 u(k,j,i+1) )                  &
2012                                   )**2 +                                      &
2013                                   v(k,j,i)**2 +                               &
2014                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
2015                                                 w(k-1,j,i)   +                &
2016                                                 w(k,j-1,i)   +                &
2017                                                 w(k,j,i) )                    &
2018                                   )**2                                        &
2019                                 ) *                                           &
2020                             v(k,j,i)
2021
2022!
2023!--             Calculate preliminary new velocity, based on pre_tend
2024                pre_v = v(k,j,i) + dt_3d * pre_tend
2025!
2026!--             Compare sign of old velocity and new preliminary velocity,
2027!--             and in case the signs are different, limit the tendency
2028                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
2029                   pre_tend = - v(k,j,i) * ddt_3d
2030                ELSE
2031                   pre_tend = pre_tend
2032                ENDIF
2033!
2034!--             Calculate final tendency
2035                tend(k,j,i) = tend(k,j,i) + pre_tend
2036             ENDDO
2037
2038
2039!
2040!--       w-component
2041          CASE ( 3 )
2042!
2043!--          Determine topography-top index on w-grid
2044             DO  k = topo_top_ind(j,i,3) + 1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
2045
2046                kk = k - topo_top_ind(j,i,3)  !- lad arrays are defined flat
2047
2048                pre_tend = 0.0_wp
2049                pre_w = 0.0_wp
2050!
2051!--             Calculate preliminary value (pre_tend) of the tendency
2052                pre_tend = - canopy_drag_coeff *                               &
2053                             (0.5_wp *                                         &
2054                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
2055                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
2056                                                 u(k,j,i+1)  +                 &
2057                                                 u(k+1,j,i)  +                 &
2058                                                 u(k+1,j,i+1) )                &
2059                                   )**2 +                                      &
2060                                   ( 0.25_wp * ( v(k,j,i)    +                 &
2061                                                 v(k,j+1,i)  +                 &
2062                                                 v(k+1,j,i)  +                 &
2063                                                 v(k+1,j+1,i) )                &
2064                                   )**2 +                                      &
2065                                   w(k,j,i)**2                                 &
2066                                 ) *                                           &
2067                             w(k,j,i)
2068!
2069!--             Calculate preliminary new velocity, based on pre_tend
2070                pre_w = w(k,j,i) + dt_3d * pre_tend
2071!
2072!--             Compare sign of old velocity and new preliminary velocity,
2073!--             and in case the signs are different, limit the tendency
2074                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
2075                   pre_tend = - w(k,j,i) * ddt_3d
2076                ELSE
2077                   pre_tend = pre_tend
2078                ENDIF
2079!
2080!--             Calculate final tendency
2081                tend(k,j,i) = tend(k,j,i) + pre_tend
2082             ENDDO
2083
2084!
2085!--       potential temperature
2086          CASE ( 4 )
2087!
2088!--          Determine topography-top index on scalar grid
2089             IF ( humidity ) THEN
2090                DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2091                   kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
2092                   tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) -       &
2093                                 pcm_latent_rate(kk,j,i)
2094                ENDDO
2095             ELSE
2096                DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2097                   kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
2098                   tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
2099                ENDDO
2100             ENDIF
2101
2102!
2103!--       humidity
2104          CASE ( 5 )
2105!
2106!--          Determine topography-top index on scalar grid
2107             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2108                kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
2109                IF ( .NOT. plant_canopy_transpiration ) THEN
2110                   ! pcm_transpiration_rate is calculated in radiation model
2111                   ! in case of plant_canopy_transpiration = .T.
2112                   ! to include also the dependecy to the radiation
2113                   ! in the plant canopy box
2114                   pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff   &
2115                                    * lad_s(kk,j,i) *                          &
2116                                    SQRT( ( 0.5_wp * ( u(k,j,i) +              &
2117                                                       u(k,j,i+1) )            &
2118                                          )**2  +                              &
2119                                          ( 0.5_wp * ( v(k,j,i) +              &
2120                                                       v(k,j+1,i) )            &
2121                                          )**2 +                               &
2122                                          ( 0.5_wp * ( w(k-1,j,i) +            &
2123                                                       w(k,j,i) )              &
2124                                          )**2                                 &
2125                                        ) *                                    &
2126                                    ( q(k,j,i) - leaf_surface_conc )
2127                ENDIF
2128
2129                tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
2130
2131             ENDDO   
2132
2133!
2134!--       sgs-tke
2135          CASE ( 6 )
2136!
2137!--          Determine topography-top index on scalar grid
2138             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2139
2140                kk = k - topo_top_ind(j,i,0)
2141                tend(k,j,i) = tend(k,j,i) -                                    &
2142                                 2.0_wp * canopy_drag_coeff *                  &
2143                                 lad_s(kk,j,i) *                               &
2144                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2145                                                    u(k,j,i+1) )               &
2146                                       )**2 +                                  & 
2147                                       ( 0.5_wp * ( v(k,j,i) +                 &
2148                                                    v(k,j+1,i) )               &
2149                                       )**2 +                                  &
2150                                       ( 0.5_wp * ( w(k,j,i) +                 &
2151                                                    w(k+1,j,i) )               &
2152                                       )**2                                    &
2153                                     ) *                                       &
2154                                 e(k,j,i)
2155             ENDDO
2156             
2157!
2158!--       scalar concentration
2159          CASE ( 7 )
2160!
2161!--          Determine topography-top index on scalar grid
2162             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2163
2164                kk = k - topo_top_ind(j,i,0)
2165                tend(k,j,i) = tend(k,j,i) -                                    &
2166                                 leaf_scalar_exch_coeff *                      &
2167                                 lad_s(kk,j,i) *                               &
2168                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2169                                                    u(k,j,i+1) )               &
2170                                       )**2  +                                 &
2171                                       ( 0.5_wp * ( v(k,j,i) +                 &
2172                                                    v(k,j+1,i) )               &
2173                                       )**2 +                                  &
2174                                       ( 0.5_wp * ( w(k-1,j,i) +               &
2175                                                    w(k,j,i) )                 &
2176                                       )**2                                    &
2177                                     ) *                                       &
2178                                 ( s(k,j,i) - leaf_surface_conc )
2179             ENDDO               
2180
2181       CASE DEFAULT
2182
2183          WRITE( message_string, * ) 'wrong component: ', component
2184          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
2185
2186       END SELECT
2187
2188    END SUBROUTINE pcm_tendency_ij
2189
2190
2191 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.