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

Last change on this file since 4808 was 4803, checked in by raasch, 4 years ago

file re-formatted to follow the PALM coding standard

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