Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |