Home

Context Navigation

source: palm/trunk/SOURCE/plant_canopy_model.f90 @ 1660

Last change on this file since 1660 was 1485, checked in by kanani, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 40.7 KB

Line
1	MODULE plant_canopy_model_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later 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-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: plant_canopy_model.f90 1485 2014-10-21 11:09:54Z gronemeier $
27	!
28	! 1484 2014-10-21 10:53:05Z kanani
29	! Changes due to new module structure of the plant canopy model:
30	! module plant_canopy_model_mod now contains a subroutine for the
31	! initialization of the canopy model (init_plant_canopy),
32	! limitation of the canopy drag (previously accounted for by calculation of
33	! a limiting canopy timestep for the determination of the maximum LES timestep
34	! in subroutine timestep) is now realized by the calculation of pre-tendencies
35	! and preliminary velocities in subroutine plant_canopy_model,
36	! some redundant MPI communication removed in subroutine init_plant_canopy
37	! (was previously in init_3d_model),
38	! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
39	! respective grid is now provided only by lad_s (e.g. in the calculation of
40	! the tendency terms or of cum_lai_hf),
41	! drag_coefficient, lai, leaf_surface_concentration,
42	! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
43	! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
44	! respectively,
45	! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
46	! USE-statements and ONLY-lists modified accordingly
47	!
48	! 1340 2014-03-25 19:45:13Z kanani
49	! REAL constants defined as wp-kind
50	!
51	! 1320 2014-03-20 08:40:49Z raasch
52	! ONLY-attribute added to USE-statements,
53	! kind-parameters added to all INTEGER and REAL declaration statements,
54	! kinds are defined in new module kinds,
55	! old module precision_kind is removed,
56	! revision history before 2012 removed,
57	! comment fields (!:) to be used for variable explanations added to
58	! all variable declaration statements
59	!
60	! 1036 2012-10-22 13:43:42Z raasch
61	! code put under GPL (PALM 3.9)
62	!
63	! 138 2007-11-28 10:03:58Z letzel
64	! Initial revision
65	!
66	! Description:
67	! ------------
68	! 1) Initialization of the canopy model, e.g. construction of leaf area density
69	! profile (subroutine init_plant_canopy).
70	! 2) Calculation of sinks and sources of momentum, heat and scalar concentration
71	! due to canopy elements (subroutine plant_canopy_model).
72	!------------------------------------------------------------------------------!
73	USE arrays_3d, &
74	ONLY: dzu, dzw, e, q, shf, tend, u, v, w, zu, zw
75
76	USE indices, &
77	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
78	nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
79
80	USE kinds
81
82
83	IMPLICIT NONE
84
85
86	CHARACTER (LEN=20) :: canopy_mode = 'block' !: canopy coverage
87
88	INTEGER(iwp) :: pch_index = 0 !: plant canopy height/top index
89	INTEGER(iwp) :: &
90	lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index)
91
92	LOGICAL :: calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func.
93	LOGICAL :: plant_canopy = .FALSE. !: switch for use of canopy model
94
95	REAL(wp) :: alpha_lad = 9999999.9_wp !: coefficient for lad calculation
96	REAL(wp) :: beta_lad = 9999999.9_wp !: coefficient for lad calculation
97	REAL(wp) :: canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter)
98	REAL(wp) :: cdc = 0.0_wp !: canopy drag coeff. (abbreviation used in equations)
99	REAL(wp) :: cthf = 0.0_wp !: canopy top heat flux
100	REAL(wp) :: dt_plant_canopy = 0.0_wp !: timestep account. for canopy drag
101	REAL(wp) :: lad_surface = 0.0_wp !: lad surface value
102	REAL(wp) :: lai_beta = 0.0_wp !: leaf area index (lai) for lad calc.
103	REAL(wp) :: &
104	leaf_scalar_exch_coeff = 0.0_wp !: canopy scalar exchange coeff.
105	REAL(wp) :: &
106	leaf_surface_conc = 0.0_wp !: leaf surface concentration
107	REAL(wp) :: lsec = 0.0_wp !: leaf scalar exchange coeff.
108	REAL(wp) :: lsc = 0.0_wp !: leaf surface concentration
109
110	REAL(wp) :: &
111	lad_vertical_gradient(10) = 0.0_wp !: lad gradient
112	REAL(wp) :: &
113	lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m)
114
115	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !: leaf area density
116	REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !: preliminary lad
117
118	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
119	canopy_heat_flux !: canopy heat flux
120	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !: cumulative lai for heatflux calc.
121	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !: lad on scalar-grid
122
123
124	SAVE
125
126
127	PRIVATE
128	PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, &
129	canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad, &
130	lad_s, lad_surface, lad_vertical_gradient, &
131	lad_vertical_gradient_level, lad_vertical_gradient_level_ind, &
132	lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec, &
133	pch_index, plant_canopy, plant_canopy_model
134
135
136	INTERFACE init_plant_canopy
137	MODULE PROCEDURE init_plant_canopy
138	END INTERFACE init_plant_canopy
139
140	INTERFACE plant_canopy_model
141	MODULE PROCEDURE plant_canopy_model
142	MODULE PROCEDURE plant_canopy_model_ij
143	END INTERFACE plant_canopy_model
144
145
146
147
148	CONTAINS
149
150
151	!------------------------------------------------------------------------------!
152	! Description:
153	! ------------
154	!-- Initialization of the plant canopy model
155	!------------------------------------------------------------------------------!
156	SUBROUTINE init_plant_canopy
157
158
159	USE control_parameters, &
160	ONLY: dz, ocean, passive_scalar
161
162
163	IMPLICIT NONE
164
165	INTEGER(iwp) :: i !: running index
166	INTEGER(iwp) :: j !: running index
167	INTEGER(iwp) :: k !: running index
168
169	REAL(wp) :: int_bpdf !: vertical integral for lad-profile construction
170	REAL(wp) :: dzh !: vertical grid spacing in units of canopy height
171	REAL(wp) :: gradient !: gradient for lad-profile construction
172	REAL(wp) :: canopy_height !: canopy height for lad-profile construction
173
174	!
175	!-- Allocate one-dimensional arrays for the computation of the
176	!-- leaf area density (lad) profile
177	ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
178	lad = 0.0_wp
179	pre_lad = 0.0_wp
180
181	!
182	!-- Compute the profile of leaf area density used in the plant
183	!-- canopy model. The profile can either be constructed from
184	!-- prescribed vertical gradients of the leaf area density or by
185	!-- using a beta probability density function (see e.g. Markkanen et al.,
186	!-- 2003: Boundary-Layer Meteorology, 106, 437-459)
187	IF ( .NOT. calc_beta_lad_profile ) THEN
188
189	!
190	!-- Use vertical gradients for lad-profile construction
191	i = 1
192	gradient = 0.0_wp
193
194	IF ( .NOT. ocean ) THEN
195
196	lad(0) = lad_surface
197	lad_vertical_gradient_level_ind(1) = 0
198
199	DO k = 1, pch_index
200	IF ( i < 11 ) THEN
201	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
202	lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
203	gradient = lad_vertical_gradient(i)
204	lad_vertical_gradient_level_ind(i) = k - 1
205	i = i + 1
206	ENDIF
207	ENDIF
208	IF ( gradient /= 0.0_wp ) THEN
209	IF ( k /= 1 ) THEN
210	lad(k) = lad(k-1) + dzu(k) * gradient
211	ELSE
212	lad(k) = lad_surface + dzu(k) * gradient
213	ENDIF
214	ELSE
215	lad(k) = lad(k-1)
216	ENDIF
217	ENDDO
218
219	ENDIF
220
221	!
222	!-- In case of no given leaf area density gradients, choose a vanishing
223	!-- gradient. This information is used for the HEADER and the RUN_CONTROL
224	!-- file.
225	IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
226	lad_vertical_gradient_level(1) = 0.0_wp
227	ENDIF
228
229	ELSE
230
231	!
232	!-- Use beta function for lad-profile construction
233	int_bpdf = 0.0_wp
234	canopy_height = pch_index * dz
235
236	DO k = nzb, pch_index
237	int_bpdf = int_bpdf + &
238	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
239	( ( 1.0_wp - ( zw(k) / canopy_height ) )*( beta_lad-1.0_wp ) ) &
240	( ( zw(k+1)-zw(k) ) / canopy_height ) )
241	ENDDO
242
243	!
244	!-- Preliminary lad profile (defined on w-grid)
245	DO k = nzb, pch_index
246	pre_lad(k) = lai_beta * &
247	( ( ( zw(k) / canopy_height )*( alpha_lad-1.0_wp ) ) &
248	( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) / &
249	int_bpdf &
250	) / canopy_height
251	ENDDO
252
253	!
254	!-- Final lad profile (defined on scalar-grid level, since most prognostic
255	!-- quantities are defined there, hence, less interpolation is required
256	!-- when calculating the canopy tendencies)
257	lad(0) = pre_lad(0)
258	DO k = nzb+1, pch_index
259	lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
260	ENDDO
261
262	ENDIF
263
264	!
265	!-- Allocate 3D-array for the leaf area density (lad_s). In case of a
266	!-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
267	!-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
268	ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
269
270	IF ( cthf /= 0.0_wp ) THEN
271	ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
272	canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
273	ENDIF
274
275	!
276	!-- Initialize canopy parameters cdc (canopy drag coefficient),
277	!-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
278	!-- with the prescribed values
279	cdc = canopy_drag_coeff
280	lsec = leaf_scalar_exch_coeff
281	lsc = leaf_surface_conc
282
283	!
284	!-- Initialization of the canopy coverage in the model domain:
285	!-- Setting the parameter canopy_mode = 'block' initializes a canopy, which
286	!-- fully covers the domain surface
287	SELECT CASE ( TRIM( canopy_mode ) )
288
289	CASE( 'block' )
290
291	DO i = nxlg, nxrg
292	DO j = nysg, nyng
293	lad_s(:,j,i) = lad(:)
294	ENDDO
295	ENDDO
296
297	CASE DEFAULT
298
299	!
300	!-- The DEFAULT case is reached either if the parameter
301	!-- canopy mode contains a wrong character string or if the
302	!-- user has coded a special case in the user interface.
303	!-- There, the subroutine user_init_plant_canopy checks
304	!-- which of these two conditions applies.
305	CALL user_init_plant_canopy
306
307	END SELECT
308
309	!
310	!-- Initialization of the canopy heat source distribution
311	IF ( cthf /= 0.0_wp ) THEN
312	!
313	!-- Piecewise calculation of the leaf area index by vertical
314	!-- integration of the leaf area density
315	cum_lai_hf(:,:,:) = 0.0_wp
316	DO i = nxlg, nxrg
317	DO j = nysg, nyng
318	DO k = pch_index-1, 0, -1
319	IF ( k == pch_index-1 ) THEN
320	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
321	( 0.5_wp * lad_s(k+1,j,i) * &
322	( zw(k+1) - zu(k+1) ) ) + &
323	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
324	lad_s(k,j,i) ) + &
325	lad_s(k+1,j,i) ) * &
326	( zu(k+1) - zw(k) ) )
327	ELSE
328	cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + &
329	( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + &
330	lad_s(k+1,j,i) ) + &
331	lad_s(k+1,j,i) ) * &
332	( zw(k+1) - zu(k+1) ) ) + &
333	( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + &
334	lad_s(k,j,i) ) + &
335	lad_s(k+1,j,i) ) * &
336	( zu(k+1) - zw(k) ) )
337	ENDIF
338	ENDDO
339	ENDDO
340	ENDDO
341
342	!
343	!-- Calculation of the upward kinematic vertical heat flux within the
344	!-- canopy
345	DO i = nxlg, nxrg
346	DO j = nysg, nyng
347	DO k = 0, pch_index
348	canopy_heat_flux(k,j,i) = cthf * &
349	exp( -0.6_wp * cum_lai_hf(k,j,i) )
350	ENDDO
351	ENDDO
352	ENDDO
353
354	!
355	!-- The surface heat flux is set to the surface value of the calculated
356	!-- in-canopy heat flux distribution
357	shf(:,:) = canopy_heat_flux(0,:,:)
358
359	ENDIF
360
361
362
363	END SUBROUTINE init_plant_canopy
364
365
366
367	!------------------------------------------------------------------------------!
368	! Description:
369	! ------------
370	!-- Calculation of the tendency terms, accounting for the effect of the plant
371	!-- canopy on momentum and scalar quantities.
372	!--
373	!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
374	!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
375	!-- The lad on the w-grid is vertically interpolated from the surrounding
376	!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
377	!-- k = pch_index. Here, the lad is zero.
378	!--
379	!-- The canopy drag must be limited (previously accounted for by calculation of
380	!-- a limiting canopy timestep for the determination of the maximum LES timestep
381	!-- in subroutine timestep), since it is physically impossible that the canopy
382	!-- drag alone can locally change the sign of a velocity component. This
383	!-- limitation is realized by calculating preliminary tendencies and velocities.
384	!-- It is subsequently checked if the preliminary new velocity has a different
385	!-- sign than the current velocity. If so, the tendency is limited in a way that
386	!-- the velocity can at maximum be reduced to zero by the canopy drag.
387	!--
388	!--
389	!-- Call for all grid points
390	!------------------------------------------------------------------------------!
391	SUBROUTINE plant_canopy_model( component )
392
393
394	USE control_parameters, &
395	ONLY: dt_3d, message_string
396
397	USE kinds
398
399	IMPLICIT NONE
400
401	INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e)
402	INTEGER(iwp) :: i !: running index
403	INTEGER(iwp) :: j !: running index
404	INTEGER(iwp) :: k !: running index
405
406	REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d)
407	REAL(wp) :: lad_local !: local lad value
408	REAL(wp) :: pre_tend !: preliminary tendency
409	REAL(wp) :: pre_u !: preliminary u-value
410	REAL(wp) :: pre_v !: preliminary v-value
411	REAL(wp) :: pre_w !: preliminary w-value
412
413
414	ddt_3d = 1.0_wp / dt_3d
415
416	!
417	!-- Compute drag for the three velocity components and the SGS-TKE:
418	SELECT CASE ( component )
419
420	!
421	!-- u-component
422	CASE ( 1 )
423	DO i = nxlu, nxr
424	DO j = nys, nyn
425	DO k = nzb_u_inner(j,i)+1, pch_index
426
427	!
428	!-- In order to create sharp boundaries of the plant canopy,
429	!-- the lad on the u-grid at index (k,j,i) is equal to
430	!-- lad_s(k,j,i), rather than being interpolated from the
431	!-- surrounding lad_s, because this would yield smaller lad
432	!-- at the canopy boundaries than inside of the canopy.
433	!-- For the same reason, the lad at the rightmost(i+1)canopy
434	!-- boundary on the u-grid equals lad_s(k,j,i).
435	lad_local = lad_s(k,j,i)
436	IF ( lad_local == 0.0_wp .AND. &
437	lad_s(k,j,i-1) > 0.0_wp ) THEN
438	lad_local = lad_s(k,j,i-1)
439	ENDIF
440
441	pre_tend = 0.0_wp
442	pre_u = 0.0_wp
443	!
444	!-- Calculate preliminary value (pre_tend) of the tendency
445	pre_tend = - cdc * &
446	lad_local * &
447	SQRT( u(k,j,i)**2 + &
448	( 0.25_wp * ( v(k,j,i-1) + &
449	v(k,j,i) + &
450	v(k,j+1,i) + &
451	v(k,j+1,i-1) ) &
452	)**2 + &
453	( 0.25_wp * ( w(k-1,j,i-1) + &
454	w(k-1,j,i) + &
455	w(k,j,i-1) + &
456	w(k,j,i) ) &
457	)**2 &
458	) * &
459	u(k,j,i)
460
461	!
462	!-- Calculate preliminary new velocity, based on pre_tend
463	pre_u = u(k,j,i) + dt_3d * pre_tend
464	!
465	!-- Compare sign of old velocity and new preliminary velocity,
466	!-- and in case the signs are different, limit the tendency
467	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
468	pre_tend = - u(k,j,i) * ddt_3d
469	ELSE
470	pre_tend = pre_tend
471	ENDIF
472	!
473	!-- Calculate final tendency
474	tend(k,j,i) = tend(k,j,i) + pre_tend
475
476	ENDDO
477	ENDDO
478	ENDDO
479
480	!
481	!-- v-component
482	CASE ( 2 )
483	DO i = nxl, nxr
484	DO j = nysv, nyn
485	DO k = nzb_v_inner(j,i)+1, pch_index
486
487	!
488	!-- In order to create sharp boundaries of the plant canopy,
489	!-- the lad on the v-grid at index (k,j,i) is equal to
490	!-- lad_s(k,j,i), rather than being interpolated from the
491	!-- surrounding lad_s, because this would yield smaller lad
492	!-- at the canopy boundaries than inside of the canopy.
493	!-- For the same reason, the lad at the northmost(j+1) canopy
494	!-- boundary on the v-grid equals lad_s(k,j,i).
495	lad_local = lad_s(k,j,i)
496	IF ( lad_local == 0.0_wp .AND. &
497	lad_s(k,j-1,i) > 0.0_wp ) THEN
498	lad_local = lad_s(k,j-1,i)
499	ENDIF
500
501	pre_tend = 0.0_wp
502	pre_v = 0.0_wp
503	!
504	!-- Calculate preliminary value (pre_tend) of the tendency
505	pre_tend = - cdc * &
506	lad_local * &
507	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
508	u(k,j-1,i+1) + &
509	u(k,j,i) + &
510	u(k,j,i+1) ) &
511	)**2 + &
512	v(k,j,i)**2 + &
513	( 0.25_wp * ( w(k-1,j-1,i) + &
514	w(k-1,j,i) + &
515	w(k,j-1,i) + &
516	w(k,j,i) ) &
517	)**2 &
518	) * &
519	v(k,j,i)
520
521	!
522	!-- Calculate preliminary new velocity, based on pre_tend
523	pre_v = v(k,j,i) + dt_3d * pre_tend
524	!
525	!-- Compare sign of old velocity and new preliminary velocity,
526	!-- and in case the signs are different, limit the tendency
527	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
528	pre_tend = - v(k,j,i) * ddt_3d
529	ELSE
530	pre_tend = pre_tend
531	ENDIF
532	!
533	!-- Calculate final tendency
534	tend(k,j,i) = tend(k,j,i) + pre_tend
535
536	ENDDO
537	ENDDO
538	ENDDO
539
540	!
541	!-- w-component
542	CASE ( 3 )
543	DO i = nxl, nxr
544	DO j = nys, nyn
545	DO k = nzb_w_inner(j,i)+1, pch_index-1
546
547	pre_tend = 0.0_wp
548	pre_w = 0.0_wp
549	!
550	!-- Calculate preliminary value (pre_tend) of the tendency
551	pre_tend = - cdc * &
552	(0.5_wp * &
553	( lad_s(k+1,j,i) + lad_s(k,j,i) )) * &
554	SQRT( ( 0.25_wp * ( u(k,j,i) + &
555	u(k,j,i+1) + &
556	u(k+1,j,i) + &
557	u(k+1,j,i+1) ) &
558	)**2 + &
559	( 0.25_wp * ( v(k,j,i) + &
560	v(k,j+1,i) + &
561	v(k+1,j,i) + &
562	v(k+1,j+1,i) ) &
563	)**2 + &
564	w(k,j,i)**2 &
565	) * &
566	w(k,j,i)
567	!
568	!-- Calculate preliminary new velocity, based on pre_tend
569	pre_w = w(k,j,i) + dt_3d * pre_tend
570	!
571	!-- Compare sign of old velocity and new preliminary velocity,
572	!-- and in case the signs are different, limit the tendency
573	IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN
574	pre_tend = - w(k,j,i) * ddt_3d
575	ELSE
576	pre_tend = pre_tend
577	ENDIF
578	!
579	!-- Calculate final tendency
580	tend(k,j,i) = tend(k,j,i) + pre_tend
581
582	ENDDO
583	ENDDO
584	ENDDO
585
586	!
587	!-- potential temperature
588	CASE ( 4 )
589	DO i = nxl, nxr
590	DO j = nys, nyn
591	DO k = nzb_s_inner(j,i)+1, pch_index
592	tend(k,j,i) = tend(k,j,i) + &
593	( canopy_heat_flux(k,j,i) - &
594	canopy_heat_flux(k-1,j,i) ) / dzw(k)
595	ENDDO
596	ENDDO
597	ENDDO
598
599	!
600	!-- scalar concentration
601	CASE ( 5 )
602	DO i = nxl, nxr
603	DO j = nys, nyn
604	DO k = nzb_s_inner(j,i)+1, pch_index
605	tend(k,j,i) = tend(k,j,i) - &
606	lsec * &
607	lad_s(k,j,i) * &
608	SQRT( ( 0.5_wp * ( u(k,j,i) + &
609	u(k,j,i+1) ) &
610	)**2 + &
611	( 0.5_wp * ( v(k,j,i) + &
612	v(k,j+1,i) ) &
613	)**2 + &
614	( 0.5_wp * ( w(k-1,j,i) + &
615	w(k,j,i) ) &
616	)**2 &
617	) * &
618	( q(k,j,i) - lsc )
619	ENDDO
620	ENDDO
621	ENDDO
622
623	!
624	!-- sgs-tke
625	CASE ( 6 )
626	DO i = nxl, nxr
627	DO j = nys, nyn
628	DO k = nzb_s_inner(j,i)+1, pch_index
629	tend(k,j,i) = tend(k,j,i) - &
630	2.0_wp * cdc * &
631	lad_s(k,j,i) * &
632	SQRT( ( 0.5_wp * ( u(k,j,i) + &
633	u(k,j,i+1) ) &
634	)**2 + &
635	( 0.5_wp * ( v(k,j,i) + &
636	v(k,j+1,i) ) &
637	)**2 + &
638	( 0.5_wp * ( w(k,j,i) + &
639	w(k+1,j,i) ) &
640	)**2 &
641	) * &
642	e(k,j,i)
643	ENDDO
644	ENDDO
645	ENDDO
646
647
648	CASE DEFAULT
649
650	WRITE( message_string, * ) 'wrong component: ', component
651	CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
652
653	END SELECT
654
655	END SUBROUTINE plant_canopy_model
656
657
658	!------------------------------------------------------------------------------!
659	! Description:
660	! ------------
661	!-- Calculation of the tendency terms, accounting for the effect of the plant
662	!-- canopy on momentum and scalar quantities.
663	!--
664	!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
665	!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy.
666	!-- The lad on the w-grid is vertically interpolated from the surrounding
667	!-- lad_s. The upper boundary of the canopy is defined on the w-grid at
668	!-- k = pch_index. Here, the lad is zero.
669	!--
670	!-- The canopy drag must be limited (previously accounted for by calculation of
671	!-- a limiting canopy timestep for the determination of the maximum LES timestep
672	!-- in subroutine timestep), since it is physically impossible that the canopy
673	!-- drag alone can locally change the sign of a velocity component. This
674	!-- limitation is realized by calculating preliminary tendencies and velocities.
675	!-- It is subsequently checked if the preliminary new velocity has a different
676	!-- sign than the current velocity. If so, the tendency is limited in a way that
677	!-- the velocity can at maximum be reduced to zero by the canopy drag.
678	!--
679	!--
680	!-- Call for grid point i,j
681	!------------------------------------------------------------------------------!
682	SUBROUTINE plant_canopy_model_ij( i, j, component )
683
684
685	USE control_parameters, &
686	ONLY: dt_3d, message_string
687
688	USE kinds
689
690	IMPLICIT NONE
691
692	INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e)
693	INTEGER(iwp) :: i !: running index
694	INTEGER(iwp) :: j !: running index
695	INTEGER(iwp) :: k !: running index
696
697	REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d)
698	REAL(wp) :: lad_local !: local lad value
699	REAL(wp) :: pre_tend !: preliminary tendency
700	REAL(wp) :: pre_u !: preliminary u-value
701	REAL(wp) :: pre_v !: preliminary v-value
702	REAL(wp) :: pre_w !: preliminary w-value
703
704
705	ddt_3d = 1.0_wp / dt_3d
706
707	!
708	!-- Compute drag for the three velocity components and the SGS-TKE
709	SELECT CASE ( component )
710
711	!
712	!-- u-component
713	CASE ( 1 )
714	DO k = nzb_u_inner(j,i)+1, pch_index
715
716	!
717	!-- In order to create sharp boundaries of the plant canopy,
718	!-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
719	!-- rather than being interpolated from the surrounding lad_s,
720	!-- because this would yield smaller lad at the canopy boundaries
721	!-- than inside of the canopy.
722	!-- For the same reason, the lad at the rightmost(i+1)canopy
723	!-- boundary on the u-grid equals lad_s(k,j,i).
724	lad_local = lad_s(k,j,i)
725	IF ( lad_local == 0.0_wp .AND. &
726	lad_s(k,j,i-1) > 0.0_wp ) THEN
727	lad_local = lad_s(k,j,i-1)
728	ENDIF
729
730	pre_tend = 0.0_wp
731	pre_u = 0.0_wp
732	!
733	!-- Calculate preliminary value (pre_tend) of the tendency
734	pre_tend = - cdc * &
735	lad_local * &
736	SQRT( u(k,j,i)**2 + &
737	( 0.25_wp * ( v(k,j,i-1) + &
738	v(k,j,i) + &
739	v(k,j+1,i) + &
740	v(k,j+1,i-1) ) &
741	)**2 + &
742	( 0.25_wp * ( w(k-1,j,i-1) + &
743	w(k-1,j,i) + &
744	w(k,j,i-1) + &
745	w(k,j,i) ) &
746	)**2 &
747	) * &
748	u(k,j,i)
749
750	!
751	!-- Calculate preliminary new velocity, based on pre_tend
752	pre_u = u(k,j,i) + dt_3d * pre_tend
753	!
754	!-- Compare sign of old velocity and new preliminary velocity,
755	!-- and in case the signs are different, limit the tendency
756	IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN
757	pre_tend = - u(k,j,i) * ddt_3d
758	ELSE
759	pre_tend = pre_tend
760	ENDIF
761	!
762	!-- Calculate final tendency
763	tend(k,j,i) = tend(k,j,i) + pre_tend
764	ENDDO
765
766
767	!
768	!-- v-component
769	CASE ( 2 )
770	DO k = nzb_v_inner(j,i)+1, pch_index
771
772	!
773	!-- In order to create sharp boundaries of the plant canopy,
774	!-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
775	!-- rather than being interpolated from the surrounding lad_s,
776	!-- because this would yield smaller lad at the canopy boundaries
777	!-- than inside of the canopy.
778	!-- For the same reason, the lad at the northmost(j+1)canopy
779	!-- boundary on the v-grid equals lad_s(k,j,i).
780	lad_local = lad_s(k,j,i)
781	IF ( lad_local == 0.0_wp .AND. &
782	lad_s(k,j-1,i) > 0.0_wp ) THEN
783	lad_local = lad_s(k,j-1,i)
784	ENDIF
785
786	pre_tend = 0.0_wp
787	pre_v = 0.0_wp
788	!
789	!-- Calculate preliminary value (pre_tend) of the tendency
790	pre_tend = - cdc * &
791	lad_local * &
792	SQRT( ( 0.25_wp * ( u(k,j-1,i) + &
793	u(k,j-1,i+1) + &
794	u(k,j,i) + &
795	u(k,j,i+1) ) &
796	)**2 + &
797	v(k,j,i)**2 + &
798	( 0.25_wp * ( w(k-1,j-1,i) + &
799	w(k-1,j,i) + &
800	w(k,j-1,i) + &
801	w(k,j,i) ) &
802	)**2 &
803	) * &
804	v(k,j,i)
805
806	!
807	!-- Calculate preliminary new velocity, based on pre_tend
808	pre_v = v(k,j,i) + dt_3d * pre_tend
809	!
810	!-- Compare sign of old velocity and new preliminary velocity,
811	!-- and in case the signs are different, limit the tendency
812	IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN
813	pre_tend = - v(k,j,i) * ddt_3d
814	ELSE
815	pre_tend = pre_tend
816	ENDIF
817	!
818	!-- Calculate final tendency
819	tend(k,j,i) = tend(k,j,i) + pre_tend
820	ENDDO
821
822
823	!
824	!-- w-component
825	CASE ( 3 )
826	DO k = nzb_w_inner(j,i)+1, pch_index-1
827
828	pre_tend = 0.0_wp
829	pre_w = 0.0_wp
830	!
831	!-- Calculate preliminary value (pre_tend) of the tendency
832	pre_tend = - cdc * &
833	(0.5_wp * &
834	( lad_s(k+1,j,i) + lad_s(k,j,i) )) * &
835	SQRT( ( 0.25_wp * ( u(k,j,i) + &
836	u(k,j,i+1) + &
837	u(k+1,j,i) + &
838	u(k+1,j,i+1) ) &
839	)**2 + &
840	( 0.25_wp * ( v(k,j,i) + &
841	v(k,j+1,i) + &
842	v(k+1,j,i) + &
843	v(k+1,j+1,i) ) &
844	)**2 + &
845	w(k,j,i)**2 &
846	) * &
847	w(k,j,i)
848	!
849	!-- Calculate preliminary new velocity, based on pre_tend
850	pre_w = w(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_w,w(k,j,i)) /= pre_w ) THEN
855	pre_tend = - w(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	ENDDO
863
864	!
865	!-- potential temperature
866	CASE ( 4 )
867	DO k = nzb_s_inner(j,i)+1, pch_index
868	tend(k,j,i) = tend(k,j,i) + &
869	( canopy_heat_flux(k,j,i) - &
870	canopy_heat_flux(k-1,j,i) ) / dzw(k)
871	ENDDO
872
873
874	!
875	!-- scalar concentration
876	CASE ( 5 )
877	DO k = nzb_s_inner(j,i)+1, pch_index
878	tend(k,j,i) = tend(k,j,i) - &
879	lsec * &
880	lad_s(k,j,i) * &
881	SQRT( ( 0.5_wp * ( u(k,j,i) + &
882	u(k,j,i+1) ) &
883	)**2 + &
884	( 0.5_wp * ( v(k,j,i) + &
885	v(k,j+1,i) ) &
886	)**2 + &
887	( 0.5_wp * ( w(k-1,j,i) + &
888	w(k,j,i) ) &
889	)**2 &
890	) * &
891	( q(k,j,i) - lsc )
892	ENDDO
893
894	!
895	!-- sgs-tke
896	CASE ( 6 )
897	DO k = nzb_s_inner(j,i)+1, pch_index
898	tend(k,j,i) = tend(k,j,i) - &
899	2.0_wp * cdc * &
900	lad_s(k,j,i) * &
901	SQRT( ( 0.5_wp * ( u(k,j,i) + &
902	u(k,j,i+1) ) &
903	)**2 + &
904	( 0.5_wp * ( v(k,j,i) + &
905	v(k,j+1,i) ) &
906	)**2 + &
907	( 0.5_wp * ( w(k,j,i) + &
908	w(k+1,j,i) ) &
909	)**2 &
910	) * &
911	e(k,j,i)
912	ENDDO
913
914	CASE DEFAULT
915
916	WRITE( message_string, * ) 'wrong component: ', component
917	CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
918
919	END SELECT
920
921	END SUBROUTINE plant_canopy_model_ij
922
923	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 |