Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |