Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/inifor.f90 @ 3613

Last change on this file since 3613 was 3613, checked in by eckhard, 6 years ago
inifor: Average initial profiles over the PALM, not the geostrophic, region
Property svn:keywords set to `Id`
File size: 22.0 KB

Line
1	!> @file src/inifor.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 2017-2018 Leibniz Universitaet Hannover
18	! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: inifor.f90 3613 2018-12-07 18:20:37Z eckhard $
28	! Moved version output to setup_parameters()
29	!
30	! 3557 2018-11-22 16:01:22Z eckhard
31	! Updated documentation
32	!
33	! 3537 2018-11-20 10:53:14Z eckhard
34	! Print version number on program start
35	!
36	!
37	! 3456 2018-10-30 14:29:54Z eckhard
38	! NetCDf output of internal arrays only with --debug option
39	!
40	!
41	! 3401 2018-10-23 09:33:31Z eckhard
42	! Re-compute geostrophic winds every time step
43	!
44	! 3395 2018-10-22 17:32:49Z eckhard
45	! Added main loop support for computation of geostrophic winds and surface
46	! pressure
47	! Cleaned up terminal output, show some meVssages only with --debug option
48	!
49	! 3183 2018-07-27 14:25:55Z suehring
50	! Introduced new PALM grid stretching
51	! Renamend initial-condition mode variable 'mode' to 'ic_mode'
52	! Improved log messages
53	!
54	!
55	! 3182 2018-07-27 13:36:03Z suehring
56	! Initial revision
57	!
58	!
59	!
60	! Authors:
61	! --------
62	!> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach)
63	!
64	! Description:
65	! ------------
66	!> INIFOR is an interpolation tool for generating meteorological initialization
67	!> and forcing data for the urban climate model PALM-4U. The required
68	!> meteorological fields are interpolated from output data of the mesoscale
69	!> model COSMO-DE. This is the main program file.
70	!------------------------------------------------------------------------------!
71	PROGRAM inifor
72
73	USE control
74	USE defs
75	USE grid, &
76	ONLY: setup_parameters, setup_grids, setup_variable_tables, &
77	setup_io_groups, fini_grids, fini_variables, fini_io_groups, &
78	fini_file_lists, preprocess, origin_lon, origin_lat, &
79	output_file, io_group_list, output_var_table, &
80	cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3, &
81	averaging_width_ns, averaging_width_ew, phi_n, lambda_n, &
82	lam_centre, phi_centre
83	USE io
84	USE transform, &
85	ONLY: average_profile, interpolate_2d, interpolate_3d, &
86	geostrophic_winds, extrapolate_density, extrapolate_pressure, &
87	get_surface_pressure
88	USE types
89
90	IMPLICIT NONE
91
92	INTEGER :: igroup !< loop index for IO groups
93	INTEGER :: ivar !< loop index for output variables
94	INTEGER :: iter !< loop index for time steps
95
96	REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr !< array buffer for interpolated quantities
97	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre !< density profile of the centre averaging domain
98	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_arr !< geostrophic wind in x direction
99	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_arr !< geostrophic wind in y direction
100	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north !< density profile of the northern averaging domain
101	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south !< density profile of the southern averaging domain
102	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east !< density profile of the eastern averaging domain
103	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west !< density profile of the western averaging domain
104	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north !< pressure profile of the northern averaging domain
105	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south !< pressure profile of the southern averaging domain
106	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east !< pressure profile of the eastern averaging domain
107	REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west !< pressure profile of the western averaging domain
108
109	REAL(dp), POINTER, DIMENSION(:) :: internal_arr !< pointer to the currently processed internal array (density, pressure)
110	REAL(dp), POINTER, DIMENSION(:) :: ug_vg_arr !< pointer to the currently processed geostrophic wind component
111
112	TYPE(nc_var), POINTER :: output_var !< pointer to the currently processed output variable
113	TYPE(io_group), POINTER :: group !< pointer to the currently processed IO group
114	TYPE(container), ALLOCATABLE :: input_buffer(:) !< buffer of the current IO group's input arrays
115
116	LOGICAL, SAVE :: ug_vg_have_been_computed = .FALSE. !< flag for managing geostrophic wind allocation and computation
117	LOGICAL, SAVE :: debugging_output = .TRUE. !< flag controllging output of internal variables
118
119	!> \mainpage About INIFOR
120	!> ...
121	!
122	!------------------------------------------------------------------------------
123	!- Section 1: Initialization
124	!------------------------------------------------------------------------------
125	CALL run_control('init', 'void')
126
127	!
128	!-- Initialize INIFOR's parameters from command-line interface and namelists
129	CALL setup_parameters()
130
131	!
132	!-- Initialize all grids, including interpolation neighbours and weights
133	CALL setup_grids()
134	CALL run_control('time', 'init')
135
136	!
137	!-- Initialize the netCDF output file and define dimensions
138	CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date, &
139	origin_lon, origin_lat)
140	CALL run_control('time', 'write')
141
142	!
143	!-- Set up the tables containing the input and output variables and set
144	!-- the corresponding netCDF dimensions for each output variable
145	CALL setup_variable_tables(cfg % ic_mode)
146	CALL run_control('time', 'write')
147
148	!
149	!-- Add the output variables to the netCDF output file
150	CALL setup_netcdf_variables(output_file % name, output_var_table, &
151	cfg % debug)
152
153	CALL setup_io_groups()
154	CALL run_control('time', 'init')
155
156	!------------------------------------------------------------------------------
157	!- Section 2: Main loop
158	!------------------------------------------------------------------------------
159	DO igroup = 1, SIZE(io_group_list)
160
161	group => io_group_list(igroup)
162	IF ( group % to_be_processed ) THEN
163
164	DO iter = 1, group % nt
165
166	!------------------------------------------------------------------------------
167	!- Section 2.1: Read and preprocess input data
168	!------------------------------------------------------------------------------
169	CALL read_input_variables(group, iter, input_buffer)
170	CALL run_control('time', 'read')
171
172	CALL preprocess(group, input_buffer, cosmo_grid, iter)
173	CALL run_control('time', 'comp')
174
175	!TODO: move this assertion into 'preprocess'.
176	IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) ) THEN
177	message = "Input buffers for group '" // TRIM(group % kind) // &
178	"' could not be preprocessed sucessfully."
179	CALL abort('main loop', message)
180	END IF
181
182	!------------------------------------------------------------------------------
183	!- Section 2.2: Interpolate each output variable of the group
184	!------------------------------------------------------------------------------
185	DO ivar = 1, group % nv
186
187	output_var => group % out_vars( ivar )
188
189	IF ( output_var % to_be_processed .AND. &
190	iter .LE. output_var % nt ) THEN
191
192	message = "Processing '" // TRIM(output_var % name) // &
193	"' (" // TRIM(output_var % kind) // &
194	"), iteration " // TRIM(str(iter)) //" of " // &
195	TRIM(str(output_var % nt))
196	CALL report('main loop', message)
197
198	SELECT CASE( TRIM(output_var % task) )
199
200	CASE( 'interpolate_2d' )
201
202	SELECT CASE( TRIM(output_var % kind) )
203
204	CASE( 'init soil' )
205
206	ALLOCATE( output_arr( 0:output_var % grid % nx, &
207	0:output_var % grid % ny, &
208	SIZE(output_var % grid % depths) ) )
209
210	CASE ( 'surface forcing' )
211
212	ALLOCATE( output_arr( 0:output_var % grid % nx, &
213	0:output_var % grid % ny, 1 ) )
214
215	CASE DEFAULT
216
217	message = "'" // TRIM(output_var % kind) // "' is not a soil variable"
218	CALL abort("main loop", message)
219
220	END SELECT
221	CALL run_control('time', 'alloc')
222
223	CALL interpolate_2d(input_buffer(output_var % input_id) % array(:,:,:), &
224	output_arr(:,:,:), output_var % intermediate_grid, output_var)
225	CALL run_control('time', 'comp')
226
227
228	CASE ( 'interpolate_3d' )
229
230	ALLOCATE( output_arr( 0:output_var % grid % nx, &
231	0:output_var % grid % ny, &
232	1:output_var % grid % nz ) )
233
234	CALL run_control('time', 'alloc')
235	CALL interpolate_3d( &
236	input_buffer(output_var % input_id) % array(:,:,:), &
237	output_arr(:,:,:), &
238	output_var % intermediate_grid, &
239	output_var % grid)
240	CALL run_control('time', 'comp')
241
242	CASE ( 'average profile' )
243
244	ALLOCATE( output_arr( 0:output_var % grid % nx, &
245	0:output_var % grid % ny, &
246	1:output_var % grid % nz ) )
247	CALL run_control('time', 'alloc')
248
249	CALL average_profile( &
250	input_buffer(output_var % input_id) % array(:,:,:), &
251	output_arr(0,0,:), &
252	output_var % averaging_grid)
253
254	IF ( TRIM(output_var % name) == &
255	'surface_forcing_surface_pressure' ) THEN
256
257	IF ( cfg % p0_is_set ) THEN
258	output_arr(0,0,1) = p0
259	ELSE
260	CALL get_surface_pressure( &
261	output_arr(0,0,:), rho_centre, &
262	output_var % averaging_grid)
263	END IF
264
265	END IF
266	CALL run_control('time', 'comp')
267
268	CASE ( 'internal profile' )
269
270	message = "Averaging of internal profile for variable '" //&
271	TRIM(output_var % name) // "' is not supported."
272
273	SELECT CASE (TRIM(output_var % name))
274
275	CASE('internal_density_centre')
276	ALLOCATE( rho_centre( 1:output_var % grid % nz ) )
277	internal_arr => rho_centre
278
279	CASE('internal_density_north')
280	ALLOCATE( rho_north( 1:output_var % grid % nz ) )
281	internal_arr => rho_north
282
283	CASE('internal_density_south')
284	ALLOCATE( rho_south( 1:output_var % grid % nz ) )
285	internal_arr => rho_south
286
287	CASE('internal_density_east')
288	ALLOCATE( rho_east( 1:output_var % grid % nz) )
289	internal_arr => rho_east
290
291	CASE('internal_density_west')
292	ALLOCATE( rho_west( 1:output_var % grid % nz ) )
293	internal_arr => rho_west
294
295	CASE('internal_pressure_north')
296	ALLOCATE( p_north( 1:output_var % grid % nz ) )
297	internal_arr => p_north
298
299	CASE('internal_pressure_south')
300	ALLOCATE( p_south( 1:output_var % grid % nz ) )
301	internal_arr => p_south
302
303	CASE('internal_pressure_east')
304	ALLOCATE( p_east( 1:output_var % grid % nz) )
305	internal_arr => p_east
306
307	CASE('internal_pressure_west')
308	ALLOCATE( p_west( 1:output_var % grid % nz ) )
309	internal_arr => p_west
310
311	CASE DEFAULT
312	CALL abort('main loop', message)
313
314	END SELECT
315	CALL run_control('time', 'alloc')
316
317
318	CALL average_profile( &
319	input_buffer(output_var % input_id) % array(:,:,:),&
320	internal_arr(:), &
321	output_var % averaging_grid)
322
323	SELECT CASE (TRIM(output_var % name))
324
325	CASE('internal_density_centre', &
326	'internal_density_north', &
327	'internal_density_south', &
328	'internal_density_east', &
329	'internal_density_west')
330	CALL extrapolate_density(internal_arr, &
331	output_var % averaging_grid)
332
333	CASE('internal_pressure_north')
334	CALL extrapolate_pressure(internal_arr, rho_north, &
335	output_var % averaging_grid)
336
337	CASE('internal_pressure_south')
338	CALL extrapolate_pressure(internal_arr, rho_south, &
339	output_var % averaging_grid)
340
341	CASE('internal_pressure_east')
342	CALL extrapolate_pressure(internal_arr, rho_east, &
343	output_var % averaging_grid)
344
345	CASE('internal_pressure_west')
346	CALL extrapolate_pressure(internal_arr, rho_west, &
347	output_var % averaging_grid)
348
349	CASE DEFAULT
350	CALL abort('main loop', message)
351
352	END SELECT
353
354	IF (.TRUE.) THEN
355	ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
356	output_arr(1,1,:) = internal_arr(:)
357	END IF
358	CALL run_control('time', 'comp')
359
360	!
361	!-- This case gets called twice, the first time for ug, the
362	!-- second time for vg. We compute ug and vg at the first call
363	!-- and keep vg (and ug for that matter) around for the second
364	!-- call.
365	CASE ( 'geostrophic winds' )
366
367	IF (.NOT. ug_vg_have_been_computed ) THEN
368	ALLOCATE( ug_arr(output_var % grid % nz) )
369	ALLOCATE( vg_arr(output_var % grid % nz) )
370
371	IF ( cfg % ug_is_set ) THEN
372	ug_arr = cfg % ug
373	vg_arr = cfg % vg
374	ELSE
375	CALL geostrophic_winds( p_north, p_south, p_east, &
376	p_west, rho_centre, f3, &
377	averaging_width_ew, &
378	averaging_width_ns, &
379	phi_n, lambda_n, &
380	phi_centre, lam_centre, &
381	ug_arr, vg_arr )
382	END IF
383
384	ug_vg_have_been_computed = .TRUE.
385
386	END IF
387
388	!
389	!-- Prepare output of geostrophic winds
390	SELECT CASE(TRIM(output_var % name))
391	CASE ('ls_forcing_ug')
392	ug_vg_arr => ug_arr
393	CASE ('ls_forcing_vg')
394	ug_vg_arr => vg_arr
395	END SELECT
396
397	ALLOCATE( output_arr(1,1,output_var % grid % nz) )
398	output_arr(1,1,:) = ug_vg_arr(:)
399
400	CASE ( 'average scalar' )
401
402	ALLOCATE( output_arr(1,1,1) )
403	CALL run_control('time', 'alloc')
404	output_arr(1,1,1) = p0
405	CALL run_control('time', 'comp')
406
407	CASE ( 'set profile' )
408
409	ALLOCATE( output_arr( 1, 1, 1:nz ) )
410	CALL run_control('time', 'alloc')
411
412	SELECT CASE (TRIM(output_var % name))
413
414	CASE('nudging_tau')
415	output_arr(1, 1, :) = NUDGING_TAU
416
417	CASE DEFAULT
418	message = "'" // TRIM(output_var % name) // &
419	"' is not a valid '" // TRIM(output_var % kind) //&
420	"' variable kind."
421	CALL abort('main loop', message)
422	END SELECT
423	CALL run_control('time', 'comp')
424
425	CASE('average large-scale profile')
426	message = "Averaging of large-scale forcing profiles " //&
427	"has not been implemented, yet."
428	CALL abort('main loop', message)
429
430	CASE DEFAULT
431	message = "Processing task '" // TRIM(output_var % task) //&
432	"' not recognized."
433	CALL abort('', message)
434
435	END SELECT
436	CALL run_control('time', 'comp')
437
438	!------------------------------------------------------------------------------
439	!- Section 2.3: Write current time step of current variable
440	!------------------------------------------------------------------------------
441	IF (.NOT. output_var % is_internal .OR. debugging_output) THEN
442	message = "Writing variable '" // TRIM(output_var%name) // "'."
443	CALL report('main loop', message)
444	CALL update_output(output_var, output_arr, iter, &
445	output_file, cfg)
446	CALL run_control('time', 'write')
447	END IF
448
449	IF (ALLOCATED(output_arr)) DEALLOCATE(output_arr)
450	CALL run_control('time', 'alloc')
451
452	END IF
453
454	!
455	!-- output variable loop
456	END DO
457
458	ug_vg_have_been_computed = .FALSE.
459	IF ( group % kind == 'thermodynamics' ) THEN
460	DEALLOCATE( rho_centre )
461	DEALLOCATE( ug_arr )
462	DEALLOCATE( vg_arr )
463	IF ( .NOT. cfg % ug_is_set ) THEN
464	DEALLOCATE( rho_north )
465	DEALLOCATE( rho_south )
466	DEALLOCATE( rho_east )
467	DEALLOCATE( rho_west )
468	DEALLOCATE( p_north )
469	DEALLOCATE( p_south )
470	DEALLOCATE( p_east )
471	DEALLOCATE( p_west )
472	END IF
473	END IF
474
475	!
476	!-- Keep input buffer around for averaged (radiation) and
477	!-- accumulated COSMO quantities (precipitation).
478	IF ( group % kind == 'running average' .OR. &
479	group % kind == 'accumulated' ) THEN
480	ELSE
481	CALL report('main loop', 'Deallocating input buffer', cfg % debug)
482	DEALLOCATE(input_buffer)
483	END IF
484	CALL run_control('time', 'alloc')
485
486	!
487	!-- time steps / input files loop
488	END DO
489
490	IF (ALLOCATED(input_buffer)) THEN
491	CALL report('main loop', 'Deallocating input buffer', cfg % debug)
492	DEALLOCATE(input_buffer)
493	END IF
494	CALL run_control('time', 'alloc')
495
496	ELSE
497
498	message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'"
499	IF ( ALLOCATED(group % in_var_list) ) THEN
500	message = TRIM(message) // " with input variable '" // &
501	TRIM(group % in_var_list(1) % name) // "'."
502	END IF
503
504	CALL report('main loop', message, cfg % debug)
505
506	!
507	!-- IO group % to_be_processed conditional
508	END IF
509
510	!
511	!-- IO groups loop
512	END DO
513
514	!------------------------------------------------------------------------------
515	!- Section 3: Clean up.
516	!------------------------------------------------------------------------------
517	CALL fini_file_lists()
518	CALL fini_io_groups()
519	CALL fini_variables()
520	!CALL fini_grids()
521	CALL run_control('time', 'alloc')
522	CALL run_control('report', 'void')
523
524	message = "Finished writing dynamic driver '" // TRIM(output_file % name) // &
525	"' successfully."
526	CALL report('main loop', message)
527
528
529	END PROGRAM inifor

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |