Home

Context Navigation

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

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