Home

Context Navigation

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

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