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

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

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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