Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |