Home

Context Navigation

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

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