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

Last change on this file since 4355 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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