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

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

radiation: Take external radiation input from root domain dynamic file if no dynamic input file is provided for each nested domain; radiation: Combine MPI_ALLREDUCE calls to reduce latency; land_surface + plant_canopy: Give specific error numbers; land-surface: Adjust error messages for local checks; init_3d_model: Deallocate temporary string array for netcdf-data input since it may be re-used in other modules for different input files

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