Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |