Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |