Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |