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

Last change on this file since 4515 was 4515, checked in by suehring, 14 months ago

Rename error number again since this was recently given in -r 4511

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