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

Last change on this file since 4221 was 4221, checked in by suehring, 2 years ago

Bugfix in plant-canopy model, missing intialization of heating rate

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