Home

Context Navigation

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

Last change on this file since 4172 was 4168, checked in by suehring, 5 years ago
Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources
Property svn:keywords set to `Id`
File size: 94.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |