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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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