Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |