Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |