Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |