Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |