Home

Context Navigation

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

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