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

Last change on this file since 4302 was 4302, checked in by suehring, 17 months ago

Omit tall canopy mapped on top of buildings, which may happen due to topography filtering

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