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

Last change on this file since 4281 was 4279, checked in by scharf, 5 years ago

unused variables removed

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