Home

Context Navigation

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

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