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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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