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

Last change on this file since 4341 was 4341, checked in by motisi, 4 years ago

plant_canopy_model: unification of variable names: pc_ variabels now pcm_ variables
removal of bowenratio output
canopy mode 'block' now 'homogeneous'
'read_from_file_3d' now 'read_from_file'
removal of confusing comment lines
replacement of k_wall by topo_top_ind
removal of else-statement in tendency-calculation

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