Home

Context Navigation

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

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