Home

Context Navigation

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

Last change on this file since 4842 was 4842, checked in by raasch, 3 years ago
reading of namelist file and actions in case of namelist errors revised so that statement labels and goto statements are not required any more, deprecated namelists removed
Property svn:keywords set to `Id`
File size: 127.5 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 terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2021 Leibniz Universitaet Hannover
17	! Copyright 2017-2021 Institute of Computer Science of the
18	! Czech Academy of Sciences, Prague
19	!--------------------------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: plant_canopy_model_mod.f90 4842 2021-01-14 10:42:28Z raasch $
28	! reading of namelist file and actions in case of namelist errors revised so that statement labels
29	! and goto statements are not required any more,
30	! deprecaed namelist name removed
31	!
32	! 4828 2021-01-05 11:21:41Z Giersch
33	! file re-formatted to follow the PALM coding standard
34	!
35	! 4770 2020-11-03 14:04:53Z suehring
36	! Consider basal-area density as an additional sink for momentum in the prognostic equations for
37	! momentum and SGS TKE
38	!
39	! 4768 2020-11-02 19:11:23Z suehring
40	! Enable 3D data output also with 64-bit precision
41	!
42	! 4671 2020-09-09 20:27:58Z pavelkrc
43	! Implementation of downward facing USM and LSM surfaces
44	!
45	! 4535 2020-05-15 12:07:23Z raasch
46	! bugfix for restart data format query
47	!
48	! 4525 2020-05-10 17:05:07Z raasch
49	! bugfix for reading/writing pcm_...rate_av with MPI-IO
50	!
51	! 4517 2020-05-03 14:29:30Z raasch
52	! added restart with MPI-IO for reading local arrays
53	!
54	! 4515 2020-04-30 16:37:18Z suehring
55	! Rename error number again since this was recently given in -r 4511
56	!
57	! 4514 2020-04-30 16:29:59Z suehring
58	! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this, pch_index is now
59	! output for a restart run. Therefore, define global restart routines.
60	! - Error message number renamed and check for PA0505 revised in order to also consider natural
61	! surfaces with plant-canopy.
62	!
63	! 4495 2020-04-13 20:11:20Z raasch
64	! restart data handling with MPI-IO added
65	!
66	! 4457 2020-03-11 14:20:43Z raasch
67	!
68	! use statement for exchange horiz added
69	! (salim) removed the error message PA0672 to consider PC 3d data via ascii file
70	!
71	! 4392 2020-01-31 16:14:57Z pavelkrc (resler)
72	! Make pcm_heatrate_av, pcm_latentrate_av public to allow calculation of averaged Bowen ratio in the
73	! user procedure
74	!
75	! 4381 2020-01-20 13:51:46Z suehring
76	! Give error message 313 only once
77	!
78	! 4363 2020-01-07 18:11:28Z suehring
79	! Fix for last commit
80	!
81	! 4362 2020-01-07 17:15:02Z suehring
82	! Input of plant canopy variables from static driver moved to plant-canopy model
83	!
84	! 4361 2020-01-07 12:22:38Z suehring
85	! - Remove unused arrays in pmc_rrd_local
86	! - Remove one exchange of ghost points
87	!
88	! 4360 2020-01-07 11:25:50Z suehring
89	! - Bugfix, read restart data for time-averaged pcm output quantities
90	! - Output of plant-canopy quantities will fill values
91	!
92	! 4356 2019-12-20 17:09:33Z suehring
93	! Correct single message call, local check must be given by the respective mpi rank.
94	!
95	! 4346 2019-12-18 11:55:56Z motisi
96	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
97	! information used in wall_flags_static_0
98	!
99	! 4342 2019-12-16 13:49:14Z Giersch
100	! Use statements moved to module level, ocean dependency removed, redundant variables removed
101	!
102	! 4341 2019-12-16 10:43:49Z motisi
103	! - Unification of variable names: pc_-variables now pcm_-variables (pc_latent_rate,
104	! pc_heating_rate, pc_transpiration_rate)
105	! - Removal of pcm_bowenratio output
106	! - Renamed canopy-mode 'block' to 'homogeneous'
107	! - Renamed value 'read_from_file_3d' to 'read_from_file'
108	! - Removal of confusing comment lines
109	! - Replacement of k_wall by topo_top_ind
110	! - Removal of Else-Statement in tendency-calculation
111	!
112	! 4335 2019-12-12 16:39:05Z suehring
113	! Fix for LAD at building edges also implemented in vector branch.
114	!
115	! 4331 2019-12-10 18:25:02Z suehring
116	! Typo corrected
117	!
118	! 4329 2019-12-10 15:46:36Z motisi
119	! Renamed wall_flags_0 to wall_flags_static_0
120	!
121	! 4314 2019-11-29 10:29:20Z suehring
122	! - Bugfix, plant canopy was still considered at building edges on for the u- and v-component.
123	! - Relax restriction of LAD on building tops. LAD is only omitted at locations where building grid
124	! points emerged artificially by the topography filtering.
125	!
126	! 4309 2019-11-26 18:49:59Z suehring
127	! Typo
128	!
129	! 4302 2019-11-22 13:15:56Z suehring
130	! Omit tall canopy mapped on top of buildings
131	!
132	! 4279 2019-10-29 08:48:17Z scharf
133	! unused variables removed
134	!
135	! 4258 2019-10-07 13:29:08Z scharf
136	! changed check for static driver and fixed bugs in initialization and header
137	!
138	! 4258 2019-10-07 13:29:08Z suehring
139	! Check if any LAD is prescribed when plant-canopy model is applied.
140	!
141	! 4226 2019-09-10 17:03:24Z suehring
142	! Bugfix, missing initialization of heating rate
143	!
144	! 4221 2019-09-09 08:50:35Z suehring
145	! Further bugfix in 3d data output for plant canopy
146	!
147	! 4216 2019-09-04 09:09:03Z suehring
148	! Bugfixes in 3d data output
149	!
150	! 4205 2019-08-30 13:25:00Z suehring
151	! Missing working precision + bugfix in calculation of wind speed
152	!
153	! 4188 2019-08-26 14:15:47Z suehring
154	! Minor adjustment in error number
155	!
156	! 4187 2019-08-26 12:43:15Z suehring
157	! Give specific error numbers instead of PA0999
158	!
159	! 4182 2019-08-22 15:20:23Z scharf
160	! Corrected "Former revisions" section
161	!
162	! 4168 2019-08-16 13:50:17Z suehring
163	! Replace function get_topography_top_index by topo_top_ind
164	!
165	! 4127 2019-07-30 14:47:10Z suehring
166	! Output of 3D plant canopy variables changed. It is now relative to the local terrain rather than
167	! located at the acutal vertical level in the model. This way, the vertical dimension of the output
168	! can be significantly reduced.
169	! (merge from branch resler)
170	!
171	! 3885 2019-04-11 11:29:34Z kanani
172	! Changes related to global restructuring of location messages and introduction of additional debug
173	! messages
174	!
175	! 3864 2019-04-05 09:01:56Z monakurppa
176	! unsed variables removed
177	!
178	! 3745 2019-02-15 18:57:56Z suehring
179	! Bugfix in transpiration, floating invalid when temperature becomes > 40 degrees
180	!
181	! 3744 2019-02-15 18:38:58Z suehring
182	! Some interface calls moved to module_interface + cleanup
183	!
184	! 3655 2019-01-07 16:51:22Z knoop
185	! unused variables removed
186	!
187	! 138 2007-11-28 10:03:58Z letzel
188	! Initial revision
189	!
190	! Description:
191	! ------------
192	!> 1) Initialization of the canopy model, e.g. construction of leaf area density profile
193	!> (subroutine pcm_init).
194	!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration due to canopy
195	!> elements (subroutine pcm_tendency).
196	!
197	! @todo - precalculate constant terms in pcm_calc_transpiration_rate
198	! @todo - unify variable names (pcm_, pc_, ...)
199	! @todo - get rid-off dependency on radiation model
200	!--------------------------------------------------------------------------------------------------!
201	MODULE plant_canopy_model_mod
202
203	USE arrays_3d, &
204	ONLY: dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
205
206	USE basic_constants_and_equations_mod, &
207	ONLY: c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
208
209	USE bulk_cloud_model_mod, &
210	ONLY: bulk_cloud_model, microphysics_seifert
211
212	USE control_parameters, &
213	ONLY: average_count_3d, &
214	coupling_char, &
215	debug_output, &
216	dt_3d, &
217	dz, &
218	humidity, &
219	land_surface, &
220	length, &
221	message_string, &
222	ocean_mode, &
223	passive_scalar, &
224	plant_canopy, &
225	restart_data_format_output, &
226	restart_string, &
227	urban_surface
228
229	USE grid_variables, &
230	ONLY: dx, dy
231
232	USE indices, &
233	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, nz, nzb, nzt, &
234	topo_top_ind, wall_flags_total_0
235
236	USE kinds
237
238	USE netcdf_data_input_mod, &
239	ONLY: char_fill, &
240	check_existence, &
241	close_input_file, &
242	get_attribute, &
243	get_dimension_length, &
244	get_variable, &
245	input_file_static, &
246	input_pids_static, &
247	inquire_num_variables, &
248	inquire_variable_names, &
249	num_var_pids, &
250	open_read_file, &
251	pids_id, &
252	real_3d, &
253	vars_pids
254
255	USE pegrid
256
257	USE restart_data_mpi_io_mod, &
258	ONLY: rd_mpi_io_check_array, &
259	rrd_mpi_io, &
260	wrd_mpi_io
261
262	USE surface_mod, &
263	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
264
265
266	IMPLICIT NONE
267
268	CHARACTER (LEN=30) :: canopy_mode = 'homogeneous' !< canopy coverage
269	INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index
270
271	INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
272
273	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top
274
275	LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
276	LOGICAL :: plant_canopy_transpiration = .FALSE. !< flag to switch calculation of transpiration and corresponding latent heat
277	!< for resolved plant canopy inside radiation model
278	!< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
279
280	REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation
281	REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation
282	REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
283	REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux
284	REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag
285	REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient
286	REAL(wp) :: lad_surface = 0.0_wp !< lad surface value
287	REAL(wp) :: lad_type_coef(0:10) = 1.0_wp !< multiplicative coeficients for particular types
288	!< of plant canopy (e.g. deciduous tree during winter)
289	REAL(wp) :: lad_vertical_gradient(10) = 0.0_wp !< lad gradient
290	REAL(wp) :: lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
291	REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc.
292	REAL(wp) :: leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff.
293	REAL(wp) :: leaf_surface_conc = 0.0_wp !< leaf surface concentration
294
295	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density
296	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad
297
298	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: bad_s !< basal-area density on scalar-grid
299	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.
300	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid
301	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_heating_rate !< plant canopy heating rate
302	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_heatrate_av !< array for averaging plant canopy sensible heating rate
303	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_latent_rate !< plant canopy latent heating rate
304	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_latentrate_av !< array for averaging plant canopy latent heating rate
305	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_transpiration_rate !< plant canopy transpiration rate
306	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
307
308	TYPE(real_3d) :: basal_area_density_f !< input variable for basal area density - resolved vegetation
309	TYPE(real_3d) :: leaf_area_density_f !< input variable for leaf area density - resolved vegetation
310	TYPE(real_3d) :: root_area_density_lad_f !< input variable for root area density - resolved vegetation
311
312	SAVE
313
314	PRIVATE
315
316	!
317	!-- Public functions
318	PUBLIC pcm_calc_transpiration_rate, &
319	pcm_check_data_output, &
320	pcm_check_parameters, &
321	pcm_3d_data_averaging, &
322	pcm_data_output_3d, &
323	pcm_define_netcdf_grid, &
324	pcm_header, &
325	pcm_init, &
326	pcm_parin, &
327	pcm_rrd_global, &
328	pcm_rrd_local, &
329	pcm_tendency, &
330	pcm_wrd_global, &
331	pcm_wrd_local
332
333	!
334	!-- Public variables and constants
335	PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate, pcm_latent_rate, &
336	canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index, plant_canopy_transpiration, &
337	pcm_heatrate_av, pcm_latentrate_av
338
339	INTERFACE pcm_calc_transpiration_rate
340	MODULE PROCEDURE pcm_calc_transpiration_rate
341	END INTERFACE pcm_calc_transpiration_rate
342
343	INTERFACE pcm_check_data_output
344	MODULE PROCEDURE pcm_check_data_output
345	END INTERFACE pcm_check_data_output
346
347	INTERFACE pcm_check_parameters
348	MODULE PROCEDURE pcm_check_parameters
349	END INTERFACE pcm_check_parameters
350
351	INTERFACE pcm_3d_data_averaging
352	MODULE PROCEDURE pcm_3d_data_averaging
353	END INTERFACE pcm_3d_data_averaging
354
355	INTERFACE pcm_data_output_3d
356	MODULE PROCEDURE pcm_data_output_3d
357	END INTERFACE pcm_data_output_3d
358
359	INTERFACE pcm_define_netcdf_grid
360	MODULE PROCEDURE pcm_define_netcdf_grid
361	END INTERFACE pcm_define_netcdf_grid
362
363	INTERFACE pcm_header
364	MODULE PROCEDURE pcm_header
365	END INTERFACE pcm_header
366
367	INTERFACE pcm_init
368	MODULE PROCEDURE pcm_init
369	END INTERFACE pcm_init
370
371	INTERFACE pcm_parin
372	MODULE PROCEDURE pcm_parin
373	END INTERFACE pcm_parin
374
375	INTERFACE pcm_read_plant_canopy_3d
376	MODULE PROCEDURE pcm_read_plant_canopy_3d
377	END INTERFACE pcm_read_plant_canopy_3d
378
379	INTERFACE pcm_rrd_local
380	MODULE PROCEDURE pcm_rrd_local_ftn
381	MODULE PROCEDURE pcm_rrd_local_mpi
382	END INTERFACE pcm_rrd_local
383
384	INTERFACE pcm_rrd_global
385	MODULE PROCEDURE pcm_rrd_global_ftn
386	MODULE PROCEDURE pcm_rrd_global_mpi
387	END INTERFACE pcm_rrd_global
388
389	INTERFACE pcm_tendency
390	MODULE PROCEDURE pcm_tendency
391	MODULE PROCEDURE pcm_tendency_ij
392	END INTERFACE pcm_tendency
393
394	INTERFACE pcm_wrd_local
395	MODULE PROCEDURE pcm_wrd_local
396	END INTERFACE pcm_wrd_local
397
398	INTERFACE pcm_wrd_global
399	MODULE PROCEDURE pcm_wrd_global
400	END INTERFACE pcm_wrd_global
401
402
403	CONTAINS
404
405
406	!--------------------------------------------------------------------------------------------------!
407	! Description:
408	! ------------
409	!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart with
410	!> parametrizations described in Daudet et al. (1999; Agricult. and Forest Meteorol. 97) and Ngao,
411	!> Adam and Saudreau (2017; Agricult. and Forest Meteorol 237-238). Model functions f1-f4 were
412	!> adapted from Stewart (1998; Agric. and Forest. Meteorol. 43) instead, because they are valid for
413	!> broader intervals of values. Funcion f4 used in form present in van Wijk et al. (1998; Tree
414	!> Physiology 20).
415	!>
416	!> This subroutine is called from subroutine radiation_interaction after the calculation of
417	!> radiation in plant canopy boxes.
418	!> (arrays pcbinsw and pcbinlw).
419	!>
420	!--------------------------------------------------------------------------------------------------!
421	SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
422
423	!
424	!-- Input parameters
425	INTEGER(iwp), INTENT(IN) :: i, j, k, kk !< indices of the pc gridbox
426	REAL(wp), INTENT(IN) :: pcblw !< lw radiation in gridbox (W)
427	REAL(wp), INTENT(IN) :: pcbsw !< sw radiation in gridbox (W)
428	REAL(wp), INTENT(OUT) :: pcblh !< latent heat from transpiration dT/dt (K/s)
429	REAL(wp), INTENT(OUT) :: pcbtr !< transpiration rate dq/dt (kg/kg/s)
430
431	!-- Variables and parameters for calculation of transpiration rate
432	REAL(wp), PARAMETER :: gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
433	REAL(wp), PARAMETER :: g_s_max = 0.01 !< maximum stomatal conductivity (m/s)
434	REAL(wp), PARAMETER :: m_soil = 0.4_wp !< soil water content (needs to adjust or take from LSM)
435	REAL(wp), PARAMETER :: m_wilt = 0.01_wp !< wilting point soil water content (needs to adjust or take from LSM)
436	REAL(wp), PARAMETER :: m_sat = 0.51_wp !< saturation soil water content (needs to adjust or take from LSM)
437	REAL(wp), PARAMETER :: t2_min = 0.0_wp !< minimal temperature for calculation of f2
438	REAL(wp), PARAMETER :: t2_max = 40.0_wp !< maximal temperature for calculation of f2
439
440	REAL(wp) :: d_fact
441	REAL(wp) :: e_eq
442	REAL(wp) :: e_imp
443	REAL(wp) :: evapor_rate
444	REAL(wp) :: f1
445	REAL(wp) :: f2
446	REAL(wp) :: f3
447	REAL(wp) :: f4
448	REAL(wp) :: g_b
449	REAL(wp) :: g_s
450	REAL(wp) :: rad
451	REAL(wp) :: rswc
452	REAL(wp) :: sat_press
453	REAL(wp) :: sat_press_d
454	REAL(wp) :: temp
455	REAL(wp) :: v_lad
456	REAL(wp) :: vpd
457	REAL(wp) :: wind_speed
458
459	!
460	!-- Temperature (deg C)
461	temp = pt(k,j,i) * exner(k) - degc_to_k
462	!
463	!-- Coefficient for conversion of radiation to grid to radiation to unit leaves surface
464	v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
465	!
466	!-- Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
467	!-- There are updated formulas available, kept consistent with the rest of the parametrization
468	sat_press = 610.8_wp * EXP( 17.27_wp * temp / ( temp + 237.3_wp ) )
469	!
470	!-- Saturation pressure derivative (derivative of the above)
471	sat_press_d = sat_press * 17.27_wp * 237.3_wp / ( temp + 237.3_wp )**2
472	!
473	!-- Wind speed
474	wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
475	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 + &
476	( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
477	!
478	!-- Aerodynamic conductivity (Daudet et al. (1999) eq. 14
479	g_b = 0.01_wp * wind_speed + 0.0071_wp
480	!
481	!-- Radiation flux per leaf surface unit
482	rad = pcbsw * v_lad
483	!
484	!-- First function for calculation of stomatal conductivity (radiation dependency)
485	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
486	f1 = rad * ( 1000.0_wp + 42.1_wp ) / 1000.0_wp / ( rad + 42.1_wp )
487	!
488	!-- Second function for calculation of stomatal conductivity (temperature dependency)
489	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
490	f2 = MAX( t2_min, ( temp - t2_min ) * MAX( 0.0_wp, t2_max - temp )**( ( t2_max - 16.9_wp ) / &
491	( 16.9_wp - t2_min ) ) &
492	/ ( ( 16.9_wp - t2_min ) * ( t2_max - 16.9_wp )**( ( t2_max - 16.9_wp ) / &
493	( 16.9_wp - t2_min ) ) ) )
494	!
495	!-- Water pressure deficit
496	!-- Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
497	vpd = MAX( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
498	!
499	!-- Third function for calculation of stomatal conductivity (water pressure deficit dependency)
500	!-- Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
501	!-- The coefficients of the linear dependence should better correspond to broad-leaved trees than
502	!-- the coefficients from Stewart (1988) which correspond to conifer trees.
503	vpd = MIN( MAX( vpd, 770.0_wp ), 3820.0_wp )
504	f3 = -2E-4_wp * vpd + 1.154_wp
505	!
506	!-- Fourth function for calculation of stomatal conductivity (soil moisture dependency)
507	!-- Residual soil water content
508	!-- van Wijk et al. (1998; Tree Physiology 20) eq. 7
509	!-- TODO - over LSM surface might be calculated from LSM parameters
510	rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
511	!
512	!-- van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of
513	!-- Stewart(1988))
514	f4 = MAX( 0.0_wp, MIN( 1.0_wp - 0.041_wp * EXP( 3.2_wp * rswc ), 1.0_wp - 0.041_wp ) )
515	!
516	!-- Stomatal conductivity
517	!-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
518	!-- (notation according to Ngao, Adam and Saudreau (2017) and others)
519	g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
520	!
521	!-- Decoupling factor
522	!-- Daudet et al. (1999) eq. 6
523	d_fact = ( sat_press_d / gama_psychr + 2.0_wp ) / &
524	( sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
525	!
526	!-- Equilibrium evaporation rate
527	!-- Daudet et al. (1999) eq. 4
528	e_eq = ( pcbsw + pcblw ) * v_lad * sat_press_d / &
529	gama_psychr / ( sat_press_d / gama_psychr + 2.0_wp ) / l_v
530	!
531	!-- Imposed evaporation rate
532	!-- Daudet et al. (1999) eq. 5
533	e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
534	!
535	!-- Evaporation rate
536	!-- Daudet et al. (1999) eq. 3
537	!-- (evaporation rate is limited to non-negative values)
538	evapor_rate = MAX( d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp )
539	!
540	!-- Conversion of evaporation rate to q tendency in gridbox
541	!-- dq/dt = E * LAD * V_g / (rho_air * V_g)
542	pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k) !-- = dq/dt
543	!
544	!-- latent heat from evaporation
545	pcblh = pcbtr * lv_d_cp !-- = - dT/dt
546
547	END SUBROUTINE pcm_calc_transpiration_rate
548
549
550	!--------------------------------------------------------------------------------------------------!
551	! Description:
552	! ------------
553	!> Check data output for plant canopy model
554	!--------------------------------------------------------------------------------------------------!
555	SUBROUTINE pcm_check_data_output( var, unit )
556
557	CHARACTER (LEN=*) :: unit !<
558	CHARACTER (LEN=*) :: var !<
559
560
561	SELECT CASE ( TRIM( var ) )
562
563	CASE ( 'pcm_heatrate' )
564	!
565	!-- Output of heatrate can be only done if it is explicitely set by cthf, or parametrized by
566	!-- absorption of radiation. The latter, however, is only available if radiation_interactions
567	!-- are on. Note, these are enabled if land-surface or urban-surface is switched-on. Using
568	!-- radiation_interactions_on directly is not possible since it belongs to the
569	!-- radition_model, which in turn depends on the plant-canopy model, creating circular
570	!-- dependencies.
571	IF ( cthf == 0.0_wp .AND. ( .NOT. urban_surface .AND. .NOT. land_surface ) ) THEN
572	message_string = 'output of "' // TRIM( var ) // '" requi' // &
573	'res setting of parameter cthf /= 0.0'
574	CALL message( 'pcm_check_data_output', 'PA0718', 1, 2, 0, 6, 0 )
575	ENDIF
576	unit = 'K s-1'
577
578	CASE ( 'pcm_transpirationrate' )
579	unit = 'kg kg-1 s-1'
580
581	CASE ( 'pcm_latentrate' )
582	unit = 'K s-1'
583
584	CASE ( 'pcm_bad', 'pcm_lad' )
585	unit = 'm2 m-3'
586
587
588	CASE DEFAULT
589	unit = 'illegal'
590
591	END SELECT
592
593
594	END SUBROUTINE pcm_check_data_output
595
596
597	!--------------------------------------------------------------------------------------------------!
598	! Description:
599	! ------------
600	!> Check parameters routine for plant canopy model
601	!--------------------------------------------------------------------------------------------------!
602	SUBROUTINE pcm_check_parameters
603
604	IF ( ocean_mode ) THEN
605	message_string = 'plant_canopy = .TRUE. is not allowed in the ocean'
606	CALL message( 'pcm_check_parameters', 'PA0696', 1, 2, 0, 6, 0 )
607	ENDIF
608
609	IF ( canopy_drag_coeff == 0.0_wp ) THEN
610	message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
611	'coefficient & given value is canopy_drag_coeff = 0.0'
612	CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
613	ENDIF
614
615	IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp ) .OR. &
616	beta_lad /= 9999999.9_wp .AND. alpha_lad == 9999999.9_wp ) THEN
617	message_string = 'using the beta function for the construction ' // &
618	'of the leaf area density profile requires ' // &
619	'both alpha_lad and beta_lad to be /= 9999999.9'
620	CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
621	ENDIF
622
623	IF ( calc_beta_lad_profile .AND. lai_beta == 0.0_wp ) THEN
624	message_string = 'using the beta function for the construction ' // &
625	'of the leaf area density profile requires ' // &
626	'a non-zero lai_beta, but given value is ' // &
627	'lai_beta = 0.0'
628	CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
629	ENDIF
630
631	IF ( calc_beta_lad_profile .AND. lad_surface /= 0.0_wp ) THEN
632	message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
633	'combined with beta_lad /= 9999999.9 ' // &
634	'and lad_surface /= 0.0 is not possible, ' // &
635	'use either vertical gradients or the beta ' // &
636	'function for the construction of the leaf area '// &
637	'density profile'
638	CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
639	ENDIF
640
641	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
642	message_string = 'plant_canopy = .TRUE. requires cloud_scheme /= seifert_beheng'
643	CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
644	ENDIF
645
646	END SUBROUTINE pcm_check_parameters
647
648
649	!--------------------------------------------------------------------------------------------------!
650	!
651	! Description:
652	! ------------
653	!> Subroutine for averaging 3D data
654	!--------------------------------------------------------------------------------------------------!
655	SUBROUTINE pcm_3d_data_averaging( mode, variable )
656
657	CHARACTER (LEN=*) :: mode !<
658	CHARACTER (LEN=*) :: variable !<
659
660	INTEGER(iwp) :: i !<
661	INTEGER(iwp) :: j !<
662	INTEGER(iwp) :: k !<
663
664
665	IF ( mode == 'allocate' ) THEN
666
667	SELECT CASE ( TRIM( variable ) )
668
669	CASE ( 'pcm_heatrate' )
670	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
671	ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
672	ENDIF
673	pcm_heatrate_av = 0.0_wp
674
675
676	CASE ( 'pcm_latentrate' )
677	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
678	ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
679	ENDIF
680	pcm_latentrate_av = 0.0_wp
681
682
683	CASE ( 'pcm_transpirationrate' )
684	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
685	ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
686	ENDIF
687	pcm_transpirationrate_av = 0.0_wp
688
689	CASE DEFAULT
690	CONTINUE
691
692	END SELECT
693
694	ELSE IF ( mode == 'sum' ) THEN
695
696	SELECT CASE ( TRIM( variable ) )
697
698	CASE ( 'pcm_heatrate' )
699	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
700	DO i = nxl, nxr
701	DO j = nys, nyn
702	IF ( pch_index_ji(j,i) /= 0 ) THEN
703	DO k = 0, pch_index_ji(j,i)
704	pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + &
705	pcm_heating_rate(k,j,i)
706	ENDDO
707	ENDIF
708	ENDDO
709	ENDDO
710	ENDIF
711
712
713	CASE ( 'pcm_latentrate' )
714	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
715	DO i = nxl, nxr
716	DO j = nys, nyn
717	IF ( pch_index_ji(j,i) /= 0 ) THEN
718	DO k = 0, pch_index_ji(j,i)
719	pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + &
720	pcm_latent_rate(k,j,i)
721	ENDDO
722	ENDIF
723	ENDDO
724	ENDDO
725	ENDIF
726
727
728	CASE ( 'pcm_transpirationrate' )
729	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
730	DO i = nxl, nxr
731	DO j = nys, nyn
732	IF ( pch_index_ji(j,i) /= 0 ) THEN
733	DO k = 0, pch_index_ji(j,i)
734	pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + &
735	pcm_transpiration_rate(k,j,i)
736	ENDDO
737	ENDIF
738	ENDDO
739	ENDDO
740	ENDIF
741
742	CASE DEFAULT
743	CONTINUE
744
745	END SELECT
746
747	ELSE IF ( mode == 'average' ) THEN
748
749	SELECT CASE ( TRIM( variable ) )
750
751	CASE ( 'pcm_heatrate' )
752	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
753	DO i = nxlg, nxrg
754	DO j = nysg, nyng
755	IF ( pch_index_ji(j,i) /= 0 ) THEN
756	DO k = 0, pch_index_ji(j,i)
757	pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) &
758	/ REAL( average_count_3d, KIND=wp )
759	ENDDO
760	ENDIF
761	ENDDO
762	ENDDO
763	ENDIF
764
765
766	CASE ( 'pcm_latentrate' )
767	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
768	DO i = nxlg, nxrg
769	DO j = nysg, nyng
770	IF ( pch_index_ji(j,i) /= 0 ) THEN
771	DO k = 0, pch_index_ji(j,i)
772	pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) &
773	/ REAL( average_count_3d, KIND=wp )
774	ENDDO
775	ENDIF
776	ENDDO
777	ENDDO
778	ENDIF
779
780
781	CASE ( 'pcm_transpirationrate' )
782	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
783	DO i = nxlg, nxrg
784	DO j = nysg, nyng
785	IF ( pch_index_ji(j,i) /= 0 ) THEN
786	DO k = 0, pch_index_ji(j,i)
787	pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) &
788	/ REAL( average_count_3d, KIND=wp )
789	ENDDO
790	ENDIF
791	ENDDO
792	ENDDO
793	ENDIF
794
795	END SELECT
796
797	ENDIF
798
799	END SUBROUTINE pcm_3d_data_averaging
800
801	!--------------------------------------------------------------------------------------------------!
802	!
803	! Description:
804	! ------------
805	!> Subroutine defining 3D output variables.
806	!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning that 3D output is
807	!> relative to the model surface now rather than at the actual grid point where the plant canopy is
808	!> located.
809	!--------------------------------------------------------------------------------------------------!
810	SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
811
812	CHARACTER (LEN=*) :: variable !< treated variable
813
814	INTEGER(iwp) :: av !< flag indicating instantaneous or averaged data output
815	INTEGER(iwp) :: i !< grid index x-direction
816	INTEGER(iwp) :: j !< grid index y-direction
817	INTEGER(iwp) :: k !< grid index z-direction
818	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
819	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
820
821	LOGICAL :: found !< flag indicating if variable is found
822
823	REAL(wp) :: fill_value !< fill value
824
825	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< data output array
826
827
828	found = .TRUE.
829
830	local_pf = REAL( fill_value, KIND = 4 )
831
832	SELECT CASE ( TRIM( variable ) )
833	!
834	!-- Note, to save memory arrays for heating are allocated from 0:pch_index.
835	!-- Thus, output must be relative to these array indices. Further, check whether the output is
836	!-- within the vertical output range, i.e. nzb_do:nzt_do, which is necessary as local_pf is only
837	!-- allocated for this index space. Note, plant-canopy output has a separate vertical output
838	!-- coordinate zlad, so that output is mapped down to the surface.
839	CASE ( 'pcm_heatrate' )
840	IF ( av == 0 ) THEN
841	DO i = nxl, nxr
842	DO j = nys, nyn
843	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
844	local_pf(i,j,k) = pcm_heating_rate(k,j,i)
845	ENDDO
846	ENDDO
847	ENDDO
848	ELSE
849	DO i = nxl, nxr
850	DO j = nys, nyn
851	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
852	local_pf(i,j,k) = pcm_heatrate_av(k,j,i)
853	ENDDO
854	ENDDO
855	ENDDO
856	ENDIF
857
858	CASE ( 'pcm_latentrate' )
859	IF ( av == 0 ) THEN
860	DO i = nxl, nxr
861	DO j = nys, nyn
862	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
863	local_pf(i,j,k) = pcm_latent_rate(k,j,i)
864	ENDDO
865	ENDDO
866	ENDDO
867	ELSE
868	DO i = nxl, nxr
869	DO j = nys, nyn
870	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
871	local_pf(i,j,k) = pcm_latentrate_av(k,j,i)
872	ENDDO
873	ENDDO
874	ENDDO
875	ENDIF
876
877	CASE ( 'pcm_transpirationrate' )
878	IF ( av == 0 ) THEN
879	DO i = nxl, nxr
880	DO j = nys, nyn
881	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
882	local_pf(i,j,k) = pcm_transpiration_rate(k,j,i)
883	ENDDO
884	ENDDO
885	ENDDO
886	ELSE
887	DO i = nxl, nxr
888	DO j = nys, nyn
889	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
890	local_pf(i,j,k) = pcm_transpirationrate_av(k,j,i)
891	ENDDO
892	ENDDO
893	ENDDO
894	ENDIF
895
896	CASE ( 'pcm_lad' )
897	IF ( av == 0 ) THEN
898	DO i = nxl, nxr
899	DO j = nys, nyn
900	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
901	local_pf(i,j,k) = lad_s(k,j,i)
902	ENDDO
903	ENDDO
904	ENDDO
905	ENDIF
906
907	CASE ( 'pcm_bad' )
908	IF ( av == 0 ) THEN
909	DO i = nxl, nxr
910	DO j = nys, nyn
911	DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
912	local_pf(i,j,k) = bad_s(k,j,i)
913	ENDDO
914	ENDDO
915	ENDDO
916	ENDIF
917
918	CASE DEFAULT
919	found = .FALSE.
920
921	END SELECT
922
923	END SUBROUTINE pcm_data_output_3d
924
925	!--------------------------------------------------------------------------------------------------!
926	!
927	! Description:
928	! ------------
929	!> Subroutine defining appropriate grid for netcdf variables.
930	!> It is called from subroutine netcdf.
931	!--------------------------------------------------------------------------------------------------!
932	SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
933
934	CHARACTER (LEN=*), INTENT(IN) :: var !<
935	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
936	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
937	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
938
939	LOGICAL, INTENT(OUT) :: found !<
940
941
942	found = .TRUE.
943
944	!
945	!-- Check for the grid. zpc is zu(nzb:nzb+pch_index)
946	SELECT CASE ( TRIM( var ) )
947
948	CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate' )
949	grid_x = 'x'
950	grid_y = 'y'
951	grid_z = 'zpc'
952
953	CASE DEFAULT
954	found = .FALSE.
955	grid_x = 'none'
956	grid_y = 'none'
957	grid_z = 'none'
958	END SELECT
959
960	END SUBROUTINE pcm_define_netcdf_grid
961
962
963	!--------------------------------------------------------------------------------------------------!
964	! Description:
965	! ------------
966	!> Header output for plant canopy model
967	!--------------------------------------------------------------------------------------------------!
968	SUBROUTINE pcm_header ( io )
969
970	CHARACTER (LEN=10) :: coor_chr !<
971	CHARACTER (LEN=86) :: coordinates !<
972	CHARACTER (LEN=86) :: gradients !<
973	CHARACTER (LEN=86) :: leaf_area_density !<
974	CHARACTER (LEN=86) :: slices !<
975
976	INTEGER(iwp) :: i !<
977	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
978	INTEGER(iwp) :: k !<
979
980	REAL(wp) :: canopy_height !< canopy height (in m)
981
982
983	canopy_height = zw(pch_index)
984
985	WRITE( io, 1 ) canopy_mode, canopy_height, pch_index, canopy_drag_coeff
986	IF ( passive_scalar ) THEN
987	WRITE( io, 2 ) leaf_scalar_exch_coeff, leaf_surface_conc
988	ENDIF
989
990	!
991	!-- Heat flux at the top of vegetation
992	WRITE( io, 3 ) cthf
993
994	!
995	!-- Leaf area density profile, calculated either from given vertical gradients or from beta
996	!-- probability density function.
997	IF ( .NOT. calc_beta_lad_profile ) THEN
998
999	! Building output strings, starting with surface value
1000	WRITE( leaf_area_density, '(F7.4)' ) lad_surface
1001	gradients = '------'
1002	slices = ' 0'
1003	coordinates = ' 0.0'
1004	DO i = 1, UBOUND( lad_vertical_gradient_level_ind, DIM=1 )
1005	IF ( lad_vertical_gradient_level_ind(i) /= -9999 ) THEN
1006
1007	WRITE( coor_chr, '(F7.2)' ) lad(lad_vertical_gradient_level_ind(i))
1008	leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
1009
1010	WRITE( coor_chr, '(F7.2)' ) lad_vertical_gradient(i)
1011	gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
1012
1013	WRITE( coor_chr, '(I7)' ) lad_vertical_gradient_level_ind(i)
1014	slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
1015
1016	WRITE( coor_chr, '(F7.1)' ) lad_vertical_gradient_level(i)
1017	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
1018	ELSE
1019	EXIT
1020	ENDIF
1021	ENDDO
1022
1023	WRITE( io, 4 ) TRIM( coordinates ), TRIM( leaf_area_density ), TRIM( gradients ), &
1024	TRIM( slices )
1025
1026	ELSE
1027
1028	WRITE( leaf_area_density, '(F7.4)' ) lad_surface
1029	coordinates = ' 0.0'
1030
1031	DO k = 1, pch_index
1032
1033	WRITE( coor_chr,'(F7.2)' ) lad(k)
1034	leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
1035
1036	WRITE(coor_chr,'(F7.1)') zu(k)
1037	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
1038
1039	ENDDO
1040
1041	WRITE( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ), alpha_lad, beta_lad, lai_beta
1042
1043	ENDIF
1044
1045	1 FORMAT (/ /' Vegetation canopy (drag) model:' / ' ------------------------------' // &
1046	' Canopy mode: ', A / ' Canopy height: ', F6.2, 'm (',I4,' grid points)' / &
1047	' Leaf drag coefficient: ', F6.2 /)
1048	2 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / &
1049	' Scalar concentration at leaf surfaces in kg/m**3: ', F6.2 /)
1050	3 FORMAT ( ' Predefined constant heatflux at the top of the vegetation: ', F6.2, ' K m/s')
1051	4 FORMAT (/ ' Characteristic levels of the leaf area density:' // &
1052	' Height: ', A, ' m' / &
1053	' Leaf area density: ', A, ' m2/m3' / &
1054	' Gradient: ', A, ' m2/m4' / &
1055	' Gridpoint: ', A )
1056	5 FORMAT (//' Characteristic levels of the leaf area density and coefficients:' // &
1057	' Height: ', A, ' m' / &
1058	' Leaf area density: ', A, ' m2/m3' / &
1059	' Coefficient alpha: ',F6.2 / &
1060	' Coefficient beta: ',F6.2 / &
1061	' Leaf area index: ',F6.2,' m2/m2' /)
1062
1063	END SUBROUTINE pcm_header
1064
1065
1066	!--------------------------------------------------------------------------------------------------!
1067	! Description:
1068	! ------------
1069	!> Initialization of the plant canopy model
1070	!--------------------------------------------------------------------------------------------------!
1071	SUBROUTINE pcm_init
1072
1073	USE exchange_horiz_mod, &
1074	ONLY: exchange_horiz
1075
1076	INTEGER(iwp) :: i !< running index
1077	INTEGER(iwp) :: j !< running index
1078	INTEGER(iwp) :: k !< running index
1079	INTEGER(iwp) :: m !< running index
1080
1081	LOGICAL :: lad_on_top = .FALSE. !< dummy flag to indicate that LAD is defined on a building roof
1082	LOGICAL :: bad_on_top = .FALSE. !< dummy flag to indicate that BAD is defined on a building roof
1083
1084	REAL(wp) :: canopy_height !< canopy height for lad-profile construction
1085	REAL(wp) :: gradient !< gradient for lad-profile construction
1086	REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction
1087	REAL(wp) :: lad_max !< maximum LAD value in the model domain, used to perform a check
1088
1089
1090	IF ( debug_output ) CALL debug_message( 'pcm_init', 'start' )
1091	!
1092	!-- Allocate one-dimensional arrays for the computation of the leaf area density (lad) profile
1093	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
1094	lad = 0.0_wp
1095	pre_lad = 0.0_wp
1096
1097	!
1098	!-- Set flag that indicates that the lad-profile shall be calculated by using a beta probability
1099	!-- density function
1100	IF ( alpha_lad /= 9999999.9_wp .AND. beta_lad /= 9999999.9_wp ) THEN
1101	calc_beta_lad_profile = .TRUE.
1102	ENDIF
1103
1104
1105	!
1106	!-- Compute the profile of leaf area density used in the plant canopy model. The profile can
1107	!-- either be constructed from prescribed vertical gradients of the leaf area density or by using
1108	!-- a beta probability density function (see e.g. Markkanen et al., 2003: Boundary-Layer
1109	!-- Meteorology, 106, 437-459)
1110	IF ( .NOT. calc_beta_lad_profile ) THEN
1111
1112	!
1113	!-- Use vertical gradients for lad-profile construction
1114	i = 1
1115	gradient = 0.0_wp
1116
1117	lad(0) = lad_surface
1118	lad_vertical_gradient_level_ind(1) = 0
1119
1120	DO k = 1, pch_index
1121	IF ( i < 11 ) THEN
1122	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
1123	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
1124	gradient = lad_vertical_gradient(i)
1125	lad_vertical_gradient_level_ind(i) = k - 1
1126	i = i + 1
1127	ENDIF
1128	ENDIF
1129	IF ( gradient /= 0.0_wp ) THEN
1130	IF ( k /= 1 ) THEN
1131	lad(k) = lad(k-1) + dzu(k) * gradient
1132	ELSE
1133	lad(k) = lad_surface + dzu(k) * gradient
1134	ENDIF
1135	ELSE
1136	lad(k) = lad(k-1)
1137	ENDIF
1138	ENDDO
1139
1140	!
1141	!-- In case of no given leaf area density gradients, choose a vanishing gradient. This
1142	!-- information is used for the HEADER and the RUN_CONTROL file.
1143	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
1144	lad_vertical_gradient_level(1) = 0.0_wp
1145	ENDIF
1146
1147	ELSE
1148
1149	!
1150	!-- Use beta function for lad-profile construction
1151	int_bpdf = 0.0_wp
1152	canopy_height = zw(pch_index)
1153
1154	DO k = 0, pch_index
1155	int_bpdf = int_bpdf + &
1156	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
1157	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) &
1158	* ( ( zw(k+1)-zw(k) ) / canopy_height ) )
1159	ENDDO
1160
1161	!
1162	!-- Preliminary lad profile (defined on w-grid)
1163	DO k = 0, pch_index
1164	pre_lad(k) = lai_beta * &
1165	( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) &
1166	* ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) &
1167	/ int_bpdf &
1168	) / canopy_height
1169	ENDDO
1170
1171	!
1172	!-- Final lad profile (defined on scalar-grid level, since most prognostic quantities are
1173	!-- defined there, hence, less interpolation is required when calculating the canopy
1174	!-- tendencies)
1175	lad(0) = pre_lad(0)
1176	DO k = 1, pch_index
1177	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
1178	ENDDO
1179
1180	ENDIF
1181
1182	!
1183	!-- Allocate 3D-array for the leaf-area density (lad_s) as well as for basal-area densitiy
1184	!-- (bad_s). Note, by default bad_s is zero.
1185	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1186	ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1187	bad_s = 0.0_wp
1188	!
1189	!-- Initialization of the canopy coverage in the model domain:
1190	!-- Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which fully covers
1191	!-- the domain surface
1192	SELECT CASE ( TRIM( canopy_mode ) )
1193
1194	CASE ( 'homogeneous' )
1195
1196	DO i = nxlg, nxrg
1197	DO j = nysg, nyng
1198	lad_s(:,j,i) = lad(:)
1199	ENDDO
1200	ENDDO
1201
1202	CASE ( 'read_from_file' )
1203	!
1204	!-- Read plant canopy
1205	IF ( input_pids_static ) THEN
1206	!
1207	!-- Open the static input file
1208	#if defined( __netcdf )
1209	CALL open_read_file( TRIM( input_file_static ) // &
1210	TRIM( coupling_char ), &
1211	pids_id )
1212
1213	CALL inquire_num_variables( pids_id, num_var_pids )
1214	!
1215	!-- Allocate memory to store variable names and read them
1216	ALLOCATE( vars_pids(1:num_var_pids) )
1217	CALL inquire_variable_names( pids_id, vars_pids )
1218	!
1219	!-- Read leaf area density - resolved vegetation
1220	IF ( check_existence( vars_pids, 'lad' ) ) THEN
1221	leaf_area_density_f%from_file = .TRUE.
1222	CALL get_attribute( pids_id, char_fill, &
1223	leaf_area_density_f%fill, &
1224	.FALSE., 'lad' )
1225	!
1226	!-- Inquire number of vertical vegetation layer
1227	CALL get_dimension_length( pids_id, &
1228	leaf_area_density_f%nz, &
1229	'zlad' )
1230	!
1231	!-- Allocate variable for leaf-area density
1232	ALLOCATE( leaf_area_density_f%var &
1233	(0:leaf_area_density_f%nz-1,nys:nyn,nxl:nxr) )
1234
1235	CALL get_variable( pids_id, 'lad', leaf_area_density_f%var, nxl, nxr, nys, nyn, &
1236	0, leaf_area_density_f%nz-1 )
1237
1238	ELSE
1239	leaf_area_density_f%from_file = .FALSE.
1240	ENDIF
1241	!
1242	!-- Read basal area density - resolved vegetation
1243	IF ( check_existence( vars_pids, 'bad' ) ) THEN
1244	basal_area_density_f%from_file = .TRUE.
1245	CALL get_attribute( pids_id, char_fill, &
1246	basal_area_density_f%fill, &
1247	.FALSE., 'bad' )
1248	!
1249	!-- Inquire number of vertical vegetation layer
1250	CALL get_dimension_length( pids_id, &
1251	basal_area_density_f%nz, &
1252	'zlad' )
1253	!
1254	!-- Allocate variable
1255	ALLOCATE( basal_area_density_f%var &
1256	(0:basal_area_density_f%nz-1,nys:nyn,nxl:nxr) )
1257
1258	CALL get_variable( pids_id, 'bad', basal_area_density_f%var, nxl, nxr, nys, nyn,&
1259	0, basal_area_density_f%nz-1 )
1260	ELSE
1261	basal_area_density_f%from_file = .FALSE.
1262	ENDIF
1263	!
1264	!-- Read root area density - resolved vegetation
1265	IF ( check_existence( vars_pids, 'root_area_dens_r' ) ) THEN
1266	root_area_density_lad_f%from_file = .TRUE.
1267	CALL get_attribute( pids_id, char_fill, &
1268	root_area_density_lad_f%fill, &
1269	.FALSE., 'root_area_dens_r' )
1270	!
1271	!-- Inquire number of vertical soil layers
1272	CALL get_dimension_length( pids_id, &
1273	root_area_density_lad_f%nz, &
1274	'zsoil' )
1275	!
1276	!-- Allocate variable
1277	ALLOCATE( root_area_density_lad_f%var &
1278	(0:root_area_density_lad_f%nz-1,nys:nyn,nxl:nxr) )
1279
1280	CALL get_variable( pids_id, 'root_area_dens_r', root_area_density_lad_f%var, &
1281	nxl, nxr, nys, nyn, 0, root_area_density_lad_f%nz-1 )
1282	ELSE
1283	root_area_density_lad_f%from_file = .FALSE.
1284	ENDIF
1285
1286	DEALLOCATE( vars_pids )
1287	!
1288	!-- Finally, close the input file and deallocate temporary array
1289	CALL close_input_file( pids_id )
1290	#endif
1291	ENDIF
1292
1293	!
1294	!-- Initialize LAD with data from file. If LAD is given in NetCDF file, use these values,
1295	!-- else take LAD profiles from ASCII file.
1296	!-- Please note, in NetCDF file LAD is only given up to the maximum canopy top, indicated
1297	!-- by leaf_area_density_f%nz.
1298	lad_s = 0.0_wp
1299	IF ( leaf_area_density_f%from_file ) THEN
1300	!
1301	!-- Set also pch_index, used to be the upper bound of the vertical loops. Therefore, use
1302	!-- the global top of the canopy layer.
1303	pch_index = leaf_area_density_f%nz - 1
1304
1305	DO i = nxl, nxr
1306	DO j = nys, nyn
1307	DO k = 0, leaf_area_density_f%nz - 1
1308	IF ( leaf_area_density_f%var(k,j,i) /= leaf_area_density_f%fill ) &
1309	lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
1310	ENDDO
1311	!
1312	!-- Check if resolved vegetation is mapped onto buildings.
1313	!-- In general, this is allowed and also meaningful, e.g. when trees carry across
1314	!-- roofs. However, due to the topography filtering, new building grid points can
1315	!-- emerge at locations where also plant canopy is defined. As a result, plant
1316	!-- canopy is mapped on top of roofs, with siginficant impact on the downstream
1317	!-- flow field and the nearby surface radiation. In order to avoid that plant
1318	!-- canopy is mistakenly mapped onto building roofs, check for building grid
1319	!-- points (bit 6) that emerge from the filtering (bit 4) and set LAD to zero at
1320	!-- these artificially created building grid points. This case, an informative
1321	!-- message is given.
1322	IF ( ANY( lad_s(:,j,i) /= 0.0_wp ) .AND. &
1323	ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
1324	ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) ) THEN
1325	lad_s(:,j,i) = 0.0_wp
1326	lad_on_top = .TRUE.
1327	ENDIF
1328	ENDDO
1329	ENDDO
1330	#if defined( __parallel )
1331	CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_on_top, 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr)
1332	#endif
1333	IF ( lad_on_top ) THEN
1334	WRITE( message_string, * ) &
1335	'Resolved plant-canopy is defined on top of an ' // &
1336	'artificially created building grid point(s) '// &
1337	'the filtering) - LAD/BAD profile is omitted at this / ' //&
1338	'these grid point(s).'
1339	CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
1340	ENDIF
1341	CALL exchange_horiz( lad_s, nbgp )
1342	!
1343	! ASCII file
1344	!-- Initialize canopy parameters canopy_drag_coeff, leaf_scalar_exch_coeff,
1345	!-- leaf_surface_conc from file which contains complete 3D data (separate vertical profiles
1346	!-- for each location).
1347	ELSE
1348	CALL pcm_read_plant_canopy_3d
1349	ENDIF
1350	!
1351	!-- Initialize LAD with data from file. If LAD is given in NetCDF file, use these values,
1352	!-- else take LAD profiles from ASCII file.
1353	!-- Please note, in NetCDF file LAD is only given up to the maximum canopy top, indicated
1354	!-- by basal_area_density_f%nz.
1355	bad_s = 0.0_wp
1356	IF ( basal_area_density_f%from_file ) THEN
1357	!
1358	!-- Set also pch_index, used to be the upper bound of the vertical loops. Therefore, use
1359	!-- the global top of the canopy layer.
1360	pch_index = basal_area_density_f%nz - 1
1361
1362	DO i = nxl, nxr
1363	DO j = nys, nyn
1364	DO k = 0, basal_area_density_f%nz - 1
1365	IF ( basal_area_density_f%var(k,j,i) /= basal_area_density_f%fill ) &
1366	bad_s(k,j,i) = basal_area_density_f%var(k,j,i)
1367	ENDDO
1368	!
1369	!-- Check if resolved vegetation is mapped onto buildings.
1370	!-- Please see comment for leaf_area density
1371	IF ( ANY( bad_s(:,j,i) /= 0.0_wp ) .AND. &
1372	ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
1373	ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) ) THEN
1374	bad_s(:,j,i) = 0.0_wp
1375	bad_on_top = .TRUE.
1376	ENDIF
1377	ENDDO
1378	ENDDO
1379	#if defined( __parallel )
1380	CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr)
1381	#endif
1382	IF ( bad_on_top ) THEN
1383	WRITE( message_string, * ) &
1384	'Resolved plant-canopy is defined on top of an ' // &
1385	'artificially created building grid point(s) '// &
1386	'the filtering) - LAD/BAD profile is omitted at this / ' //&
1387	'these grid point(s).'
1388	CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
1389	ENDIF
1390	CALL exchange_horiz( bad_s, nbgp )
1391	ENDIF
1392
1393	CASE DEFAULT
1394	!
1395	!-- The DEFAULT case is reached either if the parameter canopy mode contains a wrong
1396	!-- character string or if the user has coded a special case in the user interface.
1397	!-- There, the subroutine user_init_plant_canopy checks which of these two conditions
1398	!-- applies.
1399	CALL user_init_plant_canopy
1400
1401	END SELECT
1402	!
1403	!-- Check that at least one grid point has an LAD /= 0, else this may cause errors in the
1404	!-- radiation model.
1405	lad_max = MAXVAL( lad_s )
1406	#if defined( __parallel )
1407	CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX, comm2d, ierr)
1408	#endif
1409	IF ( lad_max <= 0.0_wp ) THEN
1410	message_string = 'Plant-canopy model is switched-on but no ' // &
1411	'plant canopy is present in the model domain.'
1412	CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
1413	ENDIF
1414
1415	!
1416	!-- Initialize 2D index array indicating canopy top index.
1417	ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
1418	pch_index_ji = 0
1419
1420	DO i = nxlg, nxrg
1421	DO j = nysg, nyng
1422	DO k = 0, pch_index
1423	IF ( lad_s(k,j,i) /= 0.0_wp .OR. bad_s(k,j,i) /= 0.0_wp ) pch_index_ji(j,i) = k
1424	ENDDO
1425	!
1426	!-- Check whether topography and local vegetation on top exceed height of the model domain.
1427	IF ( topo_top_ind(j,i,0) + pch_index_ji(j,i) >= nzt + 1 ) THEN
1428	message_string = 'Local vegetation height on top of ' // &
1429	'topography exceeds height of model domain.'
1430	CALL message( 'pcm_init', 'PA0674', 2, 2, myid, 6, 0 )
1431	ENDIF
1432
1433	ENDDO
1434	ENDDO
1435	!
1436	!-- Calculate global pch_index value (index of top of plant canopy from ground)
1437	pch_index = MAXVAL( pch_index_ji )
1438	!
1439	!-- Exchange pch_index from all processors
1440	#if defined( __parallel )
1441	CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr)
1442	#endif
1443	!
1444	!-- Allocation of arrays pcm_heating_rate, pcm_transpiration_rate and pcm_latent_rate
1445	ALLOCATE( pcm_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1446	pcm_heating_rate = 0.0_wp
1447
1448	IF ( humidity ) THEN
1449	ALLOCATE( pcm_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1450	pcm_transpiration_rate = 0.0_wp
1451	ALLOCATE( pcm_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1452	pcm_latent_rate = 0.0_wp
1453	ENDIF
1454	!
1455	!-- Initialization of the canopy heat source distribution due to heating of the canopy layers by
1456	!-- incoming solar radiation, in case that a non-zero
1457	!-- value is set for the canopy top heat flux (cthf), which equals the available net radiation at
1458	!-- canopy top.
1459	!-- The heat source distribution is calculated by a decaying exponential function of the downward
1460	!-- cumulative leaf area index (cum_lai_hf), assuming that the foliage inside the plant canopy is
1461	!-- heated by solar radiation penetrating the canopy layers according to the distribution of net
1462	!-- radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3, 73â96). This approach has
1463	!-- been applied e.g. by Shaw & Schumann (1992; Bound.-Layer Meteorol. 61, 47â64).
1464	!-- When using the radiation_interactions, canopy heating (pcm_heating_rate) and plant canopy
1465	!-- transpiration (pcm_transpiration_rate, pcm_latent_rate) are calculated in the RTM after the
1466	!-- calculation of radiation.
1467	IF ( cthf /= 0.0_wp ) THEN
1468
1469	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1470	!
1471	!-- Piecewise calculation of the cumulative leaf area index by vertical integration of the
1472	!-- leaf area density
1473	cum_lai_hf(:,:,:) = 0.0_wp
1474	DO i = nxlg, nxrg
1475	DO j = nysg, nyng
1476	DO k = pch_index_ji(j,i)-1, 0, -1
1477	IF ( k == pch_index_ji(j,i)-1 ) THEN
1478	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
1479	( 0.5_wp * lad_s(k+1,j,i) * &
1480	( zw(k+1) - zu(k+1) ) ) + &
1481	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
1482	lad_s(k,j,i) ) + &
1483	lad_s(k+1,j,i) ) * &
1484	( zu(k+1) - zw(k) ) )
1485	ELSE
1486	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
1487	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
1488	lad_s(k+1,j,i) ) + &
1489	lad_s(k+1,j,i) ) * &
1490	( zw(k+1) - zu(k+1) ) ) + &
1491	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
1492	lad_s(k,j,i) ) + &
1493	lad_s(k+1,j,i) ) * &
1494	( zu(k+1) - zw(k) ) )
1495	ENDIF
1496	ENDDO
1497	ENDDO
1498	ENDDO
1499
1500	!
1501	!-- In areas with canopy the surface value of the canopy heat flux distribution overrides the
1502	!-- surface heat flux (shf),
1503	!-- Start with default surface type
1504	DO m = 1, surf_def_h(0)%ns
1505	i = surf_def_h(0)%i(m)
1506	j = surf_def_h(0)%j(m)
1507	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1508	surf_def_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
1509	ENDDO
1510	!
1511	!-- Natural surfaces
1512	DO m = 1, surf_lsm_h(0)%ns
1513	i = surf_lsm_h(0)%i(m)
1514	j = surf_lsm_h(0)%j(m)
1515	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1516	surf_lsm_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
1517	ENDDO
1518	!
1519	!-- Urban surfaces
1520	DO m = 1, surf_usm_h(0)%ns
1521	i = surf_usm_h(0)%i(m)
1522	j = surf_usm_h(0)%j(m)
1523	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) &
1524	surf_usm_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
1525	ENDDO
1526	!
1527	!
1528	!-- Calculation of the heating rate (K/s) within the different layers of the plant canopy.
1529	!-- Calculation is only necessary in areas covered with canopy.
1530	!-- Within the different canopy layers the plant-canopy heating rate (pcm_heating_rate) is
1531	!-- calculated as the vertical divergence of the canopy heat fluxes at the top and bottom of
1532	!-- the respective layer.
1533	DO i = nxlg, nxrg
1534	DO j = nysg, nyng
1535	DO k = 1, pch_index_ji(j,i)
1536	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) THEN
1537	pcm_heating_rate(k,j,i) = cthf * &
1538	( EXP( -ext_coef * cum_lai_hf(k,j,i) ) - &
1539	EXP( -ext_coef * cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1540	ENDIF
1541	ENDDO
1542	ENDDO
1543	ENDDO
1544
1545	ENDIF
1546
1547	IF ( debug_output ) CALL debug_message( 'pcm_init', 'end' )
1548
1549	END SUBROUTINE pcm_init
1550
1551
1552	!--------------------------------------------------------------------------------------------------!
1553	! Description:
1554	! ------------
1555	!> Parin for &plant_canopy_parameters for plant canopy model
1556	!--------------------------------------------------------------------------------------------------!
1557	SUBROUTINE pcm_parin
1558
1559	CHARACTER(LEN=100) :: line !< dummy string that contains the current line of the parameter file
1560
1561	INTEGER(iwp) :: io_status !< status after reading the namelist file
1562
1563
1564	NAMELIST /plant_canopy_parameters/ alpha_lad, &
1565	beta_lad, &
1566	canopy_drag_coeff, &
1567	canopy_mode, &
1568	cthf, &
1569	lad_surface, &
1570	lad_type_coef, &
1571	lad_vertical_gradient, &
1572	lad_vertical_gradient_level, &
1573	lai_beta, &
1574	leaf_scalar_exch_coeff, &
1575	leaf_surface_conc, &
1576	pch_index, &
1577	plant_canopy_transpiration
1578
1579	!
1580	!-- Move to the beginning of the namelist file and try to find and read the user-defined namelist
1581	!-- plant_canopy_parameters.
1582	REWIND( 11 )
1583	READ( 11, plant_canopy_parameters, IOSTAT=io_status )
1584
1585	!
1586	!-- Action depending on the READ status
1587	IF ( io_status == 0 ) THEN
1588	!
1589	!-- plant_canopy_parameters namelist was found and read correctly. Set flag that indicates that
1590	!-- the plant-canopy model is switched on.
1591	plant_canopy = .TRUE.
1592
1593	ELSEIF ( io_status > 0 ) THEN
1594	!
1595	!-- plant_canopy_parameters namelist was found but contained errors. Print an error message
1596	!-- including the line that caused the problem.
1597	BACKSPACE( 11 )
1598	READ( 11 , '(A)' ) line
1599	CALL parin_fail_message( 'plant_canopy_parameters', line )
1600
1601	ENDIF
1602
1603	END SUBROUTINE pcm_parin
1604
1605
1606	!--------------------------------------------------------------------------------------------------!
1607	! Description:
1608	! ------------
1609	!
1610	!> Loads 3D plant canopy data from file. File format is as follows:
1611	!>
1612	!> num_levels
1613	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1614	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1615	!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1616	!> ...
1617	!>
1618	!> i.e. first line determines number of levels and further lines represent plant canopy data, one
1619	!> line per column and variable. In each data line, dtype represents variable to be set:
1620	!>
1621	!> dtype=1: leaf area density (lad_s)
1622	!> dtype=2....n: some additional plant canopy input data quantity
1623	!>
1624	!> Zeros are added automatically above num_levels until top of domain. Any non-specified (x,y)
1625	!> columns have zero values as default.
1626	!--------------------------------------------------------------------------------------------------!
1627	SUBROUTINE pcm_read_plant_canopy_3d
1628
1629	USE exchange_horiz_mod, &
1630	ONLY: exchange_horiz
1631
1632	INTEGER(iwp) :: dtype !< type of input data (1=lad)
1633	INTEGER(iwp) :: i !< running index
1634	INTEGER(iwp) :: j !< running index
1635	INTEGER(iwp) :: kk !<
1636	INTEGER(iwp) :: nzp !< number of vertical layers of plant canopy
1637	INTEGER(iwp) :: nzpltop !<
1638	INTEGER(iwp) :: nzpl !<
1639	INTEGER(iwp) :: pctype !< type of plant canopy (deciduous,non-deciduous,...)
1640
1641	REAL(wp), DIMENSION(:), ALLOCATABLE :: col !< vertical column of input data
1642
1643	!
1644	!-- Initialize lad_s array
1645	lad_s = 0.0_wp
1646
1647	!
1648	!-- Open and read plant canopy input data
1649	OPEN( 152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ), ACCESS='SEQUENTIAL', &
1650	ACTION='READ', STATUS='OLD', FORM='FORMATTED', ERR=515 )
1651	READ( 152, *, ERR=516, END=517 ) nzp !< read first line = number of vertical layers
1652	nzpltop = MIN(nzt+1, nzb+nzp-1)
1653	nzpl = nzpltop - nzb + 1 !< no. of layers to assign
1654	ALLOCATE( col(0:nzp-1) )
1655
1656	DO
1657	READ( 152, *, ERR=516, END=517 ) dtype, i, j, pctype, col(:)
1658	IF ( i < nxlg .OR. i > nxrg .OR. j < nysg .OR. j > nyng ) CYCLE
1659
1660	SELECT CASE ( dtype )
1661	CASE( 1 ) !< leaf area density
1662	!
1663	!-- This is just the pure canopy layer assumed to be grounded to a flat domain surface.
1664	!-- At locations where plant canopy sits on top of any kind of topography, the vertical
1665	!-- plant column must be "lifted", which is done in SUBROUTINE pcm_tendency.
1666	IF ( pctype < 0 .OR. pctype > 10 ) THEN !< incorrect plant canopy type
1667	WRITE( message_string, * ) 'Incorrect type of plant canopy. ' // &
1668	'Allowed values 0 <= pctype <= 10, ' // &
1669	'but pctype is ', pctype
1670	CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1671	ENDIF
1672	kk = topo_top_ind(j,i,0)
1673	lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
1674	CASE DEFAULT
1675	WRITE( message_string, '(a,i2,a)' ) &
1676	'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1677	CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1678	END SELECT
1679	ENDDO
1680
1681	515 message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1682	CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
1683
1684	516 message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1685	CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1686
1687	517 CLOSE( 152 )
1688	DEALLOCATE( col )
1689
1690	CALL exchange_horiz( lad_s, nbgp )
1691
1692	END SUBROUTINE pcm_read_plant_canopy_3d
1693
1694	!--------------------------------------------------------------------------------------------------!
1695	! Description:
1696	! ------------
1697	!> Read module-specific global restart data (Fortran binary format).
1698	!--------------------------------------------------------------------------------------------------!
1699	SUBROUTINE pcm_rrd_global_ftn( found )
1700
1701	LOGICAL, INTENT(OUT) :: found
1702
1703
1704	found = .TRUE.
1705
1706	SELECT CASE ( restart_string(1:length) )
1707
1708	CASE ( 'pch_index' )
1709	READ( 13 ) pch_index
1710
1711	CASE DEFAULT
1712
1713	found = .FALSE.
1714
1715	END SELECT
1716
1717	END SUBROUTINE pcm_rrd_global_ftn
1718
1719	!--------------------------------------------------------------------------------------------------!
1720	! Description:
1721	! ------------
1722	!> Read module-specific global restart data (MPI-IO).
1723	!--------------------------------------------------------------------------------------------------!
1724	SUBROUTINE pcm_rrd_global_mpi
1725
1726	CALL rrd_mpi_io( 'pch_index', pch_index )
1727
1728	END SUBROUTINE pcm_rrd_global_mpi
1729
1730
1731	!--------------------------------------------------------------------------------------------------!
1732	! Description:
1733	! ------------
1734	!> Read module-specific local restart data arrays (Fortran binary format).
1735	!--------------------------------------------------------------------------------------------------!
1736	SUBROUTINE pcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, &
1737	nyn_on_file, nysf, nysc, nys_on_file, found )
1738
1739	INTEGER(iwp) :: k !<
1740	INTEGER(iwp) :: nxl_on_file !<
1741	INTEGER(iwp) :: nxlc !<
1742	INTEGER(iwp) :: nxlf !<
1743	INTEGER(iwp) :: nxr_on_file !<
1744	INTEGER(iwp) :: nxrc !<
1745	INTEGER(iwp) :: nxrf !<
1746	INTEGER(iwp) :: nyn_on_file !<
1747	INTEGER(iwp) :: nync !<
1748	INTEGER(iwp) :: nynf !<
1749	INTEGER(iwp) :: nys_on_file !<
1750	INTEGER(iwp) :: nysc !<
1751	INTEGER(iwp) :: nysf !<
1752
1753	LOGICAL, INTENT(OUT) :: found
1754
1755	REAL(wp), DIMENSION( 0:pch_index, &
1756	nys_on_file-nbgp:nyn_on_file+nbgp, &
1757	nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2 !< temporary 3D array for entire vertical
1758	!< extension of canopy layer
1759	found = .TRUE.
1760
1761
1762	SELECT CASE ( restart_string(1:length) )
1763
1764	CASE ( 'pcm_heatrate_av' )
1765	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
1766	ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1767	pcm_heatrate_av = 0.0_wp
1768	ENDIF
1769	IF ( k == 1 ) READ( 13 ) tmp_3d2
1770	pcm_heatrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1771	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1772
1773	CASE ( 'pcm_latentrate_av' )
1774	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
1775	ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1776	pcm_latentrate_av = 0.0_wp
1777	ENDIF
1778	IF ( k == 1 ) READ( 13 ) tmp_3d2
1779	pcm_latentrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1780	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1781
1782	CASE ( 'pcm_transpirationrate_av' )
1783	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
1784	ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1785	pcm_transpirationrate_av = 0.0_wp
1786	ENDIF
1787	IF ( k == 1 ) READ( 13 ) tmp_3d2
1788	pcm_transpirationrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1789	tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1790
1791	CASE DEFAULT
1792
1793	found = .FALSE.
1794
1795	END SELECT
1796
1797	END SUBROUTINE pcm_rrd_local_ftn
1798
1799
1800	!--------------------------------------------------------------------------------------------------!
1801	! Description:
1802	! ------------
1803	!> Read module-specific local restart data arrays (MPI-IO).
1804	!--------------------------------------------------------------------------------------------------!
1805	SUBROUTINE pcm_rrd_local_mpi
1806
1807	IMPLICIT NONE
1808
1809	LOGICAL :: array_found !<
1810
1811	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d !< temporary array to store pcm data with
1812	!< non-standard vertical index bounds
1813
1814	!
1815	!-- Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
1816	!-- full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored after reading.
1817	CALL rd_mpi_io_check_array( 'pcm_heatrate_av' , found = array_found )
1818	IF ( array_found ) THEN
1819	IF ( .NOT. ALLOCATED( pcm_heatrate_av ) ) THEN
1820	ALLOCATE( pcm_heatrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1821	ENDIF
1822	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1823	CALL rrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
1824	pcm_heatrate_av = tmp_3d(nzb:pch_index,:,:)
1825	DEALLOCATE( tmp_3d )
1826	ENDIF
1827
1828	CALL rd_mpi_io_check_array( 'pcm_latentrate_av' , found = array_found )
1829	IF ( array_found ) THEN
1830	IF ( .NOT. ALLOCATED( pcm_latentrate_av ) ) THEN
1831	ALLOCATE( pcm_latentrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1832	ENDIF
1833	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1834	CALL rrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
1835	pcm_latentrate_av = tmp_3d(nzb:pch_index,:,:)
1836	DEALLOCATE( tmp_3d )
1837	ENDIF
1838
1839	CALL rd_mpi_io_check_array( 'pcm_transpirationrate_av' , found = array_found )
1840	IF ( array_found ) THEN
1841	IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) ) THEN
1842	ALLOCATE( pcm_transpirationrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1843	ENDIF
1844	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1845	CALL rrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
1846	pcm_transpirationrate_av = tmp_3d(nzb:pch_index,:,:)
1847	DEALLOCATE( tmp_3d )
1848	ENDIF
1849
1850	END SUBROUTINE pcm_rrd_local_mpi
1851
1852
1853	!--------------------------------------------------------------------------------------------------!
1854	! Description:
1855	! ------------
1856	!> Calculation of the tendency terms, accounting for the effect of the plant canopy on momentum and
1857	!> scalar quantities.
1858	!>
1859	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 (defined on scalar grid), as
1860	!> initialized in subroutine pcm_init.
1861	!> The lad on the w-grid is vertically interpolated from the surrounding lad_s. The upper boundary
1862	!> of the canopy is defined on the w-grid at k = pch_index. Here, the lad is zero.
1863	!>
1864	!> The canopy drag must be limited (previously accounted for by calculation of a limiting canopy
1865	!> timestep for the determination of the maximum LES timestep in subroutine timestep), since it is
1866	!> physically impossible that the canopy drag alone can locally change the sign of a velocity
1867	!> component. This limitation is realized by calculating preliminary tendencies and velocities. It
1868	!> is subsequently checked if the preliminary new velocity has a different sign than the current
1869	!> velocity. If so, the tendency is limited in a way that the velocity can at maximum be reduced to
1870	!> zero by the canopy drag.
1871	!>
1872	!>
1873	!> Call for all grid points
1874	!--------------------------------------------------------------------------------------------------!
1875	SUBROUTINE pcm_tendency( component )
1876
1877	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
1878	INTEGER(iwp) :: i !< running index
1879	INTEGER(iwp) :: j !< running index
1880	INTEGER(iwp) :: k !< running index
1881	INTEGER(iwp) :: kk !< running index for flat lad arrays
1882
1883	LOGICAL :: building_edge_e !< control flag indicating an eastward-facing building edge
1884	LOGICAL :: building_edge_n !< control flag indicating a north-facing building edge
1885	LOGICAL :: building_edge_s !< control flag indicating a south-facing building edge
1886	LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge
1887
1888	REAL(wp) :: bad_local !< local bad value
1889	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
1890	REAL(wp) :: lad_local !< local lad value
1891	REAL(wp) :: pre_tend !< preliminary tendency
1892	REAL(wp) :: pre_u !< preliminary u-value
1893	REAL(wp) :: pre_v !< preliminary v-value
1894	REAL(wp) :: pre_w !< preliminary w-value
1895
1896
1897	ddt_3d = 1.0_wp / dt_3d
1898
1899	!
1900	!-- Compute drag for the three velocity components and the SGS-TKE:
1901	SELECT CASE ( component )
1902
1903	!
1904	!-- u-component
1905	CASE ( 1 )
1906	DO i = nxlu, nxr
1907	DO j = nys, nyn
1908	!
1909	!-- Set control flags indicating east- and westward-orientated building edges. Note,
1910	!-- building_egde_w is set from the perspective of the potential rooftop grid point,
1911	!-- while building_edge_e is set from the perspective of the non-building grid point.
1912	building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) &
1913	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
1914	building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) ) &
1915	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
1916	!
1917	!-- Determine topography-top index on u-grid
1918	DO k = topo_top_ind(j,i,1)+1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
1919
1920	kk = k - topo_top_ind(j,i,1) !- lad arrays are defined flat
1921	!
1922	!-- In order to create sharp boundaries of the plant canopy, the lad on the u-grid
1923	!-- at index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from
1924	!-- the surrounding lad_s, because this would yield smaller lad at the canopy
1925	!-- boundaries than inside of the canopy.
1926	!-- For the same reason, the lad at the rightmost(i+1)canopy boundary on the
1927	!-- u-grid equals lad_s(k,j,i), which is considered in the next if-statement.
1928	!-- Note, at left-sided building edges this is not applied, here the LAD equals
1929	!-- the LAD at grid point (k,j,i), in order to avoid that LAD is mistakenly mapped
1930	!-- on top of a roof where (usually) no LAD is defined. The same is also valid for
1931	!-- bad_s.
1932	lad_local = lad_s(kk,j,i)
1933	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp &
1934	.AND. .NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1)
1935
1936	bad_local = bad_s(kk,j,i)
1937	IF ( bad_local == 0.0_wp .AND. bad_s(kk,j,i-1) > 0.0_wp &
1938	.AND. .NOT. building_edge_w ) bad_local = bad_s(kk,j,i-1)
1939	!
1940	!-- In order to avoid that LAD is mistakenly considered at right-sided building
1941	!-- edges (here the topography-top index for the u-component at index j,i is still
1942	!-- on the building while the topography top for the scalar isn't), LAD is taken
1943	!-- from grid point (j,i-1). The same is also valid for bad_s.
1944	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp &
1945	.AND. building_edge_e ) lad_local = lad_s(kk,j,i-1)
1946	IF ( bad_local > 0.0_wp .AND. bad_s(kk,j,i-1) == 0.0_wp &
1947	.AND. building_edge_e ) bad_local = bad_s(kk,j,i-1)
1948
1949	pre_tend = 0.0_wp
1950	pre_u = 0.0_wp
1951	!
1952	!-- Calculate preliminary value (pre_tend) of the tendency
1953	pre_tend = - canopy_drag_coeff * &
1954	( lad_local + bad_local ) * &
1955	SQRT( u(k,j,i)**2 + &
1956	( 0.25_wp * ( v(k,j,i-1) + &
1957	v(k,j,i) + &
1958	v(k,j+1,i) + &
1959	v(k,j+1,i-1) ) &
1960	)**2 + &
1961	( 0.25_wp * ( w(k-1,j,i-1) + &
1962	w(k-1,j,i) + &
1963	w(k,j,i-1) + &
1964	w(k,j,i) ) &
1965	)**2 &
1966	) * &
1967	u(k,j,i)
1968
1969	!
1970	!-- Calculate preliminary new velocity, based on pre_tend
1971	pre_u = u(k,j,i) + dt_3d * pre_tend
1972	!
1973	!-- Compare sign of old velocity and new preliminary velocity,
1974	!-- and in case the signs are different, limit the tendency
1975	IF ( SIGN( pre_u,u(k,j,i) ) /= pre_u ) THEN
1976	pre_tend = - u(k,j,i) * ddt_3d
1977	ENDIF
1978	!
1979	!-- Calculate final tendency
1980	tend(k,j,i) = tend(k,j,i) + pre_tend
1981
1982	ENDDO
1983	ENDDO
1984	ENDDO
1985
1986	!
1987	!-- v-component
1988	CASE ( 2 )
1989	DO i = nxl, nxr
1990	DO j = nysv, nyn
1991	!
1992	!-- Set control flags indicating north- and southward-orientated building edges.
1993	!-- Note, building_egde_s is set from the perspective of the potential rooftop grid
1994	!-- point, while building_edge_n is set from the perspective of the non-building grid
1995	!-- point.
1996	building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) &
1997	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
1998	building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) ) &
1999	.AND. .NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
2000	!
2001	!-- Determine topography-top index on v-grid
2002	DO k = topo_top_ind(j,i,2)+1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
2003
2004	kk = k - topo_top_ind(j,i,2) !- lad arrays are defined flat
2005	!
2006	!-- In order to create sharp boundaries of the plant canopy, the lad on the v-grid
2007	!-- at index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from
2008	!-- the surrounding lad_s, because this would yield smaller lad at the canopy
2009	!-- boundaries than inside of the canopy.
2010	!-- For the same reason, the lad at the northmost (j+1) canopy boundary on the
2011	!-- v-grid equals lad_s(k,j,i), which is considered in the next if-statement.
2012	!-- Note, at left-sided building edges this is not applied, here the LAD equals
2013	!-- the LAD at grid point (k,j,i), in order to avoid that LAD is mistakenly mapped
2014	!-- on top of a roof where (usually) no LAD is defined.
2015	!-- The same is also valid for bad_s.
2016	lad_local = lad_s(kk,j,i)
2017	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
2018	.AND. .NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i)
2019
2020	bad_local = bad_s(kk,j,i)
2021	IF ( bad_local == 0.0_wp .AND. bad_s(kk,j-1,i) > 0.0_wp &
2022	.AND. .NOT. building_edge_s ) bad_local = bad_s(kk,j-1,i)
2023	!
2024	!-- In order to avoid that LAD is mistakenly considered at right-sided building
2025	!-- edges (here the topography-top index for the u-component at index j,i is still
2026	!-- on the building while the topography top for the scalar isn't), LAD is taken
2027	!-- from grid point (j,i-1). The same is also valid for bad_s.
2028	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp &
2029	.AND. building_edge_n ) lad_local = lad_s(kk,j-1,i)
2030
2031	IF ( bad_local > 0.0_wp .AND. bad_s(kk,j-1,i) == 0.0_wp &
2032	.AND. building_edge_n ) bad_local = bad_s(kk,j-1,i)
2033
2034	pre_tend = 0.0_wp
2035	pre_v = 0.0_wp
2036	!
2037	!-- Calculate preliminary value (pre_tend) of the tendency
2038	pre_tend = - canopy_drag_coeff * &
2039	( lad_local + bad_local ) * &
2040	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
2041	u(k,j-1,i+1) + &
2042	u(k,j,i) + &
2043	u(k,j,i+1) ) &
2044	)**2 + &
2045	v(k,j,i)**2 + &
2046	( 0.25_wp * ( w(k-1,j-1,i) + &
2047	w(k-1,j,i) + &
2048	w(k,j-1,i) + &
2049	w(k,j,i) ) &
2050	)**2 &
2051	) * &
2052	v(k,j,i)
2053
2054	!
2055	!-- Calculate preliminary new velocity, based on pre_tend
2056	pre_v = v(k,j,i) + dt_3d * pre_tend
2057	!
2058	!-- Compare sign of old velocity and new preliminary velocity, and in case the
2059	!-- signs are different, limit the tendency.
2060	IF ( SIGN( pre_v,v(k,j,i) ) /= pre_v ) THEN
2061	pre_tend = - v(k,j,i) * ddt_3d
2062	ELSE
2063	pre_tend = pre_tend
2064	ENDIF
2065	!
2066	!-- Calculate final tendency
2067	tend(k,j,i) = tend(k,j,i) + pre_tend
2068
2069	ENDDO
2070	ENDDO
2071	ENDDO
2072
2073	!
2074	!-- w-component
2075	CASE ( 3 )
2076	DO i = nxl, nxr
2077	DO j = nys, nyn
2078	!
2079	!-- Determine topography-top index on w-grid
2080	DO k = topo_top_ind(j,i,3)+1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
2081
2082	kk = k - topo_top_ind(j,i,3) !- lad arrays are defined flat
2083
2084	pre_tend = 0.0_wp
2085	pre_w = 0.0_wp
2086	!
2087	!-- Calculate preliminary value (pre_tend) of the tendency
2088	pre_tend = - canopy_drag_coeff * &
2089	( 0.5_wp * ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) + &
2090	0.5_wp * ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) * &
2091	SQRT( ( 0.25_wp * ( u(k,j,i) + &
2092	u(k,j,i+1) + &
2093	u(k+1,j,i) + &
2094	u(k+1,j,i+1) ) &
2095	)**2 + &
2096	( 0.25_wp * ( v(k,j,i) + &
2097	v(k,j+1,i) + &
2098	v(k+1,j,i) + &
2099	v(k+1,j+1,i) ) &
2100	)**2 + &
2101	w(k,j,i)**2 &
2102	) * &
2103	w(k,j,i)
2104	!
2105	!-- Calculate preliminary new velocity, based on pre_tend
2106	pre_w = w(k,j,i) + dt_3d * pre_tend
2107	!
2108	!-- Compare sign of old velocity and new preliminary velocity, and in case the
2109	!-- signs are different, limit the tendency
2110	IF ( SIGN( pre_w,w(k,j,i) ) /= pre_w ) THEN
2111	pre_tend = - w(k,j,i) * ddt_3d
2112	ELSE
2113	pre_tend = pre_tend
2114	ENDIF
2115	!
2116	!-- Calculate final tendency
2117	tend(k,j,i) = tend(k,j,i) + pre_tend
2118
2119	ENDDO
2120	ENDDO
2121	ENDDO
2122
2123	!
2124	!-- potential temperature
2125	CASE ( 4 )
2126	IF ( humidity ) THEN
2127	DO i = nxl, nxr
2128	DO j = nys, nyn
2129	!-- Determine topography-top index on scalar-grid
2130	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2131	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2132	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
2133	ENDDO
2134	ENDDO
2135	ENDDO
2136	ELSE
2137	DO i = nxl, nxr
2138	DO j = nys, nyn
2139	!-- Determine topography-top index on scalar-grid
2140	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2141	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2142	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
2143	ENDDO
2144	ENDDO
2145	ENDDO
2146	ENDIF
2147
2148	!
2149	!-- humidity
2150	CASE ( 5 )
2151	DO i = nxl, nxr
2152	DO j = nys, nyn
2153	!
2154	!-- Determine topography-top index on scalar-grid
2155	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2156
2157	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2158
2159	IF ( .NOT. plant_canopy_transpiration ) THEN
2160	! pcm_transpiration_rate is calculated in radiation model
2161	! in case of plant_canopy_transpiration = .T.
2162	! to include also the dependecy to the radiation
2163	! in the plant canopy box
2164	pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff &
2165	* lad_s(kk,j,i) * &
2166	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2167	u(k,j,i+1) ) &
2168	)**2 + &
2169	( 0.5_wp * ( v(k,j,i) + &
2170	v(k,j+1,i) ) &
2171	)**2 + &
2172	( 0.5_wp * ( w(k-1,j,i) + &
2173	w(k,j,i) ) &
2174	)**2 &
2175	) * &
2176	( q(k,j,i) - leaf_surface_conc )
2177	ENDIF
2178
2179	tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
2180	ENDDO
2181	ENDDO
2182	ENDDO
2183
2184	!
2185	!-- sgs-tke
2186	CASE ( 6 )
2187	DO i = nxl, nxr
2188	DO j = nys, nyn
2189	!
2190	!-- Determine topography-top index on scalar-grid
2191	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2192
2193	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2194	tend(k,j,i) = tend(k,j,i) - 2.0_wp * canopy_drag_coeff * &
2195	( lad_s(kk,j,i) + bad_s(kk,j,i) ) * &
2196	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2197	u(k,j,i+1) ) &
2198	)**2 + &
2199	( 0.5_wp * ( v(k,j,i) + &
2200	v(k,j+1,i) ) &
2201	)**2 + &
2202	( 0.5_wp * ( w(k,j,i) + &
2203	w(k+1,j,i) ) &
2204	)**2 &
2205	) * &
2206	e(k,j,i)
2207	ENDDO
2208	ENDDO
2209	ENDDO
2210	!
2211	!-- scalar concentration
2212	CASE ( 7 )
2213	DO i = nxl, nxr
2214	DO j = nys, nyn
2215	!
2216	!-- Determine topography-top index on scalar-grid
2217	DO k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2218
2219	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2220	tend(k,j,i) = tend(k,j,i) - leaf_scalar_exch_coeff * &
2221	lad_s(kk,j,i) * &
2222	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2223	u(k,j,i+1) ) &
2224	)**2 + &
2225	( 0.5_wp * ( v(k,j,i) + &
2226	v(k,j+1,i) ) &
2227	)**2 + &
2228	( 0.5_wp * ( w(k-1,j,i) + &
2229	w(k,j,i) ) &
2230	)**2 &
2231	) * &
2232	( s(k,j,i) - leaf_surface_conc )
2233	ENDDO
2234	ENDDO
2235	ENDDO
2236
2237
2238
2239	CASE DEFAULT
2240
2241	WRITE( message_string, * ) 'wrong component: ', component
2242	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
2243
2244	END SELECT
2245
2246	END SUBROUTINE pcm_tendency
2247
2248
2249	!--------------------------------------------------------------------------------------------------!
2250	! Description:
2251	! ------------
2252	!> Calculation of the tendency terms, accounting for the effect of the plant canopy on momentum and
2253	!> scalar quantities.
2254	!>
2255	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 (defined on scalar grid), as
2256	!> initialized in subroutine pcm_init.
2257	!> The lad on the w-grid is vertically interpolated from the surrounding lad_s. The upper boundary
2258	!> of the canopy is defined on the w-grid at k = pch_index. Here, the lad is zero.
2259	!>
2260	!> The canopy drag must be limited (previously accounted for by calculation of a limiting canopy
2261	!> timestep for the determination of the maximum LES timestep in subroutine timestep), since it is
2262	!> physically impossible that the canopy drag alone can locally change the sign of a velocity
2263	!> component. This limitation is realized by calculating preliminary tendencies and velocities. It
2264	!> is subsequently checked if the preliminary new velocity has a different sign than the current
2265	!> velocity. If so, the tendency is limited in a way that the velocity can at maximum be reduced to
2266	!> zero by the canopy drag.
2267	!>
2268	!>
2269	!> Call for grid point i,j
2270	!--------------------------------------------------------------------------------------------------!
2271	SUBROUTINE pcm_tendency_ij( i, j, component )
2272
2273	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
2274	INTEGER(iwp) :: i !< running index
2275	INTEGER(iwp) :: j !< running index
2276	INTEGER(iwp) :: k !< running index
2277	INTEGER(iwp) :: kk !< running index for flat lad arrays
2278
2279	LOGICAL :: building_edge_e !< control flag indicating an eastward-facing building edge
2280	LOGICAL :: building_edge_n !< control flag indicating a north-facing building edge
2281	LOGICAL :: building_edge_s !< control flag indicating a south-facing building edge
2282	LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge
2283
2284	REAL(wp) :: bad_local !< local lad value
2285	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
2286	REAL(wp) :: lad_local !< local lad value
2287	REAL(wp) :: pre_tend !< preliminary tendency
2288	REAL(wp) :: pre_u !< preliminary u-value
2289	REAL(wp) :: pre_v !< preliminary v-value
2290	REAL(wp) :: pre_w !< preliminary w-value
2291
2292
2293	ddt_3d = 1.0_wp / dt_3d
2294	!
2295	!-- Compute drag for the three velocity components and the SGS-TKE
2296	SELECT CASE ( component )
2297
2298	!
2299	!-- u-component
2300	CASE ( 1 )
2301	!
2302	!-- Set control flags indicating east- and westward-orientated building edges. Note,
2303	!-- building_egde_w is set from the perspective of the potential rooftop grid point, while
2304	!-- building_edge_e is set from the perspective of the non-building grid point.
2305	building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
2306	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
2307	building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) ) .AND. &
2308	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
2309	!
2310	!-- Determine topography-top index on u-grid
2311	DO k = topo_top_ind(j,i,1) + 1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
2312
2313	kk = k - topo_top_ind(j,i,1) !- lad arrays are defined flat
2314
2315	!
2316	!-- In order to create sharp boundaries of the plant canopy, the lad on the u-grid at
2317	!-- index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from the
2318	!-- surrounding lad_s, because this would yield smaller lad at the canopy boundaries
2319	!-- than inside of the canopy.
2320	!-- For the same reason, the lad at the rightmost(i+1)canopy boundary on the u-grid
2321	!-- equals lad_s(k,j,i), which is considered in the next if-statement. Note, at
2322	!-- left-sided building edges this is not applied, here the LAD is equals the LAD at
2323	!-- grid point (k,j,i), in order to avoid that LAD is mistakenly mapped on top of a roof
2324	!-- where (usually) is no LAD is defined.
2325	!-- The same is also valid for bad_s.
2326	lad_local = lad_s(kk,j,i)
2327	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp .AND. &
2328	.NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1)
2329
2330	bad_local = bad_s(kk,j,i)
2331	IF ( bad_local == 0.0_wp .AND. bad_s(kk,j,i-1) > 0.0_wp .AND. &
2332	.NOT. building_edge_w ) bad_local = bad_s(kk,j,i-1)
2333	!
2334	!-- In order to avoid that LAD is mistakenly considered at right-sided building edges
2335	!-- (here the topography-top index for the u-component at index j,i is still on the
2336	!-- building while the topography top for the scalar isn't), LAD is taken from grid
2337	!-- point (j,i-1). The same is also valid for bad_s.
2338	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp .AND. &
2339	building_edge_e ) lad_local = lad_s(kk,j,i-1)
2340	IF ( bad_local > 0.0_wp .AND. bad_s(kk,j,i-1) == 0.0_wp .AND. &
2341	building_edge_e ) bad_local = bad_s(kk,j,i-1)
2342
2343	pre_tend = 0.0_wp
2344	pre_u = 0.0_wp
2345	!
2346	!-- Calculate preliminary value (pre_tend) of the tendency
2347	pre_tend = - canopy_drag_coeff * &
2348	( lad_local + bad_local ) * &
2349	SQRT( u(k,j,i)**2 + &
2350	( 0.25_wp * ( v(k,j,i-1) + &
2351	v(k,j,i) + &
2352	v(k,j+1,i) + &
2353	v(k,j+1,i-1) ) &
2354	)**2 + &
2355	( 0.25_wp * ( w(k-1,j,i-1) + &
2356	w(k-1,j,i) + &
2357	w(k,j,i-1) + &
2358	w(k,j,i) ) &
2359	)**2 &
2360	) * &
2361	u(k,j,i)
2362
2363	!
2364	!-- Calculate preliminary new velocity, based on pre_tend
2365	pre_u = u(k,j,i) + dt_3d * pre_tend
2366	!
2367	!-- Compare sign of old velocity and new preliminary velocity, and in case the signs are
2368	!-- different, limit the tendency.
2369	IF ( SIGN( pre_u,u(k,j,i) ) /= pre_u ) THEN
2370	pre_tend = - u(k,j,i) * ddt_3d
2371	ELSE
2372	pre_tend = pre_tend
2373	ENDIF
2374	!
2375	!-- Calculate final tendency
2376	tend(k,j,i) = tend(k,j,i) + pre_tend
2377	ENDDO
2378
2379
2380	!
2381	!-- v-component
2382	CASE ( 2 )
2383	!
2384	!-- Set control flags indicating north- and southward-orientated building edges. Note,
2385	!-- building_egde_s is set from the perspective of the potential rooftop grid point, while
2386	!-- building_edge_n is set from the perspective of the non-building grid point.
2387	building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. &
2388	.NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
2389	building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) ) .AND. &
2390	.NOT. ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )
2391	!
2392	!-- Determine topography-top index on v-grid
2393	DO k = topo_top_ind(j,i,2) + 1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
2394
2395	kk = k - topo_top_ind(j,i,2) !- lad arrays are defined flat
2396	!
2397	!-- In order to create sharp boundaries of the plant canopy, the lad on the v-grid at
2398	!-- index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from the
2399	!-- surrounding lad_s, because this would yield smaller lad at the canopy boundaries
2400	!-- than inside of the canopy.
2401	!-- For the same reason, the lad at the northmost (j+1) canopy boundary on the v-grid
2402	!-- equals lad_s(k,j,i), which is considered in the next if-statement. Note, at
2403	!-- left-sided building edges this is not applied, here the LAD is equals the LAD at
2404	!-- grid point (k,j,i), in order to avoid that LAD is mistakenly mapped on top of a roof
2405	!-- where (usually) is no LAD is defined.
2406	!-- The same is also valid for bad_s.
2407	lad_local = lad_s(kk,j,i)
2408	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp .AND. &
2409	.NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i)
2410
2411	bad_local = bad_s(kk,j,i)
2412	IF ( bad_local == 0.0_wp .AND. bad_s(kk,j-1,i) > 0.0_wp .AND. &
2413	.NOT. building_edge_s ) bad_local = bad_s(kk,j-1,i)
2414	!
2415	!-- In order to avoid that LAD is mistakenly considered at right-sided building edges
2416	!-- (here the topography-top index for the u-component at index j,i is still on the
2417	!-- building while the topography top for the scalar isn't), LAD is taken from grid
2418	!-- point (j,i-1). The same is also valid for bad_s.
2419	IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp .AND. &
2420	building_edge_n ) lad_local = lad_s(kk,j-1,i)
2421	IF ( bad_local > 0.0_wp .AND. bad_s(kk,j-1,i) == 0.0_wp .AND. &
2422	building_edge_n ) bad_local = bad_s(kk,j-1,i)
2423
2424	pre_tend = 0.0_wp
2425	pre_v = 0.0_wp
2426	!
2427	!-- Calculate preliminary value (pre_tend) of the tendency
2428	pre_tend = - canopy_drag_coeff * &
2429	( lad_local + bad_local ) * &
2430	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
2431	u(k,j-1,i+1) + &
2432	u(k,j,i) + &
2433	u(k,j,i+1) ) &
2434	)**2 + &
2435	v(k,j,i)**2 + &
2436	( 0.25_wp * ( w(k-1,j-1,i) + &
2437	w(k-1,j,i) + &
2438	w(k,j-1,i) + &
2439	w(k,j,i) ) &
2440	)**2 &
2441	) * &
2442	v(k,j,i)
2443
2444	!
2445	!-- Calculate preliminary new velocity, based on pre_tend
2446	pre_v = v(k,j,i) + dt_3d * pre_tend
2447	!
2448	!-- Compare sign of old velocity and new preliminary velocity,
2449	!-- and in case the signs are different, limit the tendency
2450	IF ( SIGN( pre_v,v(k,j,i) ) /= pre_v ) THEN
2451	pre_tend = - v(k,j,i) * ddt_3d
2452	ELSE
2453	pre_tend = pre_tend
2454	ENDIF
2455	!
2456	!-- Calculate final tendency
2457	tend(k,j,i) = tend(k,j,i) + pre_tend
2458	ENDDO
2459
2460
2461	!
2462	!-- w-component
2463	CASE ( 3 )
2464	!
2465	!-- Determine topography-top index on w-grid
2466	DO k = topo_top_ind(j,i,3) + 1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
2467
2468	kk = k - topo_top_ind(j,i,3) !- lad arrays are defined flat
2469
2470	pre_tend = 0.0_wp
2471	pre_w = 0.0_wp
2472	!
2473	!-- Calculate preliminary value (pre_tend) of the tendency
2474	pre_tend = - canopy_drag_coeff * &
2475	( 0.5_wp * ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) + &
2476	0.5_wp * ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) * &
2477	SQRT( ( 0.25_wp * ( u(k,j,i) + &
2478	u(k,j,i+1) + &
2479	u(k+1,j,i) + &
2480	u(k+1,j,i+1) ) &
2481	)**2 + &
2482	( 0.25_wp * ( v(k,j,i) + &
2483	v(k,j+1,i) + &
2484	v(k+1,j,i) + &
2485	v(k+1,j+1,i) ) &
2486	)**2 + &
2487	w(k,j,i)**2 &
2488	) * &
2489	w(k,j,i)
2490	!
2491	!-- Calculate preliminary new velocity, based on pre_tend
2492	pre_w = w(k,j,i) + dt_3d * pre_tend
2493	!
2494	!-- Compare sign of old velocity and new preliminary velocity, and in case the signs are
2495	!-- different, limit the tendency.
2496	IF ( SIGN( pre_w,w(k,j,i) ) /= pre_w ) THEN
2497	pre_tend = - w(k,j,i) * ddt_3d
2498	ELSE
2499	pre_tend = pre_tend
2500	ENDIF
2501	!
2502	!-- Calculate final tendency
2503	tend(k,j,i) = tend(k,j,i) + pre_tend
2504	ENDDO
2505
2506	!
2507	!-- potential temperature
2508	CASE ( 4 )
2509	!
2510	!-- Determine topography-top index on scalar grid
2511	IF ( humidity ) THEN
2512	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2513	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2514	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
2515	ENDDO
2516	ELSE
2517	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2518	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2519	tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
2520	ENDDO
2521	ENDIF
2522
2523	!
2524	!-- humidity
2525	CASE ( 5 )
2526	!
2527	!-- Determine topography-top index on scalar grid
2528	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2529	kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat
2530	IF ( .NOT. plant_canopy_transpiration ) THEN
2531	! pcm_transpiration_rate is calculated in radiation model
2532	! in case of plant_canopy_transpiration = .T.
2533	! to include also the dependecy to the radiation
2534	! in the plant canopy box
2535	pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff &
2536	* lad_s(kk,j,i) * &
2537	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2538	u(k,j,i+1) ) &
2539	)**2 + &
2540	( 0.5_wp * ( v(k,j,i) + &
2541	v(k,j+1,i) ) &
2542	)**2 + &
2543	( 0.5_wp * ( w(k-1,j,i) + &
2544	w(k,j,i) ) &
2545	)**2 &
2546	) * &
2547	( q(k,j,i) - leaf_surface_conc )
2548	ENDIF
2549
2550	tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
2551
2552	ENDDO
2553
2554	!
2555	!-- sgs-tke
2556	CASE ( 6 )
2557	!
2558	!-- Determine topography-top index on scalar grid
2559	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2560
2561	kk = k - topo_top_ind(j,i,0)
2562	tend(k,j,i) = tend(k,j,i) - 2.0_wp * canopy_drag_coeff * &
2563	( lad_s(kk,j,i) + bad_s(kk,j,i) ) * &
2564	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2565	u(k,j,i+1) ) &
2566	)**2 + &
2567	( 0.5_wp * ( v(k,j,i) + &
2568	v(k,j+1,i) ) &
2569	)**2 + &
2570	( 0.5_wp * ( w(k,j,i) + &
2571	w(k+1,j,i) ) &
2572	)**2 &
2573	) * &
2574	e(k,j,i)
2575	ENDDO
2576	!
2577	!-- scalar concentration
2578	CASE ( 7 )
2579	!
2580	!-- Determine topography-top index on scalar grid
2581	DO k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2582
2583	kk = k - topo_top_ind(j,i,0)
2584	tend(k,j,i) = tend(k,j,i) - leaf_scalar_exch_coeff * &
2585	lad_s(kk,j,i) * &
2586	SQRT( ( 0.5_wp * ( u(k,j,i) + &
2587	u(k,j,i+1) ) &
2588	)**2 + &
2589	( 0.5_wp * ( v(k,j,i) + &
2590	v(k,j+1,i) ) &
2591	)**2 + &
2592	( 0.5_wp * ( w(k-1,j,i) + &
2593	w(k,j,i) ) &
2594	)**2 &
2595	) * &
2596	( s(k,j,i) - leaf_surface_conc )
2597	ENDDO
2598
2599	CASE DEFAULT
2600
2601	WRITE( message_string, * ) 'wrong component: ', component
2602	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
2603
2604	END SELECT
2605
2606	END SUBROUTINE pcm_tendency_ij
2607
2608	!--------------------------------------------------------------------------------------------------!
2609	! Description:
2610	! ------------
2611	!> Subroutine writes global restart data
2612	!--------------------------------------------------------------------------------------------------!
2613	SUBROUTINE pcm_wrd_global
2614
2615	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
2616
2617	CALL wrd_write_string( 'pch_index' )
2618	WRITE( 14 ) pch_index
2619
2620	ELSE IF ( restart_data_format_output(1:3) == 'mpi' ) THEN
2621
2622	CALL wrd_mpi_io( 'pch_index', pch_index )
2623
2624	ENDIF
2625
2626	END SUBROUTINE pcm_wrd_global
2627
2628	!--------------------------------------------------------------------------------------------------!
2629	! Description:
2630	! ------------
2631	!> Subroutine writes local (subdomain) restart data
2632	!--------------------------------------------------------------------------------------------------!
2633	SUBROUTINE pcm_wrd_local
2634
2635	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d !< temporary array to store pcm data with
2636	!< non-standard vertical index bounds
2637
2638
2639	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
2640
2641	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
2642	CALL wrd_write_string( 'pcm_heatrate_av' )
2643	WRITE( 14 ) pcm_heatrate_av
2644	ENDIF
2645
2646	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
2647	CALL wrd_write_string( 'pcm_latentrate_av' )
2648	WRITE( 14 ) pcm_latentrate_av
2649	ENDIF
2650
2651	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
2652	CALL wrd_write_string( 'pcm_transpirationrate_av' )
2653	WRITE( 14 ) pcm_transpirationrate_av
2654	ENDIF
2655
2656	ELSE IF ( restart_data_format_output(1:3) == 'mpi' ) THEN
2657
2658	!
2659	!-- Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
2660	!-- full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored before writing.
2661	IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
2662	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2663	tmp_3d(nzb:pch_index,:,:) = pcm_heatrate_av
2664	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2665	CALL wrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
2666	DEALLOCATE( tmp_3d )
2667	ENDIF
2668	IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
2669	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2670	tmp_3d(nzb:pch_index,:,:) = pcm_latentrate_av
2671	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2672	CALL wrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
2673	DEALLOCATE( tmp_3d )
2674	ENDIF
2675	IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
2676	ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2677	tmp_3d(nzb:pch_index,:,:) = pcm_transpirationrate_av
2678	tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2679	CALL wrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
2680	DEALLOCATE( tmp_3d )
2681	ENDIF
2682
2683	ENDIF
2684
2685	END SUBROUTINE pcm_wrd_local
2686
2687
2688	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 |