Home

Context Navigation

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

Last change on this file since 4619 was 4569, checked in by eckhard, 4 years ago
Removed unused variable
Property svn:keywords set to `Id`
File size: 62.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |