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

Last change on this file since 4392 was 4392, checked in by pavelkrc, 20 months ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

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