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

Last change on this file since 4278 was 4278, checked in by scharf, 20 months ago

changed check for static driver and fixed bugs in initialization and header of plant canopy model

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