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

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

Plant canopy: implement fix for LAD on building edges also for the vectorl-optimized branch

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