Home

Context Navigation

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

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