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

Last change on this file since 4692 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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