Home

Context Navigation

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

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