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

Last change on this file since 4342 was 4342, checked in by Giersch, 16 months ago

last commit corrected by replacing read_from_file_3d with read_from_file in all relevant p3d files, minor changes in plant_canopy_model_mod (use statements moved to module level, ocean dependency removed, redundant variables removed)

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