Home

Context Navigation

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

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