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

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