Home

Context Navigation

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

Last change on this file since 2825 was 2770, checked in by kanani, 7 years ago
Correction of parameter check
Property svn:keywords set to `Id`
File size: 71.5 KB

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