source: palm/trunk/SOURCE/plant_canopy_model_mod.f90

Last change on this file was 4843, checked in by raasch, 4 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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