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

Last change on this file since 4331 was 4331, checked in by suehring, 16 months ago

New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output

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