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

Last change on this file since 4213 was 4205, checked in by suehring, 5 years ago

plant canopy: Missing working precision + bugfix in calculation of wind speed; surface data output: Correct x,y-coordinates of vertical surfaces in netcdf output; Change definition of azimuth angle, reference is north 0 degree

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