Home

Context Navigation

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

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