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

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

Bugfixes in 3d data output of plant canopy variables

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