Home

Context Navigation

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

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

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