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

Last change on this file since 4360 was 4360, checked in by suehring, 15 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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