Home

Context Navigation

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

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