Home

Context Navigation

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

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