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

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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