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

Last change on this file since 4381 was 4381, checked in by suehring, 17 months ago

Land-surface model: Bugfix in nested soil initialization in case no dynamic input file is present; give local error messages only onces

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