Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |