Home

Context Navigation

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

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