Home

Context Navigation

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

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