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

Last change on this file since 4631 was 4535, checked in by raasch, 4 years ago

bugfix for restart data format query

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