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

Last change on this file since 4362 was 4362, checked in by suehring, 19 months ago

Input of plant-canopy variables from static driver moved from netcdf_data_input_mod to plant-canopy model

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