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

Last change on this file since 4329 was 4329, checked in by motisi, 5 years ago

Renamed wall_flags_0 to wall_flags_static_0

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