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

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

Coupling of indoor model to atmosphere; output of indoor temperatures and waste heat; enable restarts with indoor model; bugfix plant transpiration; bugfix - missing calculation of 10cm temperature

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