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

Last change on this file since 4361 was 4361, checked in by suehring, 15 months ago

Plant-canopy: Avoid an unnecessary exchange of ghost points; remove unused arrays in pcm_rrd_local

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