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

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

Minor adjustments in error messages and error numbers. Some typos are corrected.

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