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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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