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

Last change on this file since 4258 was 4258, checked in by suehring, 21 months ago

Land-surface model: Revise limitation for soil moisture in case it exceeds its saturation value; Revise initialization of soil moisture and temperature in a nested run in case dynamic input information is available. This case, the soil within the child domains can be initialized separately; As part of this revision, migrate the netcdf input of soil temperature / moisture to this module, as well as the routine to inter/extrapolate soil profiles between different grids.; Plant-canopy: Check if any LAD is prescribed when plant-canopy model is applied, in order to avoid crashes in the radiation transfer model; Surface-layer fluxes: Initialization of Obukhov length also at vertical surfaces (if allocated); Urban-surface model: Add checks to ensure that relative fractions of walls, windowns and green surfaces sum-u to one; Revise message calls dealing with local checks

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