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

Last change on this file since 4327 was 4314, checked in by suehring, 5 years ago

plant-canopy: Bugfix, plant canopy was still considered at building edges on for the u- and v-component; Relax restriction of LAD on building tops. LAD is only omitted at locations where building grid points emerged artificially by the topography filtering

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