Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/inifor_io.f90 @ 4659

Last change on this file since 4659 was 4659, checked in by eckhard, 4 years ago
inifor: Support for COSMO cloud water and precipitation
Property svn:keywords set to `Id`
File size: 63.7 KB

Line
1	!> @file src/inifor_io.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_io.f90 4659 2020-08-31 11:21:17Z eckhard $
28	! New command-line option '--precipitation' enables output of precipitation
29	! surface forcing
30	! Only define netCDF variables in enabled IO groups
31	! Improved code formatting
32	!
33	!
34	! 4569 2020-06-19 12:46:17Z eckhard
35	! Removed unused variable
36	!
37	!
38	! 4568 2020-06-19 11:56:30Z eckhard
39	! Handle COSMO soil data with and without additional surface temperature
40	!
41	!
42	! 4553 2020-06-03 16:34:15Z eckhard
43	! Option --help now points the user to INIFOR's wiki page
44	! Remove deprecated command-line options -clon and -clat
45	!
46	!
47	! 4538 2020-05-18 13:45:35Z eckhard
48	! Make setting the vertical PALM origin mandatory
49	!
50	!
51	! 4523 2020-05-07 15:58:16Z eckhard
52	! respect integer working precision (iwp) specified in inifor_defs.f90
53	!
54	!
55	! 4481 2020-03-31 18:55:54Z maronga
56	! Pass hhl_file directly instead of entire INIFOR configuration
57	!
58	!
59	! 3997 2019-05-23 12:35:57Z eckhard
60	! Added boolean indicator for --elevation option invocation
61	! Stop INIFOR if no command-line options given
62	!
63	!
64	! 3866 2019-04-05 14:25:01Z eckhard
65	! Use PALM's working precision
66	! Improved coding style
67	!
68	!
69	! 3801 2019-03-15 17:14:25Z eckhard
70	! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
71	! Moved get_soil_layer_thickness() here from inifor_grid
72	!
73	! 3785 2019-03-06 10:41:14Z eckhard
74	! Temporariliy disabled height-based geostrophic wind averaging
75	! Improved variable naming
76	!
77	!
78	! 3764 2019-02-26 13:42:09Z eckhard
79	! Removed dependency on radiation input files
80	!
81	!
82	! 3716 2019-02-05 17:02:38Z eckhard
83	! Removed dependency on soilmoisture input files
84	!
85	!
86	! 3680 2019-01-18 14:54:12Z knoop
87	! Moved get_input_file_list() here from grid module, added check for presence of
88	! input files
89	!
90	!
91	! 3618 2018-12-10 13:25:22Z eckhard
92	! Prefixed all INIFOR modules with inifor_
93	!
94	!
95	! 3615 2018-12-10 07:21:03Z raasch
96	! bugfix: abort replaced by inifor_abort
97	!
98	! 3557 2018-11-22 16:01:22Z eckhard
99	! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
100	!
101	!
102	! 3537 2018-11-20 10:53:14Z eckhard
103	! New routine get_netcdf_dim_vector()
104	!
105	!
106	! 3534 2018-11-19 15:35:16Z raasch
107	! bugfix: INTENT attribute changed
108	!
109	! 3456 2018-10-30 14:29:54Z eckhard
110	! NetCDf output of internal arrays only with --debug option
111	!
112	!
113	! 3447 2018-10-29 15:52:54Z eckhard
114	! Removed INCLUDE statement for get_netcdf_variable()
115	! Renamed source files for compatibilty with PALM build system
116	!
117	!
118	! 3395 2018-10-22 17:32:49Z eckhard
119	! Added command-line options for configuring the computation of geostrophic
120	! winds (--averagin-mode, --averaging-angle)
121	! Added command-line option --input-prefix for setting input file prefixes all
122	! at once
123	! Added --debug option for more verbose terminal output
124	! Added option-specific *_is_set LOGICALs to indicate invocation from the
125	! command-line
126	! Improved error messages in case of empty file-name strings
127	! Improved routine naming
128	!
129	! 3183 2018-07-27 14:25:55Z suehring
130	! Introduced new PALM grid stretching
131	! Updated variable names and metadata for PIDS v1.9 compatibility
132	! Improved handling of the start date string
133	! Better compatibility with older Intel compilers:
134	! - avoiding implicit array allocation with new get_netcdf_variable()
135	! subroutine instead of function
136	! Improved command line interface:
137	! - Added configuration validation
138	! - New options to configure input file prefixes
139	! - GNU-style short and long option names
140	! - Added version and copyright output
141	!
142	!
143	! 3182 2018-07-27 13:36:03Z suehring
144	! Initial revision
145	!
146	!
147	!
148	! Authors:
149	! --------
150	!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
151	!
152	! Description:
153	! ------------
154	!> The io module contains the functions needed to read and write netCDF data in
155	!> INIFOR.
156	!------------------------------------------------------------------------------!
157	#if defined ( __netcdf )
158	MODULE inifor_io
159
160	USE inifor_control
161	USE inifor_defs, &
162	ONLY: DATE, SNAME, PATH, PI, TO_RADIANS, TO_DEGREES, VERSION, &
163	NC_DEPTH_DIM_IDX, NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME, &
164	NC_RLON_NAME, NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME, &
165	NC_POLE_LONGITUDE_NAME, RHO_L, iwp, wp, &
166	PIDS_ORIGIN_LON, PIDS_ORIGIN_LAT, PIDS_ORIGIN_Z
167	USE inifor_types
168	USE inifor_util, &
169	ONLY: add_hours_to, nearly_equal, reverse, str, real_to_str
170	USE netcdf
171
172	IMPLICIT NONE
173
174	!------------------------------------------------------------------------------!
175	! Description:
176	! ------------
177	!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
178	!> variable 'in_var%name' from the file 'in_file'. The netCDF data array is
179	!> stored in the 'buffer' array and metadata added to the respective members of
180	!> the given 'in_var'.
181	!------------------------------------------------------------------------------!
182	INTERFACE get_netcdf_variable
183	MODULE PROCEDURE get_netcdf_variable_int
184	MODULE PROCEDURE get_netcdf_variable_real
185	END INTERFACE get_netcdf_variable
186
187	PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real
188
189	CONTAINS
190
191	!------------------------------------------------------------------------------!
192	! Description:
193	! ------------
194	!> get_netcdf_variable_int() implements the integer variant for the
195	!> get_netcdf_variable interface.
196	!------------------------------------------------------------------------------!
197	SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
198
199	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
200	TYPE(nc_var), INTENT(INOUT) :: in_var
201	INTEGER(iwp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
202
203	INTEGER :: ncid
204	INTEGER, DIMENSION(3) :: start, count
205
206	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
207	nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
208
209	CALL get_input_dimensions(in_var, ncid)
210
211	CALL get_netcdf_start_and_count(in_var, start, count)
212	CALL log_runtime('time', 'read')
213
214	ALLOCATE( buffer( count(1), count(2), count(3) ) )
215	CALL log_runtime('time', 'alloc')
216
217	CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
218	start = start, &
219	count = count ))
220
221	ELSE
222
223	message = "Failed to read '" // TRIM(in_var%name) // &
224	"' from file '" // TRIM(in_file) // "'."
225	CALL inifor_abort('get_netcdf_variable', message)
226
227	ENDIF
228
229	CALL check(nf90_close(ncid))
230	CALL log_runtime('time', 'read')
231
232	END SUBROUTINE get_netcdf_variable_int
233
234
235	!------------------------------------------------------------------------------!
236	! Description:
237	! ------------
238	!> get_netcdf_variable_real() implements the real variant for the
239	!> get_netcdf_variable interface.
240	!------------------------------------------------------------------------------!
241	SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
242
243	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
244	TYPE(nc_var), INTENT(INOUT) :: in_var
245	REAL(wp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
246
247	INTEGER :: ncid
248	INTEGER, DIMENSION(3) :: start, count
249
250	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
251	nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
252
253	CALL get_input_dimensions(in_var, ncid)
254
255	CALL get_netcdf_start_and_count(in_var, start, count)
256	CALL log_runtime('time', 'read')
257
258	ALLOCATE( buffer( count(1), count(2), count(3) ) )
259	CALL log_runtime('time', 'alloc')
260
261	CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
262	start = start, &
263	count = count ))
264
265	ELSE
266
267	message = "Failed to read '" // TRIM(in_var%name) // &
268	"' from file '" // TRIM(in_file) // "'."
269	CALL inifor_abort('get_netcdf_variable', message)
270
271	ENDIF
272
273	CALL check(nf90_close(ncid))
274	CALL log_runtime('time', 'read')
275
276	END SUBROUTINE get_netcdf_variable_real
277
278
279	!------------------------------------------------------------------------------!
280	! Description:
281	! ------------
282	!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
283	!> netCDF file 'filename'.
284	!------------------------------------------------------------------------------!
285	SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
286
287	CHARACTER(LEN=*), INTENT(IN) :: filename
288	CHARACTER(LEN=*), INTENT(IN) :: coordname
289	REAL(wp), ALLOCATABLE, INTENT(INOUT) :: coords(:)
290
291	INTEGER :: ncid, varid, dimlen
292	INTEGER :: dimids(NF90_MAX_VAR_DIMS)
293
294	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
295	nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR ) THEN
296
297	CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
298	CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
299
300	ALLOCATE(coords(dimlen))
301	CALL check(nf90_get_var( ncid, varid, coords))
302
303	ELSE
304
305	message = "Failed to read '" // TRIM(coordname) // &
306	"' from file '" // TRIM(filename) // "'."
307	CALL inifor_abort('get_netcdf_dim_vector', message)
308
309	ENDIF
310
311	END SUBROUTINE get_netcdf_dim_vector
312
313
314	!------------------------------------------------------------------------------!
315	! Description:
316	! ------------
317	!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
318	!> by 'in_var%name' into 'in_var' data structure.
319	!------------------------------------------------------------------------------!
320	SUBROUTINE get_input_dimensions(in_var, ncid)
321
322	TYPE(nc_var), INTENT(INOUT) :: in_var
323	INTEGER, INTENT(IN) :: ncid
324
325	INTEGER :: i
326
327	CALL check(nf90_get_att( ncid, in_var%varid, "long_name", &
328	in_var%long_name))
329
330	CALL check(nf90_get_att( ncid, in_var%varid, "units", &
331	in_var%units ))
332
333	CALL check(nf90_inquire_variable( ncid, in_var%varid, &
334	ndims = in_var%ndim, &
335	dimids = in_var%dimids ))
336
337	DO i = 1, in_var%ndim
338	CALL check(nf90_inquire_dimension( ncid, in_var%dimids(i), &
339	name = in_var%dimname(i), &
340	len = in_var%dimlen(i) ))
341	ENDDO
342
343	END SUBROUTINE get_input_dimensions
344
345
346	!------------------------------------------------------------------------------!
347	! Description:
348	! ------------
349	!> get_netcdf_start_and_count() gets the start position and element counts for
350	!> the given netCDF file. This information is used in get_netcdf_variable_int()
351	!> and _real() for reading input variables..
352	!------------------------------------------------------------------------------!
353	SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
354
355	TYPE(nc_var), INTENT(INOUT) :: in_var
356	INTEGER, DIMENSION(3), INTENT(OUT) :: start, count
357
358	INTEGER :: ndim
359
360	IF ( in_var%ndim .LT. 2 .OR. in_var%ndim .GT. 4 ) THEN
361
362	message = "Failed reading NetCDF variable " // &
363	TRIM(in_var%name) // " with " // &
364	TRIM(str(INT(in_var%ndim, kind=iwp))) // &
365	" dimensions because only two- and and three-dimensional" // &
366	" variables are supported."
367	CALL inifor_abort('get_netcdf_start_and_count', message)
368
369	ENDIF
370
371	start = (/ 1, 1, 1 /)
372	IF ( TRIM(in_var%name) .EQ. 'T_SO' .AND. &
373	in_var%has_redundant_first_level ) THEN
374
375	!
376	!-- Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
377	in_var%dimlen(3) = in_var%dimlen(3) - 1
378
379	!
380	!-- Start reading from second level, e.g. depth = 0.005 instead of 0.0
381	start(3) = 2
382	ENDIF
383
384	IF (in_var%ndim .EQ. 2) THEN
385	in_var%dimlen(3) = 1
386	ENDIF
387
388	ndim = MIN(in_var%ndim, 3)
389	count = (/ 1, 1, 1 /)
390	count(1:ndim) = in_var%dimlen(1:ndim)
391
392	END SUBROUTINE get_netcdf_start_and_count
393
394
395	!------------------------------------------------------------------------------!
396	! Description:
397	! ------------
398	!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
399	!> output.
400	!------------------------------------------------------------------------------!
401	SUBROUTINE netcdf_define_variable(var, ncid)
402
403	TYPE(nc_var), INTENT(INOUT) :: var
404	INTEGER, INTENT(IN) :: ncid
405
406	CALL check(nf90_def_var(ncid, var%name, NF90_FLOAT, var%dimids(1:var%ndim), var%varid))
407	CALL check(nf90_put_att(ncid, var%varid, "long_name", var%long_name))
408	CALL check(nf90_put_att(ncid, var%varid, "units", var%units))
409	IF ( var%lod .GE. 0 ) THEN
410	CALL check(nf90_put_att(ncid, var%varid, "lod", var%lod))
411	ENDIF
412	CALL check(nf90_put_att(ncid, var%varid, "source", var%source))
413	CALL check(nf90_put_att(ncid, var%varid, "_FillValue", NF90_FILL_REAL))
414
415	END SUBROUTINE netcdf_define_variable
416
417
418	!------------------------------------------------------------------------------!
419	! Description:
420	! ------------
421	!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
422	!> them in the given the 'var' data structure. This information is used later
423	!> for writing output variables in update_output().
424	!------------------------------------------------------------------------------!
425	SUBROUTINE netcdf_get_dimensions(var, ncid)
426
427	TYPE(nc_var), INTENT(INOUT) :: var
428	INTEGER, INTENT(IN) :: ncid
429	INTEGER :: i
430	CHARACTER(SNAME) :: null
431
432	DO i = 1, var%ndim
433	CALL check(nf90_inquire_dimension(ncid, var%dimids(i), &
434	name = null, &
435	len = var%dimlen(i) ) )
436	ENDDO
437
438	END SUBROUTINE netcdf_get_dimensions
439
440
441	!------------------------------------------------------------------------------!
442	! Description:
443	! ------------
444	!> This routine parses and interpretes the command-line options and stores the
445	!> resulting settings in the 'cfg' data structure.
446	!------------------------------------------------------------------------------!
447	SUBROUTINE parse_command_line_arguments( cfg )
448
449	TYPE(inifor_config), INTENT(INOUT) :: cfg
450
451	CHARACTER(LEN=PATH) :: option, arg
452	INTEGER :: arg_count, i
453
454	cfg%flow_prefix_is_set = .FALSE.
455	cfg%input_prefix_is_set = .FALSE.
456	cfg%p0_is_set = .FALSE.
457	cfg%radiation_prefix_is_set = .FALSE.
458	cfg%soil_prefix_is_set = .FALSE.
459	cfg%precipitation_prefix_is_set = .FALSE.
460	cfg%process_precipitation = .FALSE.
461	cfg%static_driver_is_set = .FALSE.
462	cfg%ug_defined_by_user = .FALSE.
463	cfg%vg_defined_by_user = .FALSE.
464	cfg%z0_is_set = .FALSE.
465
466	arg_count = COMMAND_ARGUMENT_COUNT()
467	IF (arg_count .GT. 0) THEN
468
469	i = 1
470	DO WHILE (i .LE. arg_count)
471
472	CALL GET_COMMAND_ARGUMENT( i, option )
473
474	SELECT CASE( TRIM(option) )
475
476	CASE( '--averaging-mode' )
477	CALL get_option_argument( i, arg )
478	cfg%averaging_mode = TRIM(arg)
479
480	CASE( '-date', '-d', '--date' )
481	CALL get_option_argument( i, arg )
482	cfg%start_date = TRIM(arg)
483
484	CASE( '--debug' )
485	cfg%debug = .TRUE.
486
487	CASE( '-z0', '-z', '--elevation' )
488	cfg%z0_is_set = .TRUE.
489	CALL get_option_argument( i, arg )
490	READ(arg, *) cfg%z0
491
492	CASE( '-p0', '-r', '--surface-pressure' )
493	cfg%p0_is_set = .TRUE.
494	CALL get_option_argument( i, arg )
495	READ(arg, *) cfg%p0
496
497	CASE( '-ug', '-u', '--geostrophic-u' )
498	cfg%ug_defined_by_user = .TRUE.
499	CALL get_option_argument( i, arg )
500	READ(arg, *) cfg%ug
501
502	CASE( '-vg', '-v', '--geostrophic-v' )
503	cfg%vg_defined_by_user = .TRUE.
504	CALL get_option_argument( i, arg )
505	READ(arg, *) cfg%vg
506
507	CASE( '-path', '-p', '--path' )
508	CALL get_option_argument( i, arg )
509	cfg%input_path = TRIM(arg)
510
511	CASE( '-hhl', '-l', '--hhl-file' )
512	CALL get_option_argument( i, arg )
513	cfg%hhl_file = TRIM(arg)
514
515	CASE( '--input-prefix')
516	CALL get_option_argument( i, arg )
517	cfg%input_prefix = TRIM(arg)
518	cfg%input_prefix_is_set = .TRUE.
519
520	CASE( '-a', '--averaging-angle' )
521	CALL get_option_argument( i, arg )
522	READ(arg, *) cfg%averaging_angle
523
524	CASE( '-static', '-t', '--static-driver' )
525	cfg%static_driver_is_set = .TRUE.
526	CALL get_option_argument( i, arg )
527	cfg%static_driver_file = TRIM(arg)
528
529	CASE( '-soil', '-s', '--soil-file')
530	CALL get_option_argument( i, arg )
531	cfg%soiltyp_file = TRIM(arg)
532
533	CASE( '--flow-prefix')
534	CALL get_option_argument( i, arg )
535	cfg%flow_prefix = TRIM(arg)
536	cfg%flow_prefix_is_set = .TRUE.
537
538	CASE( '--radiation-prefix')
539	CALL get_option_argument( i, arg )
540	cfg%radiation_prefix = TRIM(arg)
541	cfg%radiation_prefix_is_set = .TRUE.
542
543	CASE( '--soil-prefix')
544	CALL get_option_argument( i, arg )
545	cfg%soil_prefix = TRIM(arg)
546	cfg%soil_prefix_is_set = .TRUE.
547
548	CASE( '--precipitation-prefix')
549	CALL get_option_argument( i, arg )
550	cfg%precipitation_prefix = TRIM(arg)
551	cfg%precipitation_prefix_is_set = .TRUE.
552
553	CASE( '--precipitation')
554	cfg%process_precipitation = .TRUE.
555
556	CASE( '-o', '--output' )
557	CALL get_option_argument( i, arg )
558	cfg%output_file = TRIM(arg)
559
560	CASE( '-n', '--namelist' )
561	CALL get_option_argument( i, arg )
562	cfg%namelist_file = TRIM(arg)
563
564	CASE( '-mode', '-i', '--init-mode' )
565	CALL get_option_argument( i, arg )
566	cfg%ic_mode = TRIM(arg)
567
568	CASE( '-f', '--forcing-mode' )
569	CALL get_option_argument( i, arg )
570	cfg%bc_mode = TRIM(arg)
571
572	CASE( '--version' )
573	CALL print_version
574	STOP
575
576	CASE( '--help' )
577	CALL print_version
578	PRINT *, ""
579	PRINT *, &
580	"For documentation and a list of available command-line options " // NEW_LINE( " " ) // &
581	" please visit https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/inifor."
582	STOP
583
584	CASE DEFAULT
585	message = "unknown option '" // TRIM(option) // "'."
586	CALL inifor_abort('parse_command_line_arguments', message)
587
588	END SELECT
589
590	i = i + 1
591
592	ENDDO
593
594	ELSE
595
596	CALL print_version
597	CALL report( 'parse_command_line_arguments', 'No arguments present, exiting.' )
598	STOP
599
600	ENDIF
601
602	END SUBROUTINE parse_command_line_arguments
603
604
605
606	SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
607	step_hour, input_path, prefix, suffix, &
608	file_list )
609
610	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
611	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
612	INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
613	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
614
615	INTEGER(iwp) :: number_of_intervals, hour, i
616	CHARACTER(LEN=DATE) :: date_string
617
618	number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
619	ALLOCATE( file_list(number_of_intervals + 1) )
620
621	DO i = 0, number_of_intervals
622
623	hour = start_hour + i * step_hour
624	date_string = add_hours_to(start_date_string, hour)
625
626	file_list(i+1) = TRIM(input_path) // TRIM(prefix) // &
627	TRIM(date_string) // TRIM(suffix) // '.nc'
628
629	ENDDO
630
631	END SUBROUTINE get_datetime_file_list
632
633	!------------------------------------------------------------------------------!
634	! Description:
635	! ------------
636	!> Establish a list of files based on the given start and end times and file
637	!> prefixes and suffixes.
638	!------------------------------------------------------------------------------!
639	SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour, &
640	step_hour, input_path, prefix, suffix, &
641	file_list )
642
643	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
644	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
645	INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
646	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
647
648	CALL get_datetime_file_list( start_date_string, start_hour, end_hour, &
649	step_hour, input_path, prefix, suffix, &
650	file_list )
651
652	END SUBROUTINE get_input_file_list
653
654
655	!------------------------------------------------------------------------------!
656	! Description:
657	! ------------
658	!> Abort INIFOR if the given file is not present.
659	!------------------------------------------------------------------------------!
660	LOGICAL FUNCTION file_is_present( filename, file_kind, message )
661
662	CHARACTER(LEN=*), INTENT(IN) :: filename, file_kind
663	CHARACTER(LEN=*), INTENT(OUT) :: message
664
665	file_is_present = file_present(filename)
666
667	IF (.NOT. file_is_present) THEN
668
669	IF (LEN( TRIM( filename ) ) == 0) THEN
670	message = "No name was given for the " // TRIM( file_kind ) // " file."
671	ELSE
672	message = "The " // TRIM( file_kind ) // " file '" // &
673	TRIM( filename ) // "' was not found."
674	ENDIF
675
676	ENDIF
677
678	END FUNCTION file_is_present
679
680
681	!------------------------------------------------------------------------------!
682	! Description:
683	! ------------
684	!> Abort INIFOR if the given file is not present.
685	!------------------------------------------------------------------------------!
686	SUBROUTINE verify_file( file_name, file_kind )
687
688	CHARACTER(LEN=*), INTENT(IN) :: file_name, file_kind
689
690	IF (.NOT. file_is_present( file_name, file_kind, message )) THEN
691
692	CALL inifor_abort( 'verify_file', message )
693
694	ENDIF
695
696	message = "Set up input file name '" // TRIM( file_name ) // "'"
697	CALL report( 'verify_file', message )
698
699	END SUBROUTINE verify_file
700	!------------------------------------------------------------------------------!
701	! Description:
702	! ------------
703	!> Get the argument of the i'th command line option, which is at the location
704	!> i+1 of the argument list.
705	!------------------------------------------------------------------------------!
706	SUBROUTINE get_option_argument( i, arg )
707	CHARACTER(LEN=PATH), INTENT(INOUT) :: arg
708	INTEGER, INTENT(INOUT) :: i
709
710	i = i + 1
711	CALL GET_COMMAND_ARGUMENT( i, arg )
712
713	END SUBROUTINE
714
715
716	!------------------------------------------------------------------------------!
717	! Description:
718	! ------------
719	!> Checks the INIFOR configuration 'cfg' for plausibility.
720	!------------------------------------------------------------------------------!
721	SUBROUTINE validate_config( cfg )
722	TYPE(inifor_config), INTENT(IN) :: cfg
723
724	CALL verify_file( cfg%hhl_file, 'HHL' )
725	CALL verify_file( cfg%namelist_file, 'NAMELIST' )
726	CALL verify_file( cfg%soiltyp_file, 'SOILTYP' )
727
728	!
729	!-- Only check optional static driver file name, if it has been given.
730	IF (TRIM( cfg%static_driver_file ) .NE. '') THEN
731	CALL verify_file( cfg%static_driver_file, 'static driver' )
732	ENDIF
733
734	SELECT CASE( TRIM( cfg%ic_mode ) )
735	CASE( 'profile', 'volume' )
736	CASE DEFAULT
737	message = "Initialization mode '" // TRIM( cfg%ic_mode ) //&
738	"' is not supported. " //&
739	"Please select either 'profile' or 'volume', " //&
740	"or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
741	"to the latter."
742	CALL inifor_abort( 'validate_config', message )
743	END SELECT
744
745	SELECT CASE( TRIM(cfg%bc_mode) )
746	CASE( 'real', 'ideal')
747	CASE DEFAULT
748	message = "Forcing mode '" // TRIM( cfg%bc_mode ) //&
749	"' is not supported. " //&
750	"Please select either 'real' or 'ideal', " //&
751	"or omit the -f/--forcing-mode option entirely, which corresponds "//&
752	"to the latter."
753	CALL inifor_abort( 'validate_config', message )
754	END SELECT
755
756	SELECT CASE( TRIM( cfg%averaging_mode ) )
757	CASE( 'level' )
758	CASE( 'height' )
759	message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
760	"' is currently not supported. " //&
761	"Please use level-based averaging by selecting 'level', " //&
762	"or by omitting the --averaging-mode option entirely."
763	CALL inifor_abort( 'validate_config', message )
764	CASE DEFAULT
765	message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
766	"' is not supported. " //&
767	! "Please select either 'height' or 'level', " //&
768	! "or omit the --averaging-mode option entirely, which corresponds "//&
769	! "to the latter."
770	"Please use level-based averaging by selecting 'level', " //&
771	"or by omitting the --averaging-mode option entirely."
772	CALL inifor_abort( 'validate_config', message )
773	END SELECT
774
775	IF ( cfg%ug_defined_by_user .NEQV. cfg%vg_defined_by_user ) THEN
776	message = "You specified only one component of the geostrophic " // &
777	"wind. Please specify either both or none."
778	CALL inifor_abort( 'validate_config', message )
779	ENDIF
780
781	IF ( .NOT. cfg%static_driver_is_set .AND. .NOT. cfg%z0_is_set ) THEN
782	message = &
783	"The vertical origin of the PALM grid has not been defined. " // NEW_LINE( " " ) // &
784	"Please specify the right value for your setup by either " // NEW_LINE( " " ) // &
785	" - using the command-line option --elevation <height above sea level>, or by" // NEW_LINE( " " ) // &
786	" - specifying a static driver file using --static <filename> in order to use " // NEW_LINE( " " ) // &
787	" use the value of origin_z (and origin_lon and origin_lat) specifed therein."
788	CALL inifor_abort( 'validate_config', message )
789	ENDIF
790
791	END SUBROUTINE validate_config
792
793
794	SUBROUTINE get_cosmo_grid( hhl_file, soil_file, rlon, rlat, hhl, hfl, depths, &
795	d_depth, d_depth_rho_inv, phi_n, lambda_n, &
796	phi_equat, &
797	lonmin_cosmo, lonmax_cosmo, &
798	latmin_cosmo, latmax_cosmo, &
799	nlon, nlat, nlev, ndepths )
800
801	CHARACTER(LEN=PATH), INTENT(IN) :: hhl_file !< path to file containing the HHL variable (height of half layers)
802	CHARACTER(LEN=PATH), INTENT(IN) :: soil_file !< path to one of the soil input files for reading soil layer depths
803	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlon !< longitudes of COSMO-DE's rotated-pole grid
804	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlat !< latitudes of COSMO-DE's rotated-pole grid
805	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
806	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
807	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: depths !< COSMO-DE's TERRA-ML soil layer depths
808	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth
809	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth_rho_inv
810	REAL(wp), INTENT(OUT) :: phi_n
811	REAL(wp), INTENT(OUT) :: phi_equat
812	REAL(wp), INTENT(OUT) :: lambda_n
813	REAL(wp), INTENT(OUT) :: lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
814	REAL(wp), INTENT(OUT) :: lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
815	REAL(wp), INTENT(OUT) :: latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
816	REAL(wp), INTENT(OUT) :: latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
817	INTEGER(iwp), INTENT(OUT) :: nlon, nlat, nlev, ndepths
818
819	TYPE(nc_var) :: cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
820	INTEGER(iwp) :: k
821
822	!
823	!-- Read in COSMO's heights of half layers (vertical cell faces)
824	cosmo_var%name = NC_HHL_NAME
825	CALL get_netcdf_variable( hhl_file, cosmo_var, hhl )
826	CALL get_netcdf_dim_vector( hhl_file, NC_RLON_NAME, rlon )
827	CALL get_netcdf_dim_vector( hhl_file, NC_RLAT_NAME, rlat )
828	CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
829	CALL log_runtime( 'time', 'read' )
830
831	CALL reverse( hhl )
832	nlon = SIZE( hhl, 1 )
833	nlat = SIZE( hhl, 2 )
834	nlev = SIZE( hhl, 3 )
835	ndepths = SIZE( depths )
836
837	CALL log_runtime( 'time', 'comp' )
838
839	ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
840	ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
841	CALL log_runtime('time', 'alloc')
842
843	CALL get_soil_layer_thickness( depths, d_depth )
844	d_depth_rho_inv = 1.0_wp / ( d_depth * RHO_L )
845
846	!
847	!-- Compute COSMO's heights of full layers (cell centres)
848	DO k = 1, nlev-1
849	hfl(:,:,k) = 0.5_wp * ( hhl(:,:,k) + &
850	hhl(:,:,k+1) )
851	ENDDO
852	!
853	!-- COSMO rotated pole coordinates
854	phi_n = TO_RADIANS &
855	* get_netcdf_variable_attribute( hhl_file, &
856	NC_ROTATED_POLE_NAME, &
857	NC_POLE_LATITUDE_NAME )
858
859	lambda_n = TO_RADIANS &
860	* get_netcdf_variable_attribute( hhl_file, &
861	NC_ROTATED_POLE_NAME, &
862	NC_POLE_LONGITUDE_NAME )
863
864	phi_equat = 90.0_wp * TO_RADIANS - phi_n
865
866	lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
867	lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
868	latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
869	latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
870	CALL log_runtime('time', 'comp')
871
872	END SUBROUTINE get_cosmo_grid
873
874
875	!------------------------------------------------------------------------------!
876	! Description:
877	! ------------
878	!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
879	!> TERRA_ML's [1]) soil layer boundaries follow the rule
880	!>
881	!> depth(0) = 0.0, and
882	!> depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
883	!>
884	!> and full levels are defined as the midpoints between two layer boundaries,
885	!> all except the first layer thicknesses equal the depth of the midpoint.
886	!>
887	!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
888	!> Physical Parameterization*, Sect. 11 TERRA_ML.
889	!> http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
890	!>
891	!> Input parameters:
892	!> -----------------
893	!>
894	!> depths: array of full soil layer depths (cell midpoints)
895	!>
896	!>
897	!> Output parameters:
898	!> ------------------
899	!>
900	!> d_depth: array of soil layer thicknesses
901	!>
902	!------------------------------------------------------------------------------!
903	SUBROUTINE get_soil_layer_thickness( depths, d_depth )
904
905	REAL(wp), INTENT(IN) :: depths(:)
906	REAL(wp), INTENT(OUT) :: d_depth(:)
907
908	d_depth(:) = depths(:)
909	d_depth(1) = 2.0_wp * depths(1)
910
911	END SUBROUTINE get_soil_layer_thickness
912	!------------------------------------------------------------------------------!
913	! Description:
914	! ------------
915	!> Check whether the given file is present on the filesystem.
916	!------------------------------------------------------------------------------!
917	LOGICAL FUNCTION file_present( filename )
918	CHARACTER(LEN=PATH), INTENT(IN) :: filename
919
920	INQUIRE( FILE=filename, EXIST=file_present )
921
922	END FUNCTION file_present
923
924
925	!------------------------------------------------------------------------------!
926	! Description:
927	! ------------
928	!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
929	!> file.
930	!>
931	!> Besides writing metadata, such as global attributes, coordinates, variables,
932	!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
933	!> writes the actual data.
934	!------------------------------------------------------------------------------!
935	SUBROUTINE setup_netcdf_dimensions( output_file, palm_grid, &
936	start_date_string, origin_lon, origin_lat )
937
938	TYPE(nc_file), INTENT(INOUT) :: output_file
939	TYPE(grid_definition), INTENT(IN) :: palm_grid
940	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
941	REAL(wp), INTENT(IN) :: origin_lon, origin_lat
942
943	CHARACTER (LEN=8) :: date_string
944	CHARACTER (LEN=10) :: time_string
945	CHARACTER (LEN=5) :: zone_string
946	CHARACTER (LEN=SNAME) :: history_string
947	INTEGER :: nx, ny, nz, nt
948	INTEGER :: ncid, dimids(3), dimvarids(3)
949	REAL(wp) :: z0
950
951	message = "Initializing PALM-4U dynamic driver file '" // &
952	TRIM(output_file%name) // "' and setting up dimensions."
953	CALL report('setup_netcdf_dimensions', message)
954
955	!
956	!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
957	#if defined( __netcdf4 )
958	CALL check(nf90_create(TRIM(output_file%name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
959	#else
960	CALL check(nf90_create(TRIM(output_file%name), NF90_CLOBBER, ncid))
961	#endif
962
963	!------------------------------------------------------------------------------
964	!- Section 1: Define NetCDF dimensions and coordinates
965	!------------------------------------------------------------------------------
966	nt = SIZE(output_file%time)
967	nx = palm_grid%nx
968	ny = palm_grid%ny
969	nz = palm_grid%nz
970	z0 = palm_grid%z0
971
972
973	!
974	!------------------------------------------------------------------------------
975	!- Section 2: Write global NetCDF attributes
976	!------------------------------------------------------------------------------
977	CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
978	history_string = &
979	'Created on '// date_string // &
980	' at ' // time_string(1:2) // ':' // time_string(3:4) // &
981	' (UTC' // zone_string // ')'
982
983	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...'))
984	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach'))
985	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
986	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string)))
987	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--'))
988	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--'))
989	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
990	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
991	!
992	!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
993	!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
994	!-- FIXME: Standard v1.9., origin_z)
995	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)'))))
996	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
997	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--'))
998
999	!
1000	!
1001	!------------------------------------------------------------------------------
1002	!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
1003	!------------------------------------------------------------------------------
1004	!
1005	!-- reset dimids first
1006	dimids = (/0, 0, 0/)
1007	CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
1008	CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
1009	CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
1010	!
1011	!-- save dimids for later
1012	output_file%dimids_scl = dimids
1013
1014	!
1015	!-- reset dimvarids first
1016	dimvarids = (/0, 0, 0/)
1017	CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
1018	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
1019	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1020
1021	CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
1022	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
1023	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1024
1025	CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
1026	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
1027	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1028	!
1029	!-- save dimvarids for later
1030	output_file%dimvarids_scl = dimvarids
1031
1032	!
1033	!-- overwrite third dimid with the one of depth
1034	CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid%depths), dimids(3)) )
1035	!
1036	!-- save dimids for later
1037	output_file%dimids_soil = dimids
1038
1039	!
1040	!-- overwrite third dimvarid with the one of depth
1041	CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file%dimids_soil(3), dimvarids(3)))
1042	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
1043	CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
1044	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1045	!
1046	!-- save dimvarids for later
1047	output_file%dimvarids_soil = dimvarids
1048	!
1049	!------------------------------------------------------------------------------
1050	!- Section 2b: Define dimensions for cell faces/velocities
1051	!------------------------------------------------------------------------------
1052	!
1053	!-- reset dimids first
1054	dimids = (/0, 0, 0/)
1055	CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
1056	CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
1057	CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
1058	!
1059	!-- save dimids for later
1060	output_file%dimids_vel = dimids
1061
1062	!
1063	!-- reset dimvarids first
1064	dimvarids = (/0, 0, 0/)
1065	CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
1066	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
1067	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1068
1069	CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
1070	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
1071	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1072
1073	CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
1074	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
1075	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1076	!
1077	!-- save dimvarids for later
1078	output_file%dimvarids_vel = dimvarids
1079
1080	!
1081	!------------------------------------------------------------------------------
1082	!- Section 2c: Define time dimension
1083	!------------------------------------------------------------------------------
1084	CALL check(nf90_def_dim(ncid, "time", nt, output_file%dimid_time) )
1085	CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
1086	output_file%dimid_time, &
1087	output_file%dimvarid_time))
1088	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "standard_name", "time"))
1089	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "long_name", "time"))
1090	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "units", &
1091	"seconds since " // start_date_string // " UTC"))
1092
1093	CALL check(nf90_enddef(ncid))
1094
1095	!
1096	!------------------------------------------------------------------------------
1097	!- Section 3: Write grid coordinates
1098	!------------------------------------------------------------------------------
1099	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(1), palm_grid%x))
1100	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(2), palm_grid%y))
1101	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(3), palm_grid%z))
1102
1103	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(1), palm_grid%xu))
1104	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(2), palm_grid%yv))
1105	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(3), palm_grid%zw))
1106
1107	!
1108	!-- TODO Read in soil depths from input file before this.
1109	CALL check(nf90_put_var(ncid, output_file%dimvarids_soil(3), palm_grid%depths))
1110
1111	!
1112	!-- Write time vector
1113	CALL check(nf90_put_var(ncid, output_file%dimvarid_time, output_file%time))
1114
1115	!
1116	!-- Close the file
1117	CALL check(nf90_close(ncid))
1118
1119	END SUBROUTINE setup_netcdf_dimensions
1120
1121
1122	!------------------------------------------------------------------------------!
1123	! Description:
1124	! ------------
1125	!> Defines the netCDF variables to be written to the dynamic driver file
1126	!------------------------------------------------------------------------------!
1127	SUBROUTINE setup_netcdf_variables(filename, io_group_list)
1128
1129	CHARACTER (LEN=*), INTENT(IN) :: filename
1130	TYPE(io_group), INTENT(IN), TARGET :: io_group_list(:)
1131
1132	TYPE(io_group), POINTER :: group
1133	TYPE(nc_var), POINTER :: var
1134	INTEGER(iwp) :: group_idx, var_idx, n_var
1135	INTEGER :: ncid
1136	LOGICAL :: to_be_written
1137
1138	message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
1139	CALL report('setup_netcdf_variables', message)
1140
1141	CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
1142	CALL check(nf90_redef(ncid))
1143
1144	n_var = 0
1145	DO group_idx = 1, SIZE( io_group_list )
1146
1147	group => io_group_list(group_idx)
1148	DO var_idx = 1, SIZE( group%out_vars )
1149
1150	to_be_written = .FALSE.
1151
1152	IF (ALLOCATED( group%out_vars )) THEN
1153	var => group%out_vars(var_idx)
1154	n_var = n_var + 1
1155
1156	to_be_written = ( &
1157	group%to_be_processed .AND. var%to_be_processed &
1158	.AND. .NOT. var%is_internal &
1159	)
1160	ENDIF
1161
1162	IF ( to_be_written ) THEN
1163	message = " variable #" // TRIM(str(n_var)) // " '" // TRIM(var%name) // "'."
1164	CALL report('setup_netcdf_variables', message)
1165
1166	CALL netcdf_define_variable(var, ncid)
1167	CALL netcdf_get_dimensions(var, ncid)
1168	ENDIF
1169
1170	ENDDO
1171
1172	ENDDO
1173
1174	CALL check(nf90_enddef(ncid))
1175	CALL check(nf90_close(ncid))
1176
1177	message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
1178	CALL report('setup_netcdf_variables', message)
1179
1180	END SUBROUTINE setup_netcdf_variables
1181
1182
1183	!------------------------------------------------------------------------------!
1184	! Description:
1185	! ------------
1186	!> This routine reads and returns all input variables of the given IO group
1187	!> It accomodates the data by allocating a container variable that represents a
1188	!> list of arrays of the same length as the groups variable list. (This list
1189	!> will typically contain one or two items.) After the container, its members
1190	!> are allocated one by one with the appropriate, possibly different,
1191	!> dimensions.
1192	!>
1193	!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
1194	!> record netCDF IDs in the 'in_var_list()' member variable.
1195	!------------------------------------------------------------------------------!
1196	SUBROUTINE read_input_variables(group, iter, buffer)
1197	TYPE(io_group), INTENT(INOUT), TARGET :: group
1198	INTEGER(iwp), INTENT(IN) :: iter
1199	TYPE(container), ALLOCATABLE, INTENT(INOUT) :: buffer(:)
1200	INTEGER(iwp) :: hour, buf_id
1201	TYPE(nc_var), POINTER :: input_var
1202	CHARACTER(LEN=PATH), POINTER :: input_file
1203	INTEGER(iwp) :: ivar, nbuffers
1204
1205	message = "Reading data for I/O group '" // TRIM(group%in_var_list(1)%name) // "'."
1206	CALL report('read_input_variables', message)
1207
1208	input_file => group%in_files(iter)
1209
1210	!
1211	!------------------------------------------------------------------------------
1212	!- Section 1: Load input buffers for accumulated variables
1213	!------------------------------------------------------------------------------
1214	!
1215	!-- radiation budgets, precipitation
1216	IF ( group%kind == 'running average' .OR. &
1217	group%kind == 'accumulated' ) THEN
1218
1219	IF ( SIZE( group%in_var_list ) .GT. 1 ) THEN
1220	message = "I/O groups may not contain more than one " // &
1221	"accumulated variable. Group '" // TRIM(group%kind) //&
1222	"' contains " // &
1223	TRIM( str( SIZE( group%in_var_list, kind=iwp ) ) ) // "."
1224	CALL inifor_abort( 'read_input_variables \| accumulation', message )
1225	ENDIF
1226
1227	!
1228	!-- use two buffer arrays
1229	nbuffers = 2
1230	IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) )
1231
1232	!
1233	!-- hour of the day
1234	hour = iter - 1
1235	!
1236	!-- chose correct buffer array
1237	buf_id = select_buffer(hour)
1238
1239	CALL log_runtime('time', 'read')
1240	IF ( ALLOCATED(buffer(buf_id)%array) ) DEALLOCATE(buffer(buf_id)%array)
1241	CALL log_runtime('time', 'alloc')
1242
1243	input_var => group%in_var_list(1)
1244	CALL get_netcdf_variable(input_file, input_var, buffer(buf_id)%array)
1245	CALL report('read_input_variables', "Read accumulated " // TRIM(group%in_var_list(1)%name))
1246
1247	IF ( input_var%is_upside_down ) CALL reverse(buffer(buf_id)%array)
1248	CALL log_runtime('time', 'comp')
1249
1250	!------------------------------------------------------------------------------
1251	!- Section 2: Load input buffers for normal I/O groups
1252	!------------------------------------------------------------------------------
1253	ELSE
1254
1255	!
1256	!-- Allocate one input buffer per output quantity. If more quantities
1257	!-- have to be computed than input variables exist in this group,
1258	!-- allocate more buffers. For instance, in the thermodynamics group,
1259	!-- there are three input variabels (PP, T, Qv) and four quantities
1260	!-- necessart (P, Theta, Rho, qv) for the corresponding output fields
1261	!-- (p0, Theta, qv, ug, and vg)
1262	ALLOCATE( buffer(group%n_output_quantities) )
1263	CALL log_runtime('time', 'alloc')
1264
1265	!
1266	!-- Read in all input variables, leave extra buffers-if any-untouched.
1267	DO ivar = 1, group%n_inputs
1268
1269	input_var => group%in_var_list(ivar)
1270
1271	IF ( input_var%to_be_processed ) THEN
1272	! Check wheather P or PP is present in input file
1273	IF (input_var%name == 'P') THEN
1274	input_var%name = TRIM( get_pressure_varname(input_file) )
1275	CALL log_runtime('time', 'read')
1276	ENDIF
1277
1278	CALL get_netcdf_variable(input_file, input_var, buffer(ivar)%array)
1279
1280	IF ( input_var%is_upside_down ) CALL reverse(buffer(ivar)%array)
1281	CALL log_runtime('time', 'comp')
1282	ENDIF
1283
1284	ENDDO
1285	ENDIF
1286
1287	END SUBROUTINE read_input_variables
1288
1289
1290	!------------------------------------------------------------------------------!
1291	! Description:
1292	! ------------
1293	!> Select the appropriate buffer ID for accumulated COSMO input variables
1294	!> depending on the current hour.
1295	!------------------------------------------------------------------------------!
1296	INTEGER(iwp) FUNCTION select_buffer(hour)
1297	INTEGER(iwp), INTENT(IN) :: hour
1298	INTEGER(iwp) :: step
1299
1300	select_buffer = 0_iwp
1301	step = MODULO(hour, 3_iwp) + 1_iwp
1302
1303	SELECT CASE(step)
1304	CASE(1, 3)
1305	select_buffer = 1_iwp
1306	CASE(2)
1307	select_buffer = 2_iwp
1308	CASE DEFAULT
1309	message = "Invalid step '" // TRIM(str(step))
1310	CALL inifor_abort('select_buffer', message)
1311	END SELECT
1312	END FUNCTION select_buffer
1313
1314
1315	!------------------------------------------------------------------------------!
1316	! Description:
1317	! ------------
1318	!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
1319	!> perturbation, 'PP', and returns the appropriate string.
1320	!------------------------------------------------------------------------------!
1321	CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
1322	CHARACTER(LEN=*) :: input_file
1323	INTEGER :: ncid, varid
1324
1325	CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
1326	IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR ) THEN
1327
1328	var = 'P'
1329
1330	ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR ) THEN
1331
1332	var = 'PP'
1333	CALL report('get_pressure_var', 'Using PP instead of P')
1334
1335	ELSE
1336
1337	message = "Failed to read '" // TRIM(var) // &
1338	"' from file '" // TRIM(input_file) // "'."
1339	CALL inifor_abort('get_pressure_var', message)
1340
1341	ENDIF
1342
1343	CALL check(nf90_close(ncid))
1344
1345	END FUNCTION get_pressure_varname
1346
1347
1348	!------------------------------------------------------------------------------!
1349	! Description:
1350	! ------------
1351	!> Read the given global attribute form the given netCDF file.
1352	!------------------------------------------------------------------------------!
1353	FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
1354
1355	CHARACTER(LEN=*), INTENT(IN) :: filename, attribute
1356	REAL(wp) :: attribute_value
1357
1358	INTEGER :: ncid
1359
1360	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
1361
1362	CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
1363	CALL check(nf90_close(ncid))
1364
1365	ELSE
1366
1367	message = "Failed to read '" // TRIM(attribute) // &
1368	"' from file '" // TRIM(filename) // "'."
1369	CALL inifor_abort('get_netcdf_attribute', message)
1370
1371	ENDIF
1372
1373	END FUNCTION get_netcdf_attribute
1374
1375
1376	!------------------------------------------------------------------------------!
1377	! Description:
1378	! ------------
1379	!> Read the attribute of the given variable form the given netCDF file.
1380	!------------------------------------------------------------------------------!
1381	FUNCTION get_netcdf_variable_attribute(filename, varname, attribute) &
1382	RESULT(attribute_value)
1383
1384	CHARACTER(LEN=*), INTENT(IN) :: filename, varname, attribute
1385	REAL(wp) :: attribute_value
1386
1387	INTEGER :: ncid, varid
1388
1389	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
1390
1391	CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1392	CALL check( nf90_get_att( ncid, varid, TRIM( attribute ), &
1393	attribute_value ) )
1394	CALL check( nf90_close( ncid ) )
1395
1396	ELSE
1397
1398	message = "Failed to read '" // TRIM( varname ) // ":" // &
1399	TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
1400	CALL inifor_abort('get_netcdf_variable_attribute', message)
1401
1402	ENDIF
1403
1404	END FUNCTION get_netcdf_variable_attribute
1405
1406	!------------------------------------------------------------------------------!
1407	! Description:
1408	! ------------
1409	!> Updates the dynamic driver with the interpolated field of the current
1410	!> variable at the current time step.
1411	!------------------------------------------------------------------------------!
1412	SUBROUTINE update_output(var, array, iter, output_file, cfg)
1413	TYPE(nc_var), INTENT(IN) :: var
1414	REAL(wp), INTENT(IN) :: array(:,:,:)
1415	INTEGER(iwp), INTENT(IN) :: iter
1416	TYPE(nc_file), INTENT(IN) :: output_file
1417	TYPE(inifor_config) :: cfg
1418
1419	INTEGER :: ncid, ndim, start(4), count(4)
1420	LOGICAL :: var_is_time_dependent
1421
1422	var_is_time_dependent = ( &
1423	var%dimids( var%ndim ) == output_file%dimid_time &
1424	)
1425
1426	!
1427	!-- Skip time dimension for output
1428	ndim = var%ndim
1429	IF ( var_is_time_dependent ) ndim = var%ndim - 1
1430
1431	start(:) = (/1,1,1,1/)
1432	start(ndim+1) = iter
1433	count(1:ndim) = var%dimlen(1:ndim)
1434
1435	CALL check(nf90_open(output_file%name, NF90_WRITE, ncid))
1436
1437	!
1438	!-- Reduce dimension of output array according to variable kind
1439	SELECT CASE (TRIM(var%kind))
1440
1441	CASE ( 'init scalar profile', 'init u profile', 'init v profile', &
1442	'init w profile' )
1443
1444	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
1445
1446	CASE ( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
1447
1448	CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
1449
1450	CASE( 'left scalar', 'right scalar', 'left w', 'right w' )
1451
1452	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
1453	start=start(1:ndim+1), &
1454	count=count(1:ndim) ) )
1455
1456
1457	IF (.NOT. SIZE(array, 2) .EQ. var%dimlen(1)) THEN
1458	PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
1459	TRIM(var%name), " (", var%dimlen(1), &
1460	") does not match the dimension of the output array (", &
1461	SIZE(array, 2), ")."
1462	STOP
1463	ENDIF
1464
1465
1466	CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
1467
1468	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
1469	start=start(1:ndim+1), &
1470	count=count(1:ndim) ) )
1471
1472
1473	CASE( 'surface forcing', 'top scalar', 'top w' )
1474
1475	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
1476	start=start(1:ndim+1), &
1477	count=count(1:ndim) ) )
1478
1479	CASE ( 'left u', 'right u', 'left v', 'right v' )
1480
1481	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
1482	start=start(1:ndim+1), &
1483	count=count(1:ndim) ) )
1484
1485	CASE ( 'north u', 'south u', 'north v', 'south v' )
1486
1487	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
1488	start=start(1:ndim+1), &
1489	count=count(1:ndim) ) )
1490
1491	CASE ( 'top u', 'top v' )
1492
1493	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
1494	start=start(1:ndim+1), &
1495	count=count(1:ndim) ) )
1496
1497	CASE ( 'time series' )
1498
1499	CALL check(nf90_put_var( ncid, var%varid, array(1,1,1), &
1500	start=start(1:ndim+1) ) )
1501
1502	CASE ( 'constant scalar profile', 'geostrophic' )
1503
1504	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
1505	start=start(1:ndim+1), &
1506	count=count(1:ndim) ) )
1507
1508	CASE ( 'internal profile' )
1509
1510	IF ( cfg%debug ) THEN
1511	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
1512	start=start(1:ndim+1), &
1513	count=count(1:ndim) ) )
1514	ENDIF
1515
1516	CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
1517
1518	message = "Doing nothing in terms of writing large-scale forings."
1519	CALL report('update_output', message)
1520
1521	CASE DEFAULT
1522
1523	message = "Variable kind '" // TRIM(var%kind) // &
1524	"' not recognized."
1525	CALL inifor_abort('update_output', message)
1526
1527	END SELECT
1528
1529	CALL check(nf90_close(ncid))
1530
1531	END SUBROUTINE update_output
1532
1533
1534	!------------------------------------------------------------------------------!
1535	! Description:
1536	! ------------
1537	!> Checks the status of a netCDF API call and aborts if an error occured
1538	!------------------------------------------------------------------------------!
1539	SUBROUTINE check(status)
1540
1541	INTEGER, INTENT(IN) :: status
1542
1543	IF (status /= nf90_noerr) THEN
1544	message = "NetCDF API call failed with error: " // &
1545	TRIM( nf90_strerror(status) )
1546	CALL inifor_abort('io.check', message)
1547	ENDIF
1548
1549	END SUBROUTINE check
1550
1551
1552	!------------------------------------------------------------------------------!
1553	! Description:
1554	! ------------
1555	!> Setup the origin of the PALM coordinate system based on what is given in the
1556	!> INIFOR namelist file and via an optional static driver.
1557	!------------------------------------------------------------------------------!
1558	SUBROUTINE set_palm_origin( cfg, namelist_longitude, namelist_latitude, &
1559	origin_lon, origin_lat, z0 )
1560
1561	TYPE(inifor_config), INTENT(IN) :: cfg
1562	REAL(wp), INTENT(IN) :: namelist_longitude, namelist_latitude
1563	REAL(wp), INTENT(OUT) :: origin_lon, origin_lat, z0
1564
1565	message = 'Reading PALM-4U origin from'
1566	IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
1567
1568	origin_lon = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LON ) )
1569	origin_lat = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LAT ) )
1570	z0 = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_Z ) )
1571
1572	message = TRIM(message) // " static driver file '" &
1573	// TRIM( cfg%static_driver_file ) // "'"
1574
1575
1576	ELSE
1577
1578	origin_lon = namelist_longitude
1579	origin_lat = namelist_latitude
1580
1581	message = TRIM( message ) // " namlist file '" &
1582	// TRIM( cfg%namelist_file ) // "'"
1583
1584	ENDIF
1585	origin_lon = origin_lon * TO_RADIANS
1586	origin_lat = origin_lat * TO_RADIANS
1587
1588	CALL report('set_palm_origin', message)
1589
1590	IF ( cfg%z0_is_set ) THEN
1591	z0 = cfg%z0
1592	IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
1593	message = 'You specified both --static-driver/-t and --elevation/-z0. ' // &
1594	'Using the command line value (' // TRIM( real_to_str_f( cfg%z0 ) ) // &
1595	') instead of static driver value (' // TRIM( real_to_str_f( z0 ) ) // ').'
1596	CALL warn( 'set_palm_origin', message )
1597	ENDIF
1598	ENDIF
1599
1600	END SUBROUTINE set_palm_origin
1601
1602
1603	!------------------------------------------------------------------------------!
1604	! Description:
1605	! ------------
1606	! This function is meant to check weather a COSMO soil variable has an
1607	! additional and redunant surface value at depth = 0.0. For instance operational
1608	! DWD COSMO output contains the surface temperature in T_SO as a copy of the
1609	! values in the first soil layer.
1610	!------------------------------------------------------------------------------!
1611	LOGICAL FUNCTION has_surface_value( soil_var, filename )
1612
1613	TYPE(nc_var), INTENT(IN) :: soil_var
1614	CHARACTER(LEN=*), INTENT(IN) :: filename
1615
1616	REAL(wp), ALLOCATABLE :: depths(:)
1617
1618	CALL get_dimension_vector_of_variable( &
1619	soil_var%name, &
1620	dim_idx = NC_DEPTH_DIM_IDX, &
1621	filename = filename, &
1622	dim_vector = depths &
1623	)
1624
1625	has_surface_value = nearly_equal( depths(1), 0.0_wp, 10 * EPSILON(1.0_wp) )
1626
1627	END FUNCTION has_surface_value
1628
1629
1630	!------------------------------------------------------------------------------!
1631	! Description:
1632	! ------------
1633	! This routine reads the dim_idx-th dimension vector of the variable varname
1634	! from netCDF file filename. It is used for finding the depth coordinate vector
1635	! of COSMO soil variables without knowing its name.
1636	!------------------------------------------------------------------------------!
1637	SUBROUTINE get_dimension_vector_of_variable( varname, dim_idx, filename, dim_vector )
1638	CHARACTER(LEN=*), INTENT(IN) :: varname, filename
1639	INTEGER, INTENT(IN) :: dim_idx
1640
1641	REAL(wp), INTENT(OUT), ALLOCATABLE :: dim_vector(:)
1642
1643	INTEGER :: dimids(NF90_MAX_VAR_DIMS)
1644	INTEGER :: varid
1645	CHARACTER(LEN=NF90_MAX_NAME) :: dimname
1646
1647	INTEGER :: ncid
1648
1649	IF ( nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
1650
1651	!
1652	!-- get id of variable varname
1653	CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1654
1655	!
1656	!-- get dimension ids of variable with varid
1657	CALL check( nf90_inquire_variable( ncid, varid, dimids = dimids ) )
1658
1659	!
1660	!-- get name of dim_idx-th dimension variable
1661	CALL check( nf90_inquire_dimension( ncid, dimids(dim_idx), name = dimname ) )
1662	CALL check( nf90_close( ncid ) )
1663
1664	ELSE
1665
1666	message = "Failed to open file '" // TRIM(filename) // "'."
1667	CALL inifor_abort('get_netcdf_variable', message)
1668
1669	ENDIF
1670
1671	! get dimension vector with dimname
1672	CALL get_netcdf_dim_vector( filename, dimname, dim_vector )
1673
1674	END SUBROUTINE get_dimension_vector_of_variable
1675
1676
1677	LOGICAL FUNCTION netcdf_variable_present_in_file( varname, filename )
1678	CHARACTER(LEN=SNAME), INTENT(IN) :: varname
1679	CHARACTER(LEN=PATH), INTENT(IN) :: filename
1680
1681	INTEGER :: ncid, varid
1682
1683	netcdf_variable_present_in_file = ( &
1684	nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
1685	nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR .AND. &
1686	nf90_close( ncid ) .EQ. NF90_NOERR &
1687	)
1688
1689	END FUNCTION netcdf_variable_present_in_file
1690
1691
1692	END MODULE inifor_io
1693	#endif

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |