Home

Context Navigation

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

Last change on this file since 2209 was 2209, checked in by kanani, 7 years ago
small bugfix, formatting and new PCM output
Property svn:keywords set to `Id`
File size: 63.2 KB

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