Home

Context Navigation

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

Last change on this file since 2062 was 2025, checked in by kanani, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 59.4 KB

Line
1	!> @file plant_canopy_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2016 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: plant_canopy_model_mod.f90 2025 2016-10-12 16:44:03Z maronga $
27	!
28	! 2024 2016-10-12 16:42:37Z kanani
29	! Added missing lad_s initialization
30	!
31	! 2011 2016-09-19 17:29:57Z kanani
32	! Renamed canopy_heat_flux to pc_heating_rate, since the original meaning/
33	! calculation of the quantity has changed, related to the urban surface model
34	! and similar future applications.
35	!
36	! 2007 2016-08-24 15:47:17Z kanani
37	! Added SUBROUTINE pcm_read_plant_canopy_3d for reading 3d plant canopy data
38	! from file (new case canopy_mode=read_from_file_3d) in the course of
39	! introduction of urban surface model,
40	! introduced variable ext_coef,
41	! resorted SUBROUTINEs to alphabetical order
42	!
43	! 2000 2016-08-20 18:09:15Z knoop
44	! Forced header and separation lines into 80 columns
45	!
46	! 1960 2016-07-12 16:34:24Z suehring
47	! Separate humidity and passive scalar
48	!
49	! 1953 2016-06-21 09:28:42Z suehring
50	! Bugfix, lad_s and lad must be public
51	!
52	! 1826 2016-04-07 12:01:39Z maronga
53	! Further modularization
54	!
55	! 1721 2015-11-16 12:56:48Z raasch
56	! bugfixes: shf is reduced in areas covered with canopy only,
57	! canopy is set on top of topography
58	!
59	! 1682 2015-10-07 23:56:08Z knoop
60	! Code annotations made doxygen readable
61	!
62	! 1484 2014-10-21 10:53:05Z kanani
63	! Changes due to new module structure of the plant canopy model:
64	! module plant_canopy_model_mod now contains a subroutine for the
65	! initialization of the canopy model (pcm_init),
66	! limitation of the canopy drag (previously accounted for by calculation of
67	! a limiting canopy timestep for the determination of the maximum LES timestep
68	! in subroutine timestep) is now realized by the calculation of pre-tendencies
69	! and preliminary velocities in subroutine pcm_tendency,
70	! some redundant MPI communication removed in subroutine pcm_init
71	! (was previously in init_3d_model),
72	! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
73	! respective grid is now provided only by lad_s (e.g. in the calculation of
74	! the tendency terms or of cum_lai_hf),
75	! drag_coefficient, lai, leaf_surface_concentration,
76	! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
77	! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
78	! respectively,
79	! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
80	! USE-statements and ONLY-lists modified accordingly
81	!
82	! 1340 2014-03-25 19:45:13Z kanani
83	! REAL constants defined as wp-kind
84	!
85	! 1320 2014-03-20 08:40:49Z raasch
86	! ONLY-attribute added to USE-statements,
87	! kind-parameters added to all INTEGER and REAL declaration statements,
88	! kinds are defined in new module kinds,
89	! old module precision_kind is removed,
90	! revision history before 2012 removed,
91	! comment fields (!:) to be used for variable explanations added to
92	! all variable declaration statements
93	!
94	! 1036 2012-10-22 13:43:42Z raasch
95	! code put under GPL (PALM 3.9)
96	!
97	! 138 2007-11-28 10:03:58Z letzel
98	! Initial revision
99	!
100	! Description:
101	! ------------
102	!> 1) Initialization of the canopy model, e.g. construction of leaf area density
103	!> profile (subroutine pcm_init).
104	!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
105	!> due to canopy elements (subroutine pcm_tendency).
106	!------------------------------------------------------------------------------!
107	MODULE plant_canopy_model_mod
108
109	USE arrays_3d, &
110	ONLY: dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw
111
112	USE indices, &
113	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
114	nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
115
116	USE kinds
117
118
119	IMPLICIT NONE
120
121
122	CHARACTER (LEN=20) :: canopy_mode = 'block' !< canopy coverage
123
124	INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index
125	INTEGER(iwp) :: &
126	lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
127
128	LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
129	LOGICAL :: plant_canopy = .FALSE. !< switch for use of canopy model
130
131	REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation
132	REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation
133	REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
134	REAL(wp) :: cdc = 0.0_wp !< canopy drag coeff. (abbreviation used in equations)
135	REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux
136	REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag
137	REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient
138	REAL(wp) :: lad_surface = 0.0_wp !< lad surface value
139	REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc.
140	REAL(wp) :: &
141	leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff.
142	REAL(wp) :: &
143	leaf_surface_conc = 0.0_wp !< leaf surface concentration
144	REAL(wp) :: lsec = 0.0_wp !< leaf scalar exchange coeff.
145	REAL(wp) :: lsc = 0.0_wp !< leaf surface concentration
146
147	REAL(wp) :: &
148	lad_vertical_gradient(10) = 0.0_wp !< lad gradient
149	REAL(wp) :: &
150	lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
151
152	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density
153	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad
154
155	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
156	pc_heating_rate !< plant canopy heating rate
157	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.
158	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid
159
160
161	SAVE
162
163
164	PRIVATE
165
166	!
167	!-- Public functions
168	PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
169
170	!
171	!-- Public variables and constants
172	PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, &
173	pch_index, plant_canopy
174
175
176
177	INTERFACE pcm_check_parameters
178	MODULE PROCEDURE pcm_check_parameters
179	END INTERFACE pcm_check_parameters
180
181	INTERFACE pcm_header
182	MODULE PROCEDURE pcm_header
183	END INTERFACE pcm_header
184
185	INTERFACE pcm_init
186	MODULE PROCEDURE pcm_init
187	END INTERFACE pcm_init
188
189	INTERFACE pcm_parin
190	MODULE PROCEDURE pcm_parin
191	END INTERFACE pcm_parin
192
193	INTERFACE pcm_read_plant_canopy_3d
194	MODULE PROCEDURE pcm_read_plant_canopy_3d
195	END INTERFACE pcm_read_plant_canopy_3d
196
197	INTERFACE pcm_tendency
198	MODULE PROCEDURE pcm_tendency
199	MODULE PROCEDURE pcm_tendency_ij
200	END INTERFACE pcm_tendency
201
202
203	CONTAINS
204
205
206	!------------------------------------------------------------------------------!
207	! Description:
208	! ------------
209	!> Check parameters routine for plant canopy model
210	!------------------------------------------------------------------------------!
211	SUBROUTINE pcm_check_parameters
212
213	USE control_parameters, &
214	ONLY: cloud_physics, message_string, microphysics_seifert
215
216
217	IMPLICIT NONE
218
219
220	IF ( canopy_drag_coeff == 0.0_wp ) THEN
221	message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
222	'coefficient & given value is canopy_drag_coeff = 0.0'
223	CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
224	ENDIF
225
226	IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp ) .OR.&
227	beta_lad /= 9999999.9_wp .AND. alpha_lad == 9999999.9_wp ) THEN
228	message_string = 'using the beta function for the construction ' // &
229	'of the leaf area density profile requires ' // &
230	'both alpha_lad and beta_lad to be /= 9999999.9'
231	CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
232	ENDIF
233
234	IF ( calc_beta_lad_profile .AND. lai_beta == 0.0_wp ) THEN
235	message_string = 'using the beta function for the construction ' // &
236	'of the leaf area density profile requires ' // &
237	'a non-zero lai_beta, but given value is ' // &
238	'lai_beta = 0.0'
239	CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
240	ENDIF
241
242	IF ( calc_beta_lad_profile .AND. lad_surface /= 0.0_wp ) THEN
243	message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
244	'and lad_surface /= 0.0 is not possible, ' // &
245	'use either vertical gradients or the beta ' // &
246	'function for the construction of the leaf area '// &
247	'density profile'
248	CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
249	ENDIF
250
251	IF ( cloud_physics .AND. microphysics_seifert ) THEN
252	message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
253	' seifert_beheng'
254	CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
255	ENDIF
256
257
258	END SUBROUTINE pcm_check_parameters
259
260
261	!------------------------------------------------------------------------------!
262	! Description:
263	! ------------
264	!> Header output for plant canopy model
265	!------------------------------------------------------------------------------!
266	SUBROUTINE pcm_header ( io )
267
268	USE control_parameters, &
269	ONLY: dz, passive_scalar
270
271
272	IMPLICIT NONE
273
274	CHARACTER (LEN=10) :: coor_chr !<
275
276	CHARACTER (LEN=86) :: coordinates !<
277	CHARACTER (LEN=86) :: gradients !<
278	CHARACTER (LEN=86) :: leaf_area_density !<
279	CHARACTER (LEN=86) :: slices !<
280
281	INTEGER(iwp) :: i !<
282	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
283	INTEGER(iwp) :: k !<
284
285	REAL(wp) :: canopy_height !< canopy height (in m)
286
287	canopy_height = pch_index * dz
288
289	WRITE ( io, 1 ) canopy_mode, canopy_height, pch_index, &
290	canopy_drag_coeff
291	IF ( passive_scalar ) THEN
292	WRITE ( io, 2 ) leaf_scalar_exch_coeff, &
293	leaf_surface_conc
294	ENDIF
295
296	!
297	!-- Heat flux at the top of vegetation
298	WRITE ( io, 3 ) cthf
299
300	!
301	!-- Leaf area density profile, calculated either from given vertical
302	!-- gradients or from beta probability density function.
303	IF ( .NOT. calc_beta_lad_profile ) THEN
304
305	!-- Building output strings, starting with surface value
306	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
307	gradients = '------'
308	slices = ' 0'
309	coordinates = ' 0.0'
310	i = 1
311	DO WHILE ( i < 11 .AND. lad_vertical_gradient_level_ind(i) &
312	/= -9999 )
313
314	WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
315	leaf_area_density = TRIM( leaf_area_density ) // ' ' // &
316	TRIM( coor_chr )
317
318	WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
319	gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
320
321	WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
322	slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
323
324	WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
325	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
326
327	i = i + 1
328	ENDDO
329
330	WRITE ( io, 4 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
331	TRIM( gradients ), TRIM( slices )
332
333	ELSE
334
335	WRITE ( leaf_area_density, '(F7.4)' ) lad_surface
336	coordinates = ' 0.0'
337
338	DO k = 1, pch_index
339
340	WRITE (coor_chr,'(F7.2)') lad(k)
341	leaf_area_density = TRIM( leaf_area_density ) // ' ' // &
342	TRIM( coor_chr )
343
344	WRITE (coor_chr,'(F7.1)') zu(k)
345	coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr )
346
347	ENDDO
348
349	WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ), &
350	alpha_lad, beta_lad, lai_beta
351
352	ENDIF
353
354	1 FORMAT (//' Vegetation canopy (drag) model:'/ &
355	' ------------------------------'// &
356	' Canopy mode: ', A / &
357	' Canopy height: ',F6.2,'m (',I4,' grid points)' / &
358	' Leaf drag coefficient: ',F6.2 /)
359	2 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / &
360	' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
361	3 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
362	' K m/s')
363	4 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
364	' Height: ',A,' m'/ &
365	' Leaf area density: ',A,' m2/m3'/ &
366	' Gradient: ',A,' m2/m4'/ &
367	' Gridpoint: ',A)
368	5 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
369	// ' Height: ',A,' m'/ &
370	' Leaf area density: ',A,' m2/m3'/ &
371	' Coefficient alpha: ',F6.2 / &
372	' Coefficient beta: ',F6.2 / &
373	' Leaf area index: ',F6.2,' m2/m2' /)
374
375	END SUBROUTINE pcm_header
376
377
378	!------------------------------------------------------------------------------!
379	! Description:
380	! ------------
381	!> Initialization of the plant canopy model
382	!------------------------------------------------------------------------------!
383	SUBROUTINE pcm_init
384
385
386	USE control_parameters, &
387	ONLY: coupling_char, dz, humidity, io_blocks, io_group, &
388	message_string, ocean, passive_scalar
389
390
391	IMPLICIT NONE
392
393	CHARACTER(10) :: pct
394
395	INTEGER(iwp) :: i !< running index
396	INTEGER(iwp) :: ii !< index
397	INTEGER(iwp) :: j !< running index
398	INTEGER(iwp) :: k !< running index
399
400	REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction
401	REAL(wp) :: dzh !< vertical grid spacing in units of canopy height
402	REAL(wp) :: gradient !< gradient for lad-profile construction
403	REAL(wp) :: canopy_height !< canopy height for lad-profile construction
404	REAL(wp) :: pcv(nzb:nzt+1) !<
405
406	!
407	!-- Allocate one-dimensional arrays for the computation of the
408	!-- leaf area density (lad) profile
409	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
410	lad = 0.0_wp
411	pre_lad = 0.0_wp
412
413	!
414	!-- Set flag that indicates that the lad-profile shall be calculated by using
415	!-- a beta probability density function
416	IF ( alpha_lad /= 9999999.9_wp .AND. beta_lad /= 9999999.9_wp ) THEN
417	calc_beta_lad_profile = .TRUE.
418	ENDIF
419
420
421	!
422	!-- Compute the profile of leaf area density used in the plant
423	!-- canopy model. The profile can either be constructed from
424	!-- prescribed vertical gradients of the leaf area density or by
425	!-- using a beta probability density function (see e.g. Markkanen et al.,
426	!-- 2003: Boundary-Layer Meteorology, 106, 437-459)
427	IF ( .NOT. calc_beta_lad_profile ) THEN
428
429	!
430	!-- Use vertical gradients for lad-profile construction
431	i = 1
432	gradient = 0.0_wp
433
434	IF ( .NOT. ocean ) THEN
435
436	lad(0) = lad_surface
437	lad_vertical_gradient_level_ind(1) = 0
438
439	DO k = 1, pch_index
440	IF ( i < 11 ) THEN
441	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
442	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
443	gradient = lad_vertical_gradient(i)
444	lad_vertical_gradient_level_ind(i) = k - 1
445	i = i + 1
446	ENDIF
447	ENDIF
448	IF ( gradient /= 0.0_wp ) THEN
449	IF ( k /= 1 ) THEN
450	lad(k) = lad(k-1) + dzu(k) * gradient
451	ELSE
452	lad(k) = lad_surface + dzu(k) * gradient
453	ENDIF
454	ELSE
455	lad(k) = lad(k-1)
456	ENDIF
457	ENDDO
458
459	ENDIF
460
461	!
462	!-- In case of no given leaf area density gradients, choose a vanishing
463	!-- gradient. This information is used for the HEADER and the RUN_CONTROL
464	!-- file.
465	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
466	lad_vertical_gradient_level(1) = 0.0_wp
467	ENDIF
468
469	ELSE
470
471	!
472	!-- Use beta function for lad-profile construction
473	int_bpdf = 0.0_wp
474	canopy_height = pch_index * dz
475
476	DO k = nzb, pch_index
477	int_bpdf = int_bpdf + &
478	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
479	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
480	beta_lad-1.0_wp ) ) &
481	* ( ( zw(k+1)-zw(k) ) / canopy_height ) )
482	ENDDO
483
484	!
485	!-- Preliminary lad profile (defined on w-grid)
486	DO k = nzb, pch_index
487	pre_lad(k) = lai_beta * &
488	( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) &
489	* ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( &
490	beta_lad-1.0_wp ) ) / int_bpdf &
491	) / canopy_height
492	ENDDO
493
494	!
495	!-- Final lad profile (defined on scalar-grid level, since most prognostic
496	!-- quantities are defined there, hence, less interpolation is required
497	!-- when calculating the canopy tendencies)
498	lad(0) = pre_lad(0)
499	DO k = nzb+1, pch_index
500	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
501	ENDDO
502
503	ENDIF
504
505	!
506	!-- Allocate 3D-array for the leaf area density (lad_s). In case of a
507	!-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
508	!-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
509	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
510
511	IF ( cthf /= 0.0_wp ) THEN
512	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
513	pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
514	ENDIF
515
516	!
517	!-- Initialize canopy parameters cdc (canopy drag coefficient),
518	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
519	!-- with the prescribed values
520	cdc = canopy_drag_coeff
521	lsec = leaf_scalar_exch_coeff
522	lsc = leaf_surface_conc
523
524	!
525	!-- Initialization of the canopy coverage in the model domain:
526	!-- Setting the parameter canopy_mode = 'block' initializes a canopy, which
527	!-- fully covers the domain surface
528	SELECT CASE ( TRIM( canopy_mode ) )
529
530	CASE( 'block' )
531
532	DO i = nxlg, nxrg
533	DO j = nysg, nyng
534	lad_s(:,j,i) = lad(:)
535	ENDDO
536	ENDDO
537
538	CASE ( 'read_from_file_3d' )
539	!
540	!-- Initialize canopy parameters cdc (canopy drag coefficient),
541	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
542	!-- from file which contains complete 3D data (separate vertical profiles for
543	!-- each location).
544	CALL pcm_read_plant_canopy_3d
545
546	CASE DEFAULT
547	!
548	!-- The DEFAULT case is reached either if the parameter
549	!-- canopy mode contains a wrong character string or if the
550	!-- user has coded a special case in the user interface.
551	!-- There, the subroutine user_init_plant_canopy checks
552	!-- which of these two conditions applies.
553	CALL user_init_plant_canopy
554
555	END SELECT
556
557	!
558	!-- Initialization of the canopy heat source distribution due to heating
559	!-- of the canopy layers by incoming solar radiation, in case that a non-zero
560	!-- value is set for the canopy top heat flux (cthf), which equals the
561	!-- available net radiation at canopy top.
562	!-- The heat source distribution is calculated by a decaying exponential
563	!-- function of the downward cumulative leaf area index (cum_lai_hf),
564	!-- assuming that the foliage inside the plant canopy is heated by solar
565	!-- radiation penetrating the canopy layers according to the distribution
566	!-- of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
567	!-- 73â96). This approach has been applied e.g. by Shaw & Schumann (1992;
568	!-- Bound.-Layer Meteorol. 61, 47â64)
569	IF ( cthf /= 0.0_wp ) THEN
570	!
571	!-- Piecewise calculation of the cumulative leaf area index by vertical
572	!-- integration of the leaf area density
573	cum_lai_hf(:,:,:) = 0.0_wp
574	DO i = nxlg, nxrg
575	DO j = nysg, nyng
576	DO k = pch_index-1, 0, -1
577	IF ( k == pch_index-1 ) THEN
578	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
579	( 0.5_wp * lad_s(k+1,j,i) * &
580	( zw(k+1) - zu(k+1) ) ) + &
581	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
582	lad_s(k,j,i) ) + &
583	lad_s(k+1,j,i) ) * &
584	( zu(k+1) - zw(k) ) )
585	ELSE
586	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
587	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
588	lad_s(k+1,j,i) ) + &
589	lad_s(k+1,j,i) ) * &
590	( zw(k+1) - zu(k+1) ) ) + &
591	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
592	lad_s(k,j,i) ) + &
593	lad_s(k+1,j,i) ) * &
594	( zu(k+1) - zw(k) ) )
595	ENDIF
596	ENDDO
597	ENDDO
598	ENDDO
599
600	!
601	!-- Calculation of the heating rate (K/s) within the different layers of
602	!-- the plant canopy
603	DO i = nxlg, nxrg
604	DO j = nysg, nyng
605	!
606	!-- Calculation only necessary in areas covered with canopy
607	IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) THEN
608	!--
609	!-- In areas with canopy the surface value of the canopy heat
610	!-- flux distribution overrides the surface heat flux (shf)
611	shf(j,i) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
612	!
613	!-- Within the different canopy layers the plant-canopy heating
614	!-- rate (pc_heating_rate) is calculated as the vertical
615	!-- divergence of the canopy heat fluxes at the top and bottom
616	!-- of the respective layer
617	DO k = 1, pch_index
618	pc_heating_rate(k,j,i) = cthf * &
619	( exp(-ext_coef*cum_lai_hf(k,j,i)) - &
620	exp(-ext_coef*cum_lai_hf(k-1,j,i)) ) / dzw(k)
621	ENDDO
622	ENDIF
623	ENDDO
624	ENDDO
625
626	ENDIF
627
628
629
630	END SUBROUTINE pcm_init
631
632
633	!------------------------------------------------------------------------------!
634	! Description:
635	! ------------
636	!> Parin for &canopy_par for plant canopy model
637	!------------------------------------------------------------------------------!
638	SUBROUTINE pcm_parin
639
640
641	IMPLICIT NONE
642
643	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
644
645	NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, &
646	canopy_mode, cthf, &
647	lad_surface, &
648	lad_vertical_gradient, &
649	lad_vertical_gradient_level, &
650	lai_beta, &
651	leaf_scalar_exch_coeff, &
652	leaf_surface_conc, pch_index
653
654	line = ' '
655
656	!
657	!-- Try to find radiation model package
658	REWIND ( 11 )
659	line = ' '
660	DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
661	READ ( 11, '(A)', END=10 ) line
662	ENDDO
663	BACKSPACE ( 11 )
664
665	!
666	!-- Read user-defined namelist
667	READ ( 11, canopy_par )
668
669	!
670	!-- Set flag that indicates that the radiation model is switched on
671	plant_canopy = .TRUE.
672
673	10 CONTINUE
674
675
676	END SUBROUTINE pcm_parin
677
678
679
680	!------------------------------------------------------------------------------!
681	! Description:
682	! ------------
683	!
684	!> Loads 3D plant canopy data from file. File format is as follows:
685	!>
686	!> num_levels
687	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
688	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
689	!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
690	!> ...
691	!>
692	!> i.e. first line determines number of levels and further lines represent plant
693	!> canopy data, one line per column and variable. In each data line,
694	!> dtype represents variable to be set:
695	!>
696	!> dtype=1: leaf area density (lad_s)
697	!> dtype=2: canopy drag coefficient (cdc)
698	!> dtype=3: leaf scalar exchange coefficient (lsec)
699	!> dtype=4: leaf surface concentration (lsc)
700	!>
701	!> Zeros are added automatically above num_levels until top of domain. Any
702	!> non-specified (x,y) columns have zero values as default.
703	!------------------------------------------------------------------------------!
704	SUBROUTINE pcm_read_plant_canopy_3d
705	USE control_parameters, &
706	ONLY: passive_scalar, message_string
707	IMPLICIT NONE
708
709	INTEGER(iwp) :: i, j, dtype, nzp, nzpltop, nzpl, kk
710	REAL(wp), DIMENSION(:), ALLOCATABLE :: col
711
712	lad_s = 0.0_wp
713	OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
714	action='READ', status='OLD', form='FORMATTED', err=515)
715	READ(152, *, err=516, end=517) nzp !< read first line = number of vertical layers
716	ALLOCATE(col(1:nzp))
717	nzpltop = MIN(nzt+1, nzb+nzp-1)
718	nzpl = nzpltop - nzb + 1 !< no. of layers to assign
719
720	DO
721	READ(152, *, err=516, end=517) dtype, i, j, col(:)
722	IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
723
724	SELECT CASE (dtype)
725	CASE( 1 ) !< leaf area density
726	!-- only lad_s has flat z-coordinate, others have regular
727	kk = nzb_s_inner(j, i)
728	lad_s(nzb:nzpltop-kk, j, i) = col(1+kk:nzpl)
729	! CASE( 2 ) !< canopy drag coefficient
730	! cdc(nzb:nzpltop, j, i) = col(1:nzpl)
731	! CASE( 3 ) !< leaf scalar exchange coefficient
732	! lsec(nzb:nzpltop, j, i) = col(1:nzpl)
733	! CASE( 4 ) !< leaf surface concentration
734	! lsc(nzb:nzpltop, j, i) = col(1:nzpl)
735	CASE DEFAULT
736	write(message_string, '(a,i2,a)') &
737	'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
738	CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
739	END SELECT
740	ENDDO
741
742	515 message_string = 'error opening file PLANT_CANOPY_DATA_3D'
743	CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
744
745	516 message_string = 'error reading file PLANT_CANOPY_DATA_3D'
746	CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
747
748	517 CLOSE(152)
749	DEALLOCATE(col)
750
751	END SUBROUTINE pcm_read_plant_canopy_3d
752
753
754
755	!------------------------------------------------------------------------------!
756	! Description:
757	! ------------
758	!> Calculation of the tendency terms, accounting for the effect of the plant
759	!> canopy on momentum and scalar quantities.
760	!>
761	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
762	!> (defined on scalar grid), as initialized in subroutine pcm_init.
763	!> The lad on the w-grid is vertically interpolated from the surrounding
764	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
765	!> k = pch_index. Here, the lad is zero.
766	!>
767	!> The canopy drag must be limited (previously accounted for by calculation of
768	!> a limiting canopy timestep for the determination of the maximum LES timestep
769	!> in subroutine timestep), since it is physically impossible that the canopy
770	!> drag alone can locally change the sign of a velocity component. This
771	!> limitation is realized by calculating preliminary tendencies and velocities.
772	!> It is subsequently checked if the preliminary new velocity has a different
773	!> sign than the current velocity. If so, the tendency is limited in a way that
774	!> the velocity can at maximum be reduced to zero by the canopy drag.
775	!>
776	!>
777	!> Call for all grid points
778	!------------------------------------------------------------------------------!
779	SUBROUTINE pcm_tendency( component )
780
781
782	USE control_parameters, &
783	ONLY: dt_3d, message_string
784
785	USE kinds
786
787	IMPLICIT NONE
788
789	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
790	INTEGER(iwp) :: i !< running index
791	INTEGER(iwp) :: j !< running index
792	INTEGER(iwp) :: k !< running index
793	INTEGER(iwp) :: kk !< running index for flat lad arrays
794
795	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
796	REAL(wp) :: lad_local !< local lad value
797	REAL(wp) :: pre_tend !< preliminary tendency
798	REAL(wp) :: pre_u !< preliminary u-value
799	REAL(wp) :: pre_v !< preliminary v-value
800	REAL(wp) :: pre_w !< preliminary w-value
801
802
803	ddt_3d = 1.0_wp / dt_3d
804
805	!
806	!-- Compute drag for the three velocity components and the SGS-TKE:
807	SELECT CASE ( component )
808
809	!
810	!-- u-component
811	CASE ( 1 )
812	DO i = nxlu, nxr
813	DO j = nys, nyn
814	DO k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
815
816	kk = k - nzb_u_inner(j,i) !- lad arrays are defined flat
817	!
818	!-- In order to create sharp boundaries of the plant canopy,
819	!-- the lad on the u-grid at index (k,j,i) is equal to
820	!-- lad_s(k,j,i), rather than being interpolated from the
821	!-- surrounding lad_s, because this would yield smaller lad
822	!-- at the canopy boundaries than inside of the canopy.
823	!-- For the same reason, the lad at the rightmost(i+1)canopy
824	!-- boundary on the u-grid equals lad_s(k,j,i).
825	lad_local = lad_s(kk,j,i)
826	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
827	THEN
828	lad_local = lad_s(kk,j,i-1)
829	ENDIF
830
831	pre_tend = 0.0_wp
832	pre_u = 0.0_wp
833	!
834	!-- Calculate preliminary value (pre_tend) of the tendency
835	pre_tend = - cdc * &
836	lad_local * &
837	SQRT( u(k,j,i)**2 + &
838	( 0.25_wp * ( v(k,j,i-1) + &
839	v(k,j,i) + &
840	v(k,j+1,i) + &
841	v(k,j+1,i-1) ) &
842	)**2 + &
843	( 0.25_wp * ( w(k-1,j,i-1) + &
844	w(k-1,j,i) + &
845	w(k,j,i-1) + &
846	w(k,j,i) ) &
847	)**2 &
848	) * &
849	u(k,j,i)
850
851	!
852	!-- Calculate preliminary new velocity, based on pre_tend
853	pre_u = u(k,j,i) + dt_3d * pre_tend
854	!
855	!-- Compare sign of old velocity and new preliminary velocity,
856	!-- and in case the signs are different, limit the tendency
857	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
858	pre_tend = - u(k,j,i) * ddt_3d
859	ELSE
860	pre_tend = pre_tend
861	ENDIF
862	!
863	!-- Calculate final tendency
864	tend(k,j,i) = tend(k,j,i) + pre_tend
865
866	ENDDO
867	ENDDO
868	ENDDO
869
870	!
871	!-- v-component
872	CASE ( 2 )
873	DO i = nxl, nxr
874	DO j = nysv, nyn
875	DO k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
876
877	kk = k - nzb_v_inner(j,i) !- lad arrays are defined flat
878	!
879	!-- In order to create sharp boundaries of the plant canopy,
880	!-- the lad on the v-grid at index (k,j,i) is equal to
881	!-- lad_s(k,j,i), rather than being interpolated from the
882	!-- surrounding lad_s, because this would yield smaller lad
883	!-- at the canopy boundaries than inside of the canopy.
884	!-- For the same reason, the lad at the northmost(j+1) canopy
885	!-- boundary on the v-grid equals lad_s(k,j,i).
886	lad_local = lad_s(kk,j,i)
887	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
888	THEN
889	lad_local = lad_s(kk,j-1,i)
890	ENDIF
891
892	pre_tend = 0.0_wp
893	pre_v = 0.0_wp
894	!
895	!-- Calculate preliminary value (pre_tend) of the tendency
896	pre_tend = - cdc * &
897	lad_local * &
898	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
899	u(k,j-1,i+1) + &
900	u(k,j,i) + &
901	u(k,j,i+1) ) &
902	)**2 + &
903	v(k,j,i)**2 + &
904	( 0.25_wp * ( w(k-1,j-1,i) + &
905	w(k-1,j,i) + &
906	w(k,j-1,i) + &
907	w(k,j,i) ) &
908	)**2 &
909	) * &
910	v(k,j,i)
911
912	!
913	!-- Calculate preliminary new velocity, based on pre_tend
914	pre_v = v(k,j,i) + dt_3d * pre_tend
915	!
916	!-- Compare sign of old velocity and new preliminary velocity,
917	!-- and in case the signs are different, limit the tendency
918	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
919	pre_tend = - v(k,j,i) * ddt_3d
920	ELSE
921	pre_tend = pre_tend
922	ENDIF
923	!
924	!-- Calculate final tendency
925	tend(k,j,i) = tend(k,j,i) + pre_tend
926
927	ENDDO
928	ENDDO
929	ENDDO
930
931	!
932	!-- w-component
933	CASE ( 3 )
934	DO i = nxl, nxr
935	DO j = nys, nyn
936	DO k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
937
938	kk = k - nzb_w_inner(j,i) !- lad arrays are defined flat
939
940	pre_tend = 0.0_wp
941	pre_w = 0.0_wp
942	!
943	!-- Calculate preliminary value (pre_tend) of the tendency
944	pre_tend = - cdc * &
945	(0.5_wp * &
946	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
947	SQRT( ( 0.25_wp * ( u(k,j,i) + &
948	u(k,j,i+1) + &
949	u(k+1,j,i) + &
950	u(k+1,j,i+1) ) &
951	)**2 + &
952	( 0.25_wp * ( v(k,j,i) + &
953	v(k,j+1,i) + &
954	v(k+1,j,i) + &
955	v(k+1,j+1,i) ) &
956	)**2 + &
957	w(k,j,i)**2 &
958	) * &
959	w(k,j,i)
960	!
961	!-- Calculate preliminary new velocity, based on pre_tend
962	pre_w = w(k,j,i) + dt_3d * pre_tend
963	!
964	!-- Compare sign of old velocity and new preliminary velocity,
965	!-- and in case the signs are different, limit the tendency
966	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
967	pre_tend = - w(k,j,i) * ddt_3d
968	ELSE
969	pre_tend = pre_tend
970	ENDIF
971	!
972	!-- Calculate final tendency
973	tend(k,j,i) = tend(k,j,i) + pre_tend
974
975	ENDDO
976	ENDDO
977	ENDDO
978
979	!
980	!-- potential temperature
981	CASE ( 4 )
982	DO i = nxl, nxr
983	DO j = nys, nyn
984	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
985	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
986	tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
987	ENDDO
988	ENDDO
989	ENDDO
990
991	!
992	!-- humidity
993	CASE ( 5 )
994	DO i = nxl, nxr
995	DO j = nys, nyn
996	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
997	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
998	tend(k,j,i) = tend(k,j,i) - &
999	lsec * &
1000	lad_s(kk,j,i) * &
1001	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1002	u(k,j,i+1) ) &
1003	)**2 + &
1004	( 0.5_wp * ( v(k,j,i) + &
1005	v(k,j+1,i) ) &
1006	)**2 + &
1007	( 0.5_wp * ( w(k-1,j,i) + &
1008	w(k,j,i) ) &
1009	)**2 &
1010	) * &
1011	( q(k,j,i) - lsc )
1012	ENDDO
1013	ENDDO
1014	ENDDO
1015
1016	!
1017	!-- sgs-tke
1018	CASE ( 6 )
1019	DO i = nxl, nxr
1020	DO j = nys, nyn
1021	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1022	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1023	tend(k,j,i) = tend(k,j,i) - &
1024	2.0_wp * cdc * &
1025	lad_s(kk,j,i) * &
1026	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1027	u(k,j,i+1) ) &
1028	)**2 + &
1029	( 0.5_wp * ( v(k,j,i) + &
1030	v(k,j+1,i) ) &
1031	)**2 + &
1032	( 0.5_wp * ( w(k,j,i) + &
1033	w(k+1,j,i) ) &
1034	)**2 &
1035	) * &
1036	e(k,j,i)
1037	ENDDO
1038	ENDDO
1039	ENDDO
1040	!
1041	!-- scalar concentration
1042	CASE ( 7 )
1043	DO i = nxl, nxr
1044	DO j = nys, nyn
1045	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1046	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1047	tend(k,j,i) = tend(k,j,i) - &
1048	lsec * &
1049	lad_s(kk,j,i) * &
1050	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1051	u(k,j,i+1) ) &
1052	)**2 + &
1053	( 0.5_wp * ( v(k,j,i) + &
1054	v(k,j+1,i) ) &
1055	)**2 + &
1056	( 0.5_wp * ( w(k-1,j,i) + &
1057	w(k,j,i) ) &
1058	)**2 &
1059	) * &
1060	( s(k,j,i) - lsc )
1061	ENDDO
1062	ENDDO
1063	ENDDO
1064
1065
1066
1067	CASE DEFAULT
1068
1069	WRITE( message_string, * ) 'wrong component: ', component
1070	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
1071
1072	END SELECT
1073
1074	END SUBROUTINE pcm_tendency
1075
1076
1077	!------------------------------------------------------------------------------!
1078	! Description:
1079	! ------------
1080	!> Calculation of the tendency terms, accounting for the effect of the plant
1081	!> canopy on momentum and scalar quantities.
1082	!>
1083	!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1084	!> (defined on scalar grid), as initialized in subroutine pcm_init.
1085	!> The lad on the w-grid is vertically interpolated from the surrounding
1086	!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1087	!> k = pch_index. Here, the lad is zero.
1088	!>
1089	!> The canopy drag must be limited (previously accounted for by calculation of
1090	!> a limiting canopy timestep for the determination of the maximum LES timestep
1091	!> in subroutine timestep), since it is physically impossible that the canopy
1092	!> drag alone can locally change the sign of a velocity component. This
1093	!> limitation is realized by calculating preliminary tendencies and velocities.
1094	!> It is subsequently checked if the preliminary new velocity has a different
1095	!> sign than the current velocity. If so, the tendency is limited in a way that
1096	!> the velocity can at maximum be reduced to zero by the canopy drag.
1097	!>
1098	!>
1099	!> Call for grid point i,j
1100	!------------------------------------------------------------------------------!
1101	SUBROUTINE pcm_tendency_ij( i, j, component )
1102
1103
1104	USE control_parameters, &
1105	ONLY: dt_3d, message_string
1106
1107	USE kinds
1108
1109	IMPLICIT NONE
1110
1111	INTEGER(iwp) :: component !< prognostic variable (u,v,w,pt,q,e)
1112	INTEGER(iwp) :: i !< running index
1113	INTEGER(iwp) :: j !< running index
1114	INTEGER(iwp) :: k !< running index
1115	INTEGER(iwp) :: kk !< running index for flat lad arrays
1116
1117	REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d)
1118	REAL(wp) :: lad_local !< local lad value
1119	REAL(wp) :: pre_tend !< preliminary tendency
1120	REAL(wp) :: pre_u !< preliminary u-value
1121	REAL(wp) :: pre_v !< preliminary v-value
1122	REAL(wp) :: pre_w !< preliminary w-value
1123
1124
1125	ddt_3d = 1.0_wp / dt_3d
1126
1127	!
1128	!-- Compute drag for the three velocity components and the SGS-TKE
1129	SELECT CASE ( component )
1130
1131	!
1132	!-- u-component
1133	CASE ( 1 )
1134	DO k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
1135
1136	kk = k - nzb_u_inner(j,i) !- lad arrays are defined flat
1137	!
1138	!-- In order to create sharp boundaries of the plant canopy,
1139	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1140	!-- rather than being interpolated from the surrounding lad_s,
1141	!-- because this would yield smaller lad at the canopy boundaries
1142	!-- than inside of the canopy.
1143	!-- For the same reason, the lad at the rightmost(i+1)canopy
1144	!-- boundary on the u-grid equals lad_s(k,j,i).
1145	lad_local = lad_s(kk,j,i)
1146	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp ) THEN
1147	lad_local = lad_s(kk,j,i-1)
1148	ENDIF
1149
1150	pre_tend = 0.0_wp
1151	pre_u = 0.0_wp
1152	!
1153	!-- Calculate preliminary value (pre_tend) of the tendency
1154	pre_tend = - cdc * &
1155	lad_local * &
1156	SQRT( u(k,j,i)**2 + &
1157	( 0.25_wp * ( v(k,j,i-1) + &
1158	v(k,j,i) + &
1159	v(k,j+1,i) + &
1160	v(k,j+1,i-1) ) &
1161	)**2 + &
1162	( 0.25_wp * ( w(k-1,j,i-1) + &
1163	w(k-1,j,i) + &
1164	w(k,j,i-1) + &
1165	w(k,j,i) ) &
1166	)**2 &
1167	) * &
1168	u(k,j,i)
1169
1170	!
1171	!-- Calculate preliminary new velocity, based on pre_tend
1172	pre_u = u(k,j,i) + dt_3d * pre_tend
1173	!
1174	!-- Compare sign of old velocity and new preliminary velocity,
1175	!-- and in case the signs are different, limit the tendency
1176	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
1177	pre_tend = - u(k,j,i) * ddt_3d
1178	ELSE
1179	pre_tend = pre_tend
1180	ENDIF
1181	!
1182	!-- Calculate final tendency
1183	tend(k,j,i) = tend(k,j,i) + pre_tend
1184	ENDDO
1185
1186
1187	!
1188	!-- v-component
1189	CASE ( 2 )
1190	DO k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
1191
1192	kk = k - nzb_v_inner(j,i) !- lad arrays are defined flat
1193	!
1194	!-- In order to create sharp boundaries of the plant canopy,
1195	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1196	!-- rather than being interpolated from the surrounding lad_s,
1197	!-- because this would yield smaller lad at the canopy boundaries
1198	!-- than inside of the canopy.
1199	!-- For the same reason, the lad at the northmost(j+1)canopy
1200	!-- boundary on the v-grid equals lad_s(k,j,i).
1201	lad_local = lad_s(kk,j,i)
1202	IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp ) THEN
1203	lad_local = lad_s(kk,j-1,i)
1204	ENDIF
1205
1206	pre_tend = 0.0_wp
1207	pre_v = 0.0_wp
1208	!
1209	!-- Calculate preliminary value (pre_tend) of the tendency
1210	pre_tend = - cdc * &
1211	lad_local * &
1212	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
1213	u(k,j-1,i+1) + &
1214	u(k,j,i) + &
1215	u(k,j,i+1) ) &
1216	)**2 + &
1217	v(k,j,i)**2 + &
1218	( 0.25_wp * ( w(k-1,j-1,i) + &
1219	w(k-1,j,i) + &
1220	w(k,j-1,i) + &
1221	w(k,j,i) ) &
1222	)**2 &
1223	) * &
1224	v(k,j,i)
1225
1226	!
1227	!-- Calculate preliminary new velocity, based on pre_tend
1228	pre_v = v(k,j,i) + dt_3d * pre_tend
1229	!
1230	!-- Compare sign of old velocity and new preliminary velocity,
1231	!-- and in case the signs are different, limit the tendency
1232	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
1233	pre_tend = - v(k,j,i) * ddt_3d
1234	ELSE
1235	pre_tend = pre_tend
1236	ENDIF
1237	!
1238	!-- Calculate final tendency
1239	tend(k,j,i) = tend(k,j,i) + pre_tend
1240	ENDDO
1241
1242
1243	!
1244	!-- w-component
1245	CASE ( 3 )
1246	DO k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
1247
1248	kk = k - nzb_w_inner(j,i) !- lad arrays are defined flat
1249
1250	pre_tend = 0.0_wp
1251	pre_w = 0.0_wp
1252	!
1253	!-- Calculate preliminary value (pre_tend) of the tendency
1254	pre_tend = - cdc * &
1255	(0.5_wp * &
1256	( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * &
1257	SQRT( ( 0.25_wp * ( u(k,j,i) + &
1258	u(k,j,i+1) + &
1259	u(k+1,j,i) + &
1260	u(k+1,j,i+1) ) &
1261	)**2 + &
1262	( 0.25_wp * ( v(k,j,i) + &
1263	v(k,j+1,i) + &
1264	v(k+1,j,i) + &
1265	v(k+1,j+1,i) ) &
1266	)**2 + &
1267	w(k,j,i)**2 &
1268	) * &
1269	w(k,j,i)
1270	!
1271	!-- Calculate preliminary new velocity, based on pre_tend
1272	pre_w = w(k,j,i) + dt_3d * pre_tend
1273	!
1274	!-- Compare sign of old velocity and new preliminary velocity,
1275	!-- and in case the signs are different, limit the tendency
1276	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
1277	pre_tend = - w(k,j,i) * ddt_3d
1278	ELSE
1279	pre_tend = pre_tend
1280	ENDIF
1281	!
1282	!-- Calculate final tendency
1283	tend(k,j,i) = tend(k,j,i) + pre_tend
1284	ENDDO
1285
1286	!
1287	!-- potential temperature
1288	CASE ( 4 )
1289	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1290	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1291	tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
1292	ENDDO
1293
1294
1295	!
1296	!-- humidity
1297	CASE ( 5 )
1298	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1299	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1300	tend(k,j,i) = tend(k,j,i) - &
1301	lsec * &
1302	lad_s(kk,j,i) * &
1303	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1304	u(k,j,i+1) ) &
1305	)**2 + &
1306	( 0.5_wp * ( v(k,j,i) + &
1307	v(k,j+1,i) ) &
1308	)**2 + &
1309	( 0.5_wp * ( w(k-1,j,i) + &
1310	w(k,j,i) ) &
1311	)**2 &
1312	) * &
1313	( q(k,j,i) - lsc )
1314	ENDDO
1315
1316	!
1317	!-- sgs-tke
1318	CASE ( 6 )
1319	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1320	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1321	tend(k,j,i) = tend(k,j,i) - &
1322	2.0_wp * cdc * &
1323	lad_s(kk,j,i) * &
1324	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1325	u(k,j,i+1) ) &
1326	)**2 + &
1327	( 0.5_wp * ( v(k,j,i) + &
1328	v(k,j+1,i) ) &
1329	)**2 + &
1330	( 0.5_wp * ( w(k,j,i) + &
1331	w(k+1,j,i) ) &
1332	)**2 &
1333	) * &
1334	e(k,j,i)
1335	ENDDO
1336
1337	!
1338	!-- scalar concentration
1339	CASE ( 7 )
1340	DO k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1341	kk = k - nzb_s_inner(j,i) !- lad arrays are defined flat
1342	tend(k,j,i) = tend(k,j,i) - &
1343	lsec * &
1344	lad_s(kk,j,i) * &
1345	SQRT( ( 0.5_wp * ( u(k,j,i) + &
1346	u(k,j,i+1) ) &
1347	)**2 + &
1348	( 0.5_wp * ( v(k,j,i) + &
1349	v(k,j+1,i) ) &
1350	)**2 + &
1351	( 0.5_wp * ( w(k-1,j,i) + &
1352	w(k,j,i) ) &
1353	)**2 &
1354	) * &
1355	( s(k,j,i) - lsc )
1356	ENDDO
1357
1358	CASE DEFAULT
1359
1360	WRITE( message_string, * ) 'wrong component: ', component
1361	CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
1362
1363	END SELECT
1364
1365	END SUBROUTINE pcm_tendency_ij
1366
1367
1368
1369	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 |