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

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

Map 3d plant-canopy output to surface (this feature was lost by the latest bugfix in 3d canopy output)

  • Property svn:keywords set to Id
File size: 89.6 KB
Line 
1!> @file plant_canopy_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Institute of Computer Science of the
19!                     Czech Academy of Sciences, Prague
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: plant_canopy_model_mod.f90 4220 2019-09-05 16:42:39Z 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       IF ( humidity )  THEN
1153          ALLOCATE( pc_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1154          pc_transpiration_rate = 0.0_wp
1155          ALLOCATE( pc_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1156          pc_latent_rate = 0.0_wp
1157       ENDIF
1158!
1159!--    Initialization of the canopy heat source distribution due to heating
1160!--    of the canopy layers by incoming solar radiation, in case that a non-zero
1161!--    value is set for the canopy top heat flux (cthf), which equals the
1162!--    available net radiation at canopy top.
1163!--    The heat source distribution is calculated by a decaying exponential
1164!--    function of the downward cumulative leaf area index (cum_lai_hf),
1165!--    assuming that the foliage inside the plant canopy is heated by solar
1166!--    radiation penetrating the canopy layers according to the distribution
1167!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
1168!--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
1169!--    Bound.-Layer Meteorol. 61, 47–64).
1170!--    When using the radiation_interactions, canopy heating (pc_heating_rate)
1171!--    and plant canopy transpiration (pc_transpiration_rate, pc_latent_rate)
1172!--    are calculated in the RTM after the calculation of radiation.
1173!--    We cannot use variable radiation_interactions here to determine the situation
1174!--    as it is assigned in init_3d_model after the call of pcm_init.
1175       IF ( cthf /= 0.0_wp )  THEN
1176
1177          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1178!
1179!--       Piecewise calculation of the cumulative leaf area index by vertical
1180!--       integration of the leaf area density
1181          cum_lai_hf(:,:,:) = 0.0_wp
1182          DO  i = nxlg, nxrg
1183             DO  j = nysg, nyng
1184                DO  k = pch_index_ji(j,i)-1, 0, -1
1185                   IF ( k == pch_index_ji(j,i)-1 )  THEN
1186                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1187                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
1188                           ( zw(k+1) - zu(k+1) ) )  +                          &
1189                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1190                                                 lad_s(k,j,i) ) +              &
1191                                      lad_s(k+1,j,i) ) *                       &
1192                           ( zu(k+1) - zw(k) ) ) 
1193                   ELSE
1194                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1195                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
1196                                                 lad_s(k+1,j,i) ) +            &
1197                                      lad_s(k+1,j,i) ) *                       &
1198                           ( zw(k+1) - zu(k+1) ) )  +                          &
1199                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1200                                                 lad_s(k,j,i) ) +              &
1201                                      lad_s(k+1,j,i) ) *                       &
1202                           ( zu(k+1) - zw(k) ) )
1203                   ENDIF
1204                ENDDO
1205             ENDDO
1206          ENDDO
1207
1208!           
1209!--       In areas with canopy the surface value of the canopy heat
1210!--       flux distribution overrides the surface heat flux (shf)
1211!--       Start with default surface type
1212          DO  m = 1, surf_def_h(0)%ns
1213             k = surf_def_h(0)%k(m)
1214             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1215                surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1216          ENDDO
1217!
1218!--       Natural surfaces
1219          DO  m = 1, surf_lsm_h%ns
1220             k = surf_lsm_h%k(m)
1221             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1222                surf_lsm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1223          ENDDO
1224!
1225!--       Urban surfaces
1226          DO  m = 1, surf_usm_h%ns
1227             k = surf_usm_h%k(m)
1228             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1229                surf_usm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1230          ENDDO
1231!
1232!
1233!--       Calculation of the heating rate (K/s) within the different layers of
1234!--       the plant canopy. Calculation is only necessary in areas covered with
1235!--       canopy.
1236!--       Within the different canopy layers the plant-canopy heating
1237!--       rate (pc_heating_rate) is calculated as the vertical
1238!--       divergence of the canopy heat fluxes at the top and bottom
1239!--       of the respective layer
1240          DO  i = nxlg, nxrg
1241             DO  j = nysg, nyng
1242                DO  k = 1, pch_index_ji(j,i)
1243                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
1244                      pc_heating_rate(k,j,i) = cthf *                          &
1245                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -           &
1246                                  exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1247                   ENDIF
1248                ENDDO
1249             ENDDO
1250          ENDDO
1251
1252       ENDIF
1253
1254       IF ( debug_output )  CALL debug_message( 'pcm_init', 'end' )
1255
1256
1257    END SUBROUTINE pcm_init
1258
1259
1260!------------------------------------------------------------------------------!
1261! Description:
1262! ------------
1263!> Parin for &plant_canopy_parameters for plant canopy model
1264!------------------------------------------------------------------------------!
1265    SUBROUTINE pcm_parin
1266
1267       USE control_parameters,                                                 &
1268           ONLY:  message_string, plant_canopy
1269
1270       IMPLICIT NONE
1271
1272       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1273       
1274       NAMELIST /plant_canopy_parameters/                                      &
1275                                  alpha_lad, beta_lad, canopy_drag_coeff,      &
1276                                  canopy_mode, cthf,                           &
1277                                  lad_surface, lad_type_coef,                  & 
1278                                  lad_vertical_gradient,                       &
1279                                  lad_vertical_gradient_level,                 &
1280                                  lai_beta,                                    &
1281                                  leaf_scalar_exch_coeff,                      &
1282                                  leaf_surface_conc, pch_index,                &
1283                                  plant_canopy_transpiration
1284
1285       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
1286                                  canopy_mode, cthf,                           &
1287                                  lad_surface, lad_type_coef,                  & 
1288                                  lad_vertical_gradient,                       &
1289                                  lad_vertical_gradient_level,                 &
1290                                  lai_beta,                                    &
1291                                  leaf_scalar_exch_coeff,                      &
1292                                  leaf_surface_conc, pch_index,                &
1293                                  plant_canopy_transpiration
1294
1295       line = ' '
1296
1297!
1298!--    Try to find radiation model package
1299       REWIND ( 11 )
1300       line = ' '
1301       DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
1302          READ ( 11, '(A)', END=12 )  line
1303       ENDDO
1304       BACKSPACE ( 11 )
1305
1306!
1307!--    Read user-defined namelist
1308       READ ( 11, plant_canopy_parameters, ERR = 10 )
1309
1310!
1311!--    Set flag that indicates that the radiation model is switched on
1312       plant_canopy = .TRUE.
1313
1314       GOTO 14
1315
1316 10    BACKSPACE( 11 )
1317       READ( 11 , '(A)') line
1318       CALL parin_fail_message( 'plant_canopy_parameters', line )
1319!
1320!--    Try to find old namelist
1321 12    REWIND ( 11 )
1322       line = ' '
1323       DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
1324          READ ( 11, '(A)', END=14 )  line
1325       ENDDO
1326       BACKSPACE ( 11 )
1327
1328!
1329!--    Read user-defined namelist
1330       READ ( 11, canopy_par, ERR = 13, END = 14 )
1331
1332       message_string = 'namelist canopy_par is deprecated and will be ' // &
1333                     'removed in near future. Please use namelist ' //      &
1334                     'plant_canopy_parameters instead'
1335       CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
1336
1337!
1338!--    Set flag that indicates that the radiation model is switched on
1339       plant_canopy = .TRUE.
1340
1341       GOTO 14
1342
1343 13    BACKSPACE( 11 )
1344       READ( 11 , '(A)') line
1345       CALL parin_fail_message( 'canopy_par', line )
1346
1347 14    CONTINUE
1348
1349
1350    END SUBROUTINE pcm_parin
1351
1352
1353
1354!------------------------------------------------------------------------------!
1355! Description:
1356! ------------
1357!
1358!> Loads 3D plant canopy data from file. File format is as follows:
1359!>
1360!> num_levels
1361!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1362!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1363!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1364!> ...
1365!>
1366!> i.e. first line determines number of levels and further lines represent plant
1367!> canopy data, one line per column and variable. In each data line,
1368!> dtype represents variable to be set:
1369!>
1370!> dtype=1: leaf area density (lad_s)
1371!> dtype=2....n: some additional plant canopy input data quantity
1372!>
1373!> Zeros are added automatically above num_levels until top of domain.  Any
1374!> non-specified (x,y) columns have zero values as default.
1375!------------------------------------------------------------------------------!
1376    SUBROUTINE pcm_read_plant_canopy_3d
1377   
1378       USE control_parameters,                                                 &
1379           ONLY:  coupling_char, message_string
1380
1381       USE indices,                                                            &
1382           ONLY:  nbgp
1383           
1384       IMPLICIT NONE
1385
1386       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
1387       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
1388       INTEGER(iwp)                        ::  i, j      !< running index
1389       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
1390       INTEGER(iwp)                        ::  nzpltop   !<
1391       INTEGER(iwp)                        ::  nzpl      !<
1392       INTEGER(iwp)                        ::  kk        !<
1393       
1394       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
1395
1396!
1397!--    Initialize lad_s array
1398       lad_s = 0.0_wp
1399       
1400!
1401!--    Open and read plant canopy input data
1402       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
1403                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
1404                 FORM='FORMATTED', ERR=515)
1405       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
1406       nzpltop = MIN(nzt+1, nzb+nzp-1)
1407       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
1408       ALLOCATE( col(0:nzp-1) )
1409
1410       DO
1411          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
1412          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
1413         
1414          SELECT CASE (dtype)
1415             CASE( 1 )   !< leaf area density
1416!
1417!--             This is just the pure canopy layer assumed to be grounded to
1418!--             a flat domain surface. At locations where plant canopy sits
1419!--             on top of any kind of topography, the vertical plant column
1420!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
1421                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
1422                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
1423                                              'Allowed values 0 <= pctype <= 10, ' //  &
1424                                              'but pctype is ', pctype
1425                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1426                ENDIF
1427                kk = topo_top_ind(j,i,0)
1428                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
1429             CASE DEFAULT
1430                WRITE(message_string, '(a,i2,a)')   &
1431                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1432                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1433          END SELECT
1434       ENDDO
1435
1436515    message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1437       CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
1438
1439516    message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1440       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1441
1442517    CLOSE(152)
1443       DEALLOCATE( col )
1444       
1445       CALL exchange_horiz( lad_s, nbgp )
1446       
1447    END SUBROUTINE pcm_read_plant_canopy_3d
1448   
1449   
1450
1451!------------------------------------------------------------------------------!
1452! Description:
1453! ------------
1454!> Calculation of the tendency terms, accounting for the effect of the plant
1455!> canopy on momentum and scalar quantities.
1456!>
1457!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1458!> (defined on scalar grid), as initialized in subroutine pcm_init.
1459!> The lad on the w-grid is vertically interpolated from the surrounding
1460!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1461!> k = pch_index. Here, the lad is zero.
1462!>
1463!> The canopy drag must be limited (previously accounted for by calculation of
1464!> a limiting canopy timestep for the determination of the maximum LES timestep
1465!> in subroutine timestep), since it is physically impossible that the canopy
1466!> drag alone can locally change the sign of a velocity component. This
1467!> limitation is realized by calculating preliminary tendencies and velocities.
1468!> It is subsequently checked if the preliminary new velocity has a different
1469!> sign than the current velocity. If so, the tendency is limited in a way that
1470!> the velocity can at maximum be reduced to zero by the canopy drag.
1471!>
1472!>
1473!> Call for all grid points
1474!------------------------------------------------------------------------------!
1475    SUBROUTINE pcm_tendency( component )
1476
1477
1478       USE control_parameters,                                                 &
1479           ONLY:  dt_3d, message_string
1480
1481       USE kinds
1482
1483       IMPLICIT NONE
1484
1485       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1486       INTEGER(iwp) ::  i         !< running index
1487       INTEGER(iwp) ::  j         !< running index
1488       INTEGER(iwp) ::  k         !< running index
1489       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
1490       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1491
1492       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1493       REAL(wp) ::  lad_local !< local lad value
1494       REAL(wp) ::  pre_tend  !< preliminary tendency
1495       REAL(wp) ::  pre_u     !< preliminary u-value
1496       REAL(wp) ::  pre_v     !< preliminary v-value
1497       REAL(wp) ::  pre_w     !< preliminary w-value
1498
1499
1500       ddt_3d = 1.0_wp / dt_3d
1501 
1502!
1503!--    Compute drag for the three velocity components and the SGS-TKE:
1504       SELECT CASE ( component )
1505
1506!
1507!--       u-component
1508          CASE ( 1 )
1509             DO  i = nxlu, nxr
1510                DO  j = nys, nyn
1511!
1512!--                Determine topography-top index on u-grid
1513                   k_wall = topo_top_ind(j,i,1)
1514                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1515
1516                      kk = k - k_wall   !- lad arrays are defined flat
1517!
1518!--                   In order to create sharp boundaries of the plant canopy,
1519!--                   the lad on the u-grid at index (k,j,i) is equal to
1520!--                   lad_s(k,j,i), rather than being interpolated from the
1521!--                   surrounding lad_s, because this would yield smaller lad
1522!--                   at the canopy boundaries than inside of the canopy.
1523!--                   For the same reason, the lad at the rightmost(i+1)canopy
1524!--                   boundary on the u-grid equals lad_s(k,j,i).
1525                      lad_local = lad_s(kk,j,i)
1526                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
1527                      THEN
1528                         lad_local = lad_s(kk,j,i-1)
1529                      ENDIF
1530
1531                      pre_tend = 0.0_wp
1532                      pre_u = 0.0_wp
1533!
1534!--                   Calculate preliminary value (pre_tend) of the tendency
1535                      pre_tend = - cdc *                                       &
1536                                   lad_local *                                 &
1537                                   SQRT( u(k,j,i)**2 +                         &
1538                                         ( 0.25_wp * ( v(k,j,i-1) +            &
1539                                                       v(k,j,i)   +            &
1540                                                       v(k,j+1,i) +            &
1541                                                       v(k,j+1,i-1) )          &
1542                                         )**2 +                                &
1543                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
1544                                                       w(k-1,j,i)   +          &
1545                                                       w(k,j,i-1)   +          &
1546                                                       w(k,j,i) )              &
1547                                         )**2                                  &
1548                                       ) *                                     &
1549                                   u(k,j,i)
1550
1551!
1552!--                   Calculate preliminary new velocity, based on pre_tend
1553                      pre_u = u(k,j,i) + dt_3d * pre_tend
1554!
1555!--                   Compare sign of old velocity and new preliminary velocity,
1556!--                   and in case the signs are different, limit the tendency
1557                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1558                         pre_tend = - u(k,j,i) * ddt_3d
1559                      ELSE
1560                         pre_tend = pre_tend
1561                      ENDIF
1562!
1563!--                   Calculate final tendency
1564                      tend(k,j,i) = tend(k,j,i) + pre_tend
1565
1566                   ENDDO
1567                ENDDO
1568             ENDDO
1569
1570!
1571!--       v-component
1572          CASE ( 2 )
1573             DO  i = nxl, nxr
1574                DO  j = nysv, nyn
1575!
1576!--                Determine topography-top index on v-grid
1577                   k_wall = topo_top_ind(j,i,2)
1578
1579                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1580
1581                      kk = k - k_wall   !- lad arrays are defined flat
1582!
1583!--                   In order to create sharp boundaries of the plant canopy,
1584!--                   the lad on the v-grid at index (k,j,i) is equal to
1585!--                   lad_s(k,j,i), rather than being interpolated from the
1586!--                   surrounding lad_s, because this would yield smaller lad
1587!--                   at the canopy boundaries than inside of the canopy.
1588!--                   For the same reason, the lad at the northmost(j+1) canopy
1589!--                   boundary on the v-grid equals lad_s(k,j,i).
1590                      lad_local = lad_s(kk,j,i)
1591                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
1592                      THEN
1593                         lad_local = lad_s(kk,j-1,i)
1594                      ENDIF
1595
1596                      pre_tend = 0.0_wp
1597                      pre_v = 0.0_wp
1598!
1599!--                   Calculate preliminary value (pre_tend) of the tendency
1600                      pre_tend = - cdc *                                       &
1601                                   lad_local *                                 &
1602                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
1603                                                       u(k,j-1,i+1) +          &
1604                                                       u(k,j,i)     +          &
1605                                                       u(k,j,i+1) )            &
1606                                         )**2 +                                &
1607                                         v(k,j,i)**2 +                         &
1608                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
1609                                                       w(k-1,j,i)   +          &
1610                                                       w(k,j-1,i)   +          &
1611                                                       w(k,j,i) )              &
1612                                         )**2                                  &
1613                                       ) *                                     &
1614                                   v(k,j,i)
1615
1616!
1617!--                   Calculate preliminary new velocity, based on pre_tend
1618                      pre_v = v(k,j,i) + dt_3d * pre_tend
1619!
1620!--                   Compare sign of old velocity and new preliminary velocity,
1621!--                   and in case the signs are different, limit the tendency
1622                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1623                         pre_tend = - v(k,j,i) * ddt_3d
1624                      ELSE
1625                         pre_tend = pre_tend
1626                      ENDIF
1627!
1628!--                   Calculate final tendency
1629                      tend(k,j,i) = tend(k,j,i) + pre_tend
1630
1631                   ENDDO
1632                ENDDO
1633             ENDDO
1634
1635!
1636!--       w-component
1637          CASE ( 3 )
1638             DO  i = nxl, nxr
1639                DO  j = nys, nyn
1640!
1641!--                Determine topography-top index on w-grid
1642                   k_wall = topo_top_ind(j,i,3)
1643
1644                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
1645
1646                      kk = k - k_wall   !- lad arrays are defined flat
1647
1648                      pre_tend = 0.0_wp
1649                      pre_w = 0.0_wp
1650!
1651!--                   Calculate preliminary value (pre_tend) of the tendency
1652                      pre_tend = - cdc *                                       &
1653                                   (0.5_wp *                                   &
1654                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
1655                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
1656                                                       u(k,j,i+1) +            &
1657                                                       u(k+1,j,i) +            &
1658                                                       u(k+1,j,i+1) )          &
1659                                         )**2 +                                &
1660                                         ( 0.25_wp * ( v(k,j,i)   +            &
1661                                                       v(k,j+1,i) +            &
1662                                                       v(k+1,j,i) +            &
1663                                                       v(k+1,j+1,i) )          &
1664                                         )**2 +                                &
1665                                         w(k,j,i)**2                           &
1666                                       ) *                                     &
1667                                   w(k,j,i)
1668!
1669!--                   Calculate preliminary new velocity, based on pre_tend
1670                      pre_w = w(k,j,i) + dt_3d * pre_tend
1671!
1672!--                   Compare sign of old velocity and new preliminary velocity,
1673!--                   and in case the signs are different, limit the tendency
1674                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1675                         pre_tend = - w(k,j,i) * ddt_3d
1676                      ELSE
1677                         pre_tend = pre_tend
1678                      ENDIF
1679!
1680!--                   Calculate final tendency
1681                      tend(k,j,i) = tend(k,j,i) + pre_tend
1682
1683                   ENDDO
1684                ENDDO
1685             ENDDO
1686
1687!
1688!--       potential temperature
1689          CASE ( 4 )
1690             IF ( humidity ) THEN
1691                DO  i = nxl, nxr
1692                   DO  j = nys, nyn
1693!--                   Determine topography-top index on scalar-grid
1694                      k_wall = topo_top_ind(j,i,0)
1695                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1696                         kk = k - k_wall   !- lad arrays are defined flat
1697                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
1698                      ENDDO
1699                   ENDDO
1700                ENDDO
1701             ELSE
1702                DO  i = nxl, nxr
1703                   DO  j = nys, nyn
1704!--                   Determine topography-top index on scalar-grid
1705                      k_wall = topo_top_ind(j,i,0)
1706                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1707                         kk = k - k_wall   !- lad arrays are defined flat
1708                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
1709                      ENDDO
1710                   ENDDO
1711                ENDDO
1712             ENDIF
1713
1714!
1715!--       humidity
1716          CASE ( 5 )
1717             DO  i = nxl, nxr
1718                DO  j = nys, nyn
1719!
1720!--                Determine topography-top index on scalar-grid
1721                   k_wall = topo_top_ind(j,i,0)
1722
1723                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1724
1725                      kk = k - k_wall   !- lad arrays are defined flat
1726
1727                      IF ( .NOT. plant_canopy_transpiration ) THEN
1728                         ! pc_transpiration_rate is calculated in radiation model
1729                         ! in case of plant_canopy_transpiration = .T.
1730                         ! to include also the dependecy to the radiation
1731                         ! in the plant canopy box
1732                         pc_transpiration_rate(kk,j,i) =  - lsec               &
1733                                          * lad_s(kk,j,i) *                    &
1734                                          SQRT( ( 0.5_wp * ( u(k,j,i) +        &
1735                                                             u(k,j,i+1) )      &
1736                                                )**2 +                         &
1737                                                ( 0.5_wp * ( v(k,j,i) +        &
1738                                                             v(k,j+1,i) )      &
1739                                                )**2 +                         &
1740                                                ( 0.5_wp * ( w(k-1,j,i) +      &
1741                                                             w(k,j,i) )        &
1742                                                )**2                           &
1743                                              ) *                              &
1744                                          ( q(k,j,i) - lsc )
1745                      ENDIF
1746
1747                      tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
1748                   ENDDO
1749                ENDDO
1750             ENDDO
1751
1752!
1753!--       sgs-tke
1754          CASE ( 6 )
1755             DO  i = nxl, nxr
1756                DO  j = nys, nyn
1757!
1758!--                Determine topography-top index on scalar-grid
1759                   k_wall = topo_top_ind(j,i,0)
1760
1761                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1762
1763                      kk = k - k_wall   !- lad arrays are defined flat
1764                      tend(k,j,i) = tend(k,j,i) -                              &
1765                                       2.0_wp * cdc *                          &
1766                                       lad_s(kk,j,i) *                         &
1767                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1768                                                          u(k,j,i+1) )         &
1769                                             )**2 +                            &
1770                                             ( 0.5_wp * ( v(k,j,i) +           &
1771                                                          v(k,j+1,i) )         &
1772                                             )**2 +                            &
1773                                             ( 0.5_wp * ( w(k,j,i) +           &
1774                                                          w(k+1,j,i) )         &
1775                                             )**2                              &
1776                                           ) *                                 &
1777                                       e(k,j,i)
1778                   ENDDO
1779                ENDDO
1780             ENDDO
1781!
1782!--       scalar concentration
1783          CASE ( 7 )
1784             DO  i = nxl, nxr
1785                DO  j = nys, nyn
1786!
1787!--                Determine topography-top index on scalar-grid
1788                   k_wall = topo_top_ind(j,i,0)
1789
1790                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1791
1792                      kk = k - k_wall   !- lad arrays are defined flat
1793                      tend(k,j,i) = tend(k,j,i) -                              &
1794                                       lsec *                                  &
1795                                       lad_s(kk,j,i) *                         &
1796                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1797                                                          u(k,j,i+1) )         &
1798                                             )**2 +                            &
1799                                             ( 0.5_wp * ( v(k,j,i) +           &
1800                                                          v(k,j+1,i) )         &
1801                                             )**2 +                            &
1802                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
1803                                                          w(k,j,i) )           &
1804                                             )**2                              &
1805                                           ) *                                 &
1806                                       ( s(k,j,i) - lsc )
1807                   ENDDO
1808                ENDDO
1809             ENDDO   
1810
1811
1812
1813          CASE DEFAULT
1814
1815             WRITE( message_string, * ) 'wrong component: ', component
1816             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1817
1818       END SELECT
1819
1820    END SUBROUTINE pcm_tendency
1821
1822
1823!------------------------------------------------------------------------------!
1824! Description:
1825! ------------
1826!> Calculation of the tendency terms, accounting for the effect of the plant
1827!> canopy on momentum and scalar quantities.
1828!>
1829!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1830!> (defined on scalar grid), as initialized in subroutine pcm_init.
1831!> The lad on the w-grid is vertically interpolated from the surrounding
1832!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1833!> k = pch_index. Here, the lad is zero.
1834!>
1835!> The canopy drag must be limited (previously accounted for by calculation of
1836!> a limiting canopy timestep for the determination of the maximum LES timestep
1837!> in subroutine timestep), since it is physically impossible that the canopy
1838!> drag alone can locally change the sign of a velocity component. This
1839!> limitation is realized by calculating preliminary tendencies and velocities.
1840!> It is subsequently checked if the preliminary new velocity has a different
1841!> sign than the current velocity. If so, the tendency is limited in a way that
1842!> the velocity can at maximum be reduced to zero by the canopy drag.
1843!>
1844!>
1845!> Call for grid point i,j
1846!------------------------------------------------------------------------------!
1847    SUBROUTINE pcm_tendency_ij( i, j, component )
1848
1849
1850       USE control_parameters,                                                 &
1851           ONLY:  dt_3d, message_string
1852
1853       USE kinds
1854
1855       IMPLICIT NONE
1856
1857       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1858       INTEGER(iwp) ::  i         !< running index
1859       INTEGER(iwp) ::  j         !< running index
1860       INTEGER(iwp) ::  k         !< running index
1861       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
1862       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1863
1864       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1865       REAL(wp) ::  lad_local !< local lad value
1866       REAL(wp) ::  pre_tend  !< preliminary tendency
1867       REAL(wp) ::  pre_u     !< preliminary u-value
1868       REAL(wp) ::  pre_v     !< preliminary v-value
1869       REAL(wp) ::  pre_w     !< preliminary w-value
1870
1871
1872       ddt_3d = 1.0_wp / dt_3d
1873!
1874!--    Compute drag for the three velocity components and the SGS-TKE
1875       SELECT CASE ( component )
1876
1877!
1878!--       u-component
1879          CASE ( 1 )
1880!
1881!--          Determine topography-top index on u-grid
1882             k_wall = topo_top_ind(j,i,1)
1883             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
1884
1885                kk = k - k_wall  !- lad arrays are defined flat
1886
1887!
1888!--             In order to create sharp boundaries of the plant canopy,
1889!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1890!--             rather than being interpolated from the surrounding lad_s,
1891!--             because this would yield smaller lad at the canopy boundaries
1892!--             than inside of the canopy.
1893!--             For the same reason, the lad at the rightmost(i+1)canopy
1894!--             boundary on the u-grid equals lad_s(k,j,i).
1895                lad_local = lad_s(kk,j,i)
1896                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
1897                   lad_local = lad_s(kk,j,i-1)
1898                ENDIF
1899
1900                pre_tend = 0.0_wp
1901                pre_u = 0.0_wp
1902!
1903!--             Calculate preliminary value (pre_tend) of the tendency
1904                pre_tend = - cdc *                                             &
1905                             lad_local *                                       &   
1906                             SQRT( u(k,j,i)**2 +                               &
1907                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
1908                                                 v(k,j,i)    +                 &
1909                                                 v(k,j+1,i)  +                 &
1910                                                 v(k,j+1,i-1) )                &
1911                                   )**2 +                                      &
1912                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
1913                                                 w(k-1,j,i)   +                &
1914                                                 w(k,j,i-1)   +                &
1915                                                 w(k,j,i) )                    &
1916                                   )**2                                        &
1917                                 ) *                                           &
1918                             u(k,j,i)
1919
1920!
1921!--             Calculate preliminary new velocity, based on pre_tend
1922                pre_u = u(k,j,i) + dt_3d * pre_tend
1923!
1924!--             Compare sign of old velocity and new preliminary velocity,
1925!--             and in case the signs are different, limit the tendency
1926                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1927                   pre_tend = - u(k,j,i) * ddt_3d
1928                ELSE
1929                   pre_tend = pre_tend
1930                ENDIF
1931!
1932!--             Calculate final tendency
1933                tend(k,j,i) = tend(k,j,i) + pre_tend
1934             ENDDO
1935
1936
1937!
1938!--       v-component
1939          CASE ( 2 )
1940!
1941!--          Determine topography-top index on v-grid
1942             k_wall = topo_top_ind(j,i,2)
1943
1944             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
1945
1946                kk = k - k_wall  !- lad arrays are defined flat
1947!
1948!--             In order to create sharp boundaries of the plant canopy,
1949!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1950!--             rather than being interpolated from the surrounding lad_s,
1951!--             because this would yield smaller lad at the canopy boundaries
1952!--             than inside of the canopy.
1953!--             For the same reason, the lad at the northmost(j+1)canopy
1954!--             boundary on the v-grid equals lad_s(k,j,i).
1955                lad_local = lad_s(kk,j,i)
1956                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
1957                   lad_local = lad_s(kk,j-1,i)
1958                ENDIF
1959
1960                pre_tend = 0.0_wp
1961                pre_v = 0.0_wp
1962!
1963!--             Calculate preliminary value (pre_tend) of the tendency
1964                pre_tend = - cdc *                                             &
1965                             lad_local *                                       &
1966                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
1967                                                 u(k,j-1,i+1) +                &
1968                                                 u(k,j,i)     +                &
1969                                                 u(k,j,i+1) )                  &
1970                                   )**2 +                                      &
1971                                   v(k,j,i)**2 +                               &
1972                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
1973                                                 w(k-1,j,i)   +                &
1974                                                 w(k,j-1,i)   +                &
1975                                                 w(k,j,i) )                    &
1976                                   )**2                                        &
1977                                 ) *                                           &
1978                             v(k,j,i)
1979
1980!
1981!--             Calculate preliminary new velocity, based on pre_tend
1982                pre_v = v(k,j,i) + dt_3d * pre_tend
1983!
1984!--             Compare sign of old velocity and new preliminary velocity,
1985!--             and in case the signs are different, limit the tendency
1986                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1987                   pre_tend = - v(k,j,i) * ddt_3d
1988                ELSE
1989                   pre_tend = pre_tend
1990                ENDIF
1991!
1992!--             Calculate final tendency
1993                tend(k,j,i) = tend(k,j,i) + pre_tend
1994             ENDDO
1995
1996
1997!
1998!--       w-component
1999          CASE ( 3 )
2000!
2001!--          Determine topography-top index on w-grid
2002             k_wall = topo_top_ind(j,i,3)
2003
2004             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
2005
2006                kk = k - k_wall  !- lad arrays are defined flat
2007
2008                pre_tend = 0.0_wp
2009                pre_w = 0.0_wp
2010!
2011!--             Calculate preliminary value (pre_tend) of the tendency
2012                pre_tend = - cdc *                                             &
2013                             (0.5_wp *                                         &
2014                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
2015                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
2016                                                 u(k,j,i+1)  +                 &
2017                                                 u(k+1,j,i)  +                 &
2018                                                 u(k+1,j,i+1) )                &
2019                                   )**2 +                                      &
2020                                   ( 0.25_wp * ( v(k,j,i)    +                 &
2021                                                 v(k,j+1,i)  +                 &
2022                                                 v(k+1,j,i)  +                 &
2023                                                 v(k+1,j+1,i) )                &
2024                                   )**2 +                                      &
2025                                   w(k,j,i)**2                                 &
2026                                 ) *                                           &
2027                             w(k,j,i)
2028!
2029!--             Calculate preliminary new velocity, based on pre_tend
2030                pre_w = w(k,j,i) + dt_3d * pre_tend
2031!
2032!--             Compare sign of old velocity and new preliminary velocity,
2033!--             and in case the signs are different, limit the tendency
2034                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
2035                   pre_tend = - w(k,j,i) * ddt_3d
2036                ELSE
2037                   pre_tend = pre_tend
2038                ENDIF
2039!
2040!--             Calculate final tendency
2041                tend(k,j,i) = tend(k,j,i) + pre_tend
2042             ENDDO
2043
2044!
2045!--       potential temperature
2046          CASE ( 4 )
2047!
2048!--          Determine topography-top index on scalar grid
2049             k_wall = topo_top_ind(j,i,0)
2050
2051             IF ( humidity ) THEN
2052                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2053                   kk = k - k_wall  !- lad arrays are defined flat
2054                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) -       &
2055                                 pc_latent_rate(kk,j,i)
2056                ENDDO
2057             ELSE
2058                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2059                   kk = k - k_wall  !- lad arrays are defined flat
2060                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
2061                ENDDO
2062             ENDIF
2063
2064!
2065!--       humidity
2066          CASE ( 5 )
2067!
2068!--          Determine topography-top index on scalar grid
2069             k_wall = topo_top_ind(j,i,0)
2070
2071             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2072                kk = k - k_wall  !- lad arrays are defined flat
2073                IF ( .NOT. plant_canopy_transpiration ) THEN
2074                   ! pc_transpiration_rate is calculated in radiation model
2075                   ! in case of plant_canopy_transpiration = .T.
2076                   ! to include also the dependecy to the radiation
2077                   ! in the plant canopy box
2078                   pc_transpiration_rate(kk,j,i) = - lsec                      &
2079                                    * lad_s(kk,j,i) *                          &
2080                                    SQRT( ( 0.5_wp * ( u(k,j,i) +              &
2081                                                       u(k,j,i+1) )            &
2082                                          )**2  +                              &
2083                                          ( 0.5_wp * ( v(k,j,i) +              &
2084                                                       v(k,j+1,i) )            &
2085                                          )**2 +                               &
2086                                          ( 0.5_wp * ( w(k-1,j,i) +            &
2087                                                       w(k,j,i) )              &
2088                                          )**2                                 &
2089                                        ) *                                    &
2090                                    ( q(k,j,i) - lsc )
2091                ENDIF
2092
2093                tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
2094
2095             ENDDO   
2096
2097!
2098!--       sgs-tke
2099          CASE ( 6 )
2100!
2101!--          Determine topography-top index on scalar grid
2102             k_wall = topo_top_ind(j,i,0)
2103
2104             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2105
2106                kk = k - k_wall
2107                tend(k,j,i) = tend(k,j,i) -                                    &
2108                                 2.0_wp * cdc *                                &
2109                                 lad_s(kk,j,i) *                               &
2110                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2111                                                    u(k,j,i+1) )               &
2112                                       )**2 +                                  & 
2113                                       ( 0.5_wp * ( v(k,j,i) +                 &
2114                                                    v(k,j+1,i) )               &
2115                                       )**2 +                                  &
2116                                       ( 0.5_wp * ( w(k,j,i) +                 &
2117                                                    w(k+1,j,i) )               &
2118                                       )**2                                    &
2119                                     ) *                                       &
2120                                 e(k,j,i)
2121             ENDDO
2122             
2123!
2124!--       scalar concentration
2125          CASE ( 7 )
2126!
2127!--          Determine topography-top index on scalar grid
2128             k_wall = topo_top_ind(j,i,0)
2129
2130             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2131
2132                kk = k - k_wall
2133                tend(k,j,i) = tend(k,j,i) -                                    &
2134                                 lsec *                                        &
2135                                 lad_s(kk,j,i) *                               &
2136                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2137                                                    u(k,j,i+1) )               &
2138                                       )**2  +                                 &
2139                                       ( 0.5_wp * ( v(k,j,i) +                 &
2140                                                    v(k,j+1,i) )               &
2141                                       )**2 +                                  &
2142                                       ( 0.5_wp * ( w(k-1,j,i) +               &
2143                                                    w(k,j,i) )                 &
2144                                       )**2                                    &
2145                                     ) *                                       &
2146                                 ( s(k,j,i) - lsc )
2147             ENDDO               
2148
2149       CASE DEFAULT
2150
2151          WRITE( message_string, * ) 'wrong component: ', component
2152          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
2153
2154       END SELECT
2155
2156    END SUBROUTINE pcm_tendency_ij
2157
2158
2159
2160 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.