Home

Context Navigation

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

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