Home

Context Navigation

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

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