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

Last change on this file since 4363 was 4363, checked in by suehring, 18 months ago

plant-canopy model: fix misplaced endif

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