Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |