libsim Versione 7.2.6
gridinfo_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18#include "config.h"
19!> Class for managing information about a single gridded georeferenced
20!! field, typically imported from an on-disk dataset like a grib file
21!! (grib_api driver) or a file in a gdal-supported format (gdal
22!! driver). This module defines a gridinfo (\a gridinfo_def TYPE)
23!! class which can contain information about a single field on a
24!! rectangular georeferenced grid, including:
25!!
26!! - geographical area and projection associated to the gridded data
27!! - time dimension descriptor
28!! - timerange (forecast, analysis, statistically processed) dimension
29!! descriptor
30!! - vertical level dimension descriptor
31!! - physical variable dimension descriptor
32!!
33!! every object contains also an identificator of the grid (\a grid_id
34!! object), carrying information about the driver used or which has to
35!! be used for import/export from/to file. The identificator should be
36!! associated to the gridinfo object at initialization time.
37!!
38!! The main methods of this class allow to:
39!!
40!! - extract (import methods) the information of the five dimension
41!! descriptors from the associated grid_id coming from a file
42!! (grib_api or gdal drivers)
43!! - insert (export method) the information of the five dimension
44!! descriptors in the associated grid_id object, for successive
45!! write to a file (grib_api driver only)
46!! - retrieve the two-dimensional field encoded into the grid_id
47!! associated to the gridinfo object
48!! - encode a desired two-dimensional field into the grid_id
49!! associated to the gridinfo object
50!!
51!! Simple example of use: \include example_vg6d_2.f90
52!! More complex example: \include example_vg6d_4.f90
53!! \ingroup volgrid6d
54MODULE gridinfo_class
55
56USE grid_class
61#ifdef HAVE_LIBGRIBAPI
62USE grib_api
63#endif
64#ifdef HAVE_LIBGDAL
65USE gdal
66#endif
70
71
72IMPLICIT NONE
73
74character (len=255),parameter:: subcategory="gridinfo_class"
75
76
77!> Object describing a single gridded message/band
78TYPE gridinfo_def
79 TYPE(griddim_def) :: griddim !< grid descriptor
80 TYPE(datetime) :: time !< time dimension descriptor
81 TYPE(vol7d_timerange) :: timerange !< timerange (forecast, analysis, statistically processed) dimension descriptor
82 TYPE(vol7d_level) :: level !< vertical level dimension descriptor
83 TYPE(volgrid6d_var) :: var !< physical variable dimension descriptor
84 TYPE(grid_id) :: gaid !< grid identificator, carrying information about the driver for importation/exportation from/to file
85 INTEGER :: category = 0 !< log4fortran category
86END TYPE gridinfo_def
87
88INTEGER, PARAMETER :: &
89 cosmo_centre(3) = (/78,80,200/), & ! emission centres using COSMO coding
90 ecmwf_centre(1) = (/98/) ! emission centres using ECMWF coding
91
92!> Constructor, it creates a new instance of the object.
93INTERFACE init
94 MODULE PROCEDURE gridinfo_init
95END INTERFACE
96
97!> Destructor, it releases every information associated with the object.
98INTERFACE delete
99 MODULE PROCEDURE gridinfo_delete
100END INTERFACE
101
102!> Clone the object, creating a new independent instance of the object
103!! exactly equal to the starting one.
104INTERFACE clone
105 MODULE PROCEDURE gridinfo_clone
106END INTERFACE
107
108!> Import information from a file or grid_id object into the gridinfo
109!! descriptors.
110INTERFACE import
111 MODULE PROCEDURE gridinfo_import, gridinfo_import_from_file
112! MODULE PROCEDURE import_time,import_timerange,import_level,import_gridinfo, &
113! import_volgrid6d_var,gridinfo_import_from_file
114END INTERFACE
115
116!> Export gridinfo descriptors information into a grid_id object.
117INTERFACE export
118 MODULE PROCEDURE gridinfo_export, gridinfo_export_to_file
119! MODULE PROCEDURE export_time,export_timerange,export_level,export_gridinfo, &
120! export_volgrid6d_var
121END INTERFACE
122
123!> Display on standard output a description of the \a gridinfo object
124!! provided.
125INTERFACE display
126 MODULE PROCEDURE gridinfo_display, gridinfov_display
127END INTERFACE
128
129!> Decode and return the data array from a grid_id object associated
130!! to a gridinfo object.
131INTERFACE decode_gridinfo
132 MODULE PROCEDURE gridinfo_decode_data
133END INTERFACE
134
135!> Encode a data array into a grid_id object associated to a gridinfo object.
136INTERFACE encode_gridinfo
137 MODULE PROCEDURE gridinfo_encode_data
138END INTERFACE
139
140#define ARRAYOF_ORIGTYPE TYPE(gridinfo_def)
141#define ARRAYOF_TYPE arrayof_gridinfo
142#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
143#include "arrayof_pre.F90"
144! from arrayof
146
147PRIVATE
150
151CONTAINS
152
153#include "arrayof_post.F90"
154
155!> Constructor, it creates a new instance of the object.
156!! All the additional parameters are optional and they will be
157!! initialised to the corresponding missing value if not provided.
158SUBROUTINE gridinfo_init(this, gaid, griddim, time, timerange, level, var, &
159 clone, categoryappend)
160TYPE(gridinfo_def),intent(out) :: this !< object to be initialized
161TYPE(grid_id),intent(in),optional :: gaid !< identificator of the grid to be described
162type(griddim_def),intent(in),optional :: griddim !< grid descriptor
163TYPE(datetime),intent(in),optional :: time !< time dimension descriptor
164TYPE(vol7d_timerange),intent(in),optional :: timerange !< timerange (forecast, analysis, statistically processed) dimension descriptor
165TYPE(vol7d_level),intent(in),optional :: level !< vertical level dimension descriptor
166TYPE(volgrid6d_var),intent(in),optional :: var !< physical variable dimension descriptor
167logical , intent(in),optional :: clone !< if provided and \c .TRUE., the \a gaid will be cloned and not simply copied into the gridinfo object
168character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
169
170character(len=512) :: a_name
171
172if (present(categoryappend))then
173 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
174else
175 call l4f_launcher(a_name,a_name_append=trim(subcategory))
176end if
177this%category=l4f_category_get(a_name)
178
179#ifdef DEBUG
180call l4f_category_log(this%category,l4f_debug,"start init gridinfo")
181#endif
182
183if (present(gaid))then
184 if (optio_log(clone))then
185 CALL copy(gaid,this%gaid)
186 else
187 this%gaid = gaid
188 endif
189else
190 this%gaid = grid_id_new()
191end if
192
193!#ifdef DEBUG
194!call l4f_category_log(this%category,L4F_DEBUG,"gaid present: "&
195! //to_char(c_e(this%gaid))//" value: "//to_char(this%gaid))
196!#endif
197
198if (present(griddim))then
199 call copy(griddim,this%griddim)
200else
201 call init(this%griddim,categoryappend=categoryappend)
202end if
203
204if (present(time))then
205 this%time=time
206else
207 call init(this%time)
208end if
209
210if (present(timerange))then
211 this%timerange=timerange
212else
213 call init(this%timerange)
214end if
215
216if (present(level))then
217 this%level=level
218else
219 call init(this%level)
220end if
221
222if (present(var))then
223 this%var=var
224else
225 call init(this%var)
226end if
227
228END SUBROUTINE gridinfo_init
229
230
231!> Destructor, it releases every information associated with the object.
232!! It releases memory and deletes the category for logging.
233SUBROUTINE gridinfo_delete(this)
234TYPE(gridinfo_def),intent(inout) :: this !< object to be deleted
235
236#ifdef DEBUG
237call l4f_category_log(this%category,l4f_debug,"start delete_gridinfo" )
238#endif
239
240call delete(this%griddim)
241call delete(this%time)
242call delete(this%timerange)
243call delete(this%level)
244call delete(this%var)
245
246#ifdef DEBUG
247call l4f_category_log(this%category,l4f_debug,"delete gaid" )
248#endif
249CALL delete(this%gaid)
250
251!chiudo il logger
252call l4f_category_delete(this%category)
253
254END SUBROUTINE gridinfo_delete
255
256
257!> Display on standard output a description of the \a gridinfo object
258!! provided. For objects imported from grib_api also the grib key
259!! names and values are printed; the set of keys returned can be
260!! controlled with the input variable namespace. Available namespaces
261!! are "ls", to get the same default keys as the grib_ls command, and
262!! "mars" to get the keys used by mars, default "ls".
263SUBROUTINE gridinfo_display(this, namespace)
264TYPE(gridinfo_def),intent(in) :: this !< object to display
265CHARACTER (len=*),OPTIONAL :: namespace !< grib_api namespace of the keys to search for, all the keys if empty, default "ls"
267#ifdef DEBUG
268call l4f_category_log(this%category,l4f_debug,"displaying gridinfo " )
269#endif
271print*,"----------------------- gridinfo display ---------------------"
272call display(this%griddim)
273call display(this%time)
274call display(this%timerange)
275call display(this%level)
276call display(this%var)
277call display(this%gaid, namespace=namespace)
278print*,"--------------------------------------------------------------"
279
280END SUBROUTINE gridinfo_display
282!> The same as gridinfo_display(), but it receives an array of
283!! gridinfo objects.
284SUBROUTINE gridinfov_display(this, namespace)
285TYPE(arrayof_gridinfo),INTENT(in) :: this !< object to display
286CHARACTER (len=*),OPTIONAL :: namespace !< grib_api namespace of the keys to search for, all the keys if empty, default "ls"
287
288INTEGER :: i
289
290print*,"----------------------- gridinfo array -----------------------"
291
292DO i = 1, this%arraysize
293
294#ifdef DEBUG
295 CALL l4f_category_log(this%array(i)%category,l4f_debug, &
296 "displaying gridinfo array, element "//t2c(i))
297#endif
298 CALL display(this%array(i), namespace)
299
300ENDDO
301print*,"--------------------------------------------------------------"
302
303END SUBROUTINE gridinfov_display
304
306!> Clone the object, creating a new independent instance of the object
307!! exactly equal to the starting one.
308SUBROUTINE gridinfo_clone(this, that, categoryappend)
309TYPE(gridinfo_def),INTENT(in) :: this !< source object
310TYPE(gridinfo_def),INTENT(out) :: that !< cloned object, it does not need to be initalized
311CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category for the cloned object
312
313CALL init(that, gaid=this%gaid, griddim=this%griddim, time=this%time, &
314 timerange=this%timerange, level=this%level, var=this%var, clone=.true., &
315 categoryappend=categoryappend)
316
317END SUBROUTINE gridinfo_clone
318
320!> Import grid_id information into a gridinfo object.
321!! This method imports into the descriptors of the gridinfo object \a
322!! this the information carried on by the grid_id object \a this%gaid,
323!! previously set, typically by reading from a file with a supported
324!! driver (e.g. grib_api or gdal). An amount of information is deduced
325!! from this%gaid and stored in the descriptors of gridinfo object \a
326!! this.
327SUBROUTINE gridinfo_import(this)
328TYPE(gridinfo_def),INTENT(inout) :: this !< gridinfo object
329
330#ifdef HAVE_LIBGRIBAPI
331INTEGER :: gaid
332#endif
333#ifdef HAVE_LIBGDAL
334TYPE(gdalrasterbandh) :: gdalid
335#endif
336
337#ifdef DEBUG
338call l4f_category_log(this%category,l4f_debug,"import from gaid")
339#endif
340
341! griddim is imported separately in grid_class
342CALL import(this%griddim, this%gaid)
343
344#ifdef HAVE_LIBGRIBAPI
345gaid = grid_id_get_gaid(this%gaid)
346IF (c_e(gaid)) CALL gridinfo_import_gribapi(this, gaid)
347#endif
348#ifdef HAVE_LIBGDAL
349gdalid = grid_id_get_gdalid(this%gaid)
350IF (gdalassociated(gdalid)) CALL gridinfo_import_gdal(this, gdalid)
351#endif
352
353END SUBROUTINE gridinfo_import
354
355
356!> Import an array of gridinfo from a file. It receives a (possibly unallocated)
357!! array of gridinfo objects which will be extended by a number of
358!! elements equal to the number of gridded messages/bands found in the
359!! file provided and it will be filled with all the data found. In
360!! case of error, the gridinfo object will not be allocated, so the
361!! success can be tested by checking \a this%arraysize.
362SUBROUTINE gridinfo_import_from_file(this, filename, categoryappend)
363TYPE(arrayof_gridinfo) :: this !< array of gridinfo objects which will be allocated/extended and into which data will be imported
364CHARACTER(len=*),INTENT(in) :: filename !< name of file to open and import, in the form [driver:]pathname
365CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
367type(gridinfo_def) :: gridinfol
368INTEGER :: ngrid, category
369CHARACTER(len=512) :: a_name
370TYPE(grid_file_id) :: input_file
371TYPE(grid_id) :: input_grid
372
373IF (PRESENT(categoryappend)) THEN
374 CALL l4f_launcher(a_name,a_name_append= &
375 trim(subcategory)//"."//trim(categoryappend))
376ELSE
377 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
378ENDIF
379category=l4f_category_get(a_name)
381#ifdef DEBUG
382CALL l4f_category_log(category,l4f_debug,"import from file")
383#endif
384
385input_file = grid_file_id_new(filename, 'r')
387ngrid = 0
388!$OMP PARALLEL DEFAULT(SHARED)
389!$OMP MASTER
390DO WHILE(.true.)
391 input_grid = grid_id_new(input_file)
392 IF (.NOT. c_e(input_grid)) EXIT
393
394 CALL l4f_category_log(category,l4f_info,"import gridinfo")
395 ngrid = ngrid + 1
396!$OMP TASK FIRSTPRIVATE(input_grid, ngrid), PRIVATE(gridinfol)
397 IF (PRESENT(categoryappend)) THEN
398 CALL init(gridinfol, gaid=input_grid, &
399 categoryappend=trim(categoryappend)//"-msg"//trim(to_char(ngrid)))
400 ELSE
401 CALL init(gridinfol, gaid=input_grid, &
402 categoryappend="msg"//trim(to_char(ngrid)))
403 ENDIF
404 CALL import(gridinfol)
405!$OMP CRITICAL
406 CALL insert(this, gridinfol)
407! gridinfol is intentionally not destroyed, since now it lives into this
408!$OMP END CRITICAL
409!$OMP END TASK
410ENDDO
411!$OMP END MASTER
412!$OMP END PARALLEL
413
414CALL packarray(this)
415
416CALL l4f_category_log(category,l4f_info, &
417 "gridinfo_import, "//t2c(ngrid)//" messages/bands imported from file "// &
418 trim(filename))
419
420! close file
421CALL delete(input_file)
422! close logger
423CALL l4f_category_delete(category)
424
425END SUBROUTINE gridinfo_import_from_file
426
427
428!> Export gridinfo descriptors information into a message/band on file.
429!! This method exports the contents of the descriptors of the gridinfo
430!! object \a this in the grid_id object \a this%gaid, previously set,
431!! for the successive write to a file. The information stored in the
432!! descriptors of gridinfo object \a this is inserted, when possible,
433!! in the grid_id object.
434SUBROUTINE gridinfo_export(this)
435TYPE(gridinfo_def),INTENT(inout) :: this !< gridinfo object
436
437#ifdef HAVE_LIBGRIBAPI
438INTEGER :: gaid
439#endif
440#ifdef HAVE_LIBGDAL
441!TYPE(gdalrasterbandh) :: gdalid
442#endif
443
444#ifdef DEBUG
445call l4f_category_log(this%category,l4f_debug,"export to gaid" )
446#endif
447
448! griddim is exported separately in grid_class
449CALL export(this%griddim, this%gaid)
450
451#ifdef HAVE_LIBGRIBAPI
452IF (grid_id_get_driver(this%gaid) == 'grib_api') THEN
453 gaid = grid_id_get_gaid(this%gaid)
454 IF (c_e(gaid)) CALL gridinfo_export_gribapi(this, gaid)
455ENDIF
456#endif
457#ifdef HAVE_LIBGDAL
458IF (grid_id_get_driver(this%gaid) == 'gdal') THEN
459!gdalid = grid_id_get_gdalid(this%gaid)
460 CALL l4f_category_log(this%category,l4f_warn,"export to gdal not implemented" )
461ENDIF
462#endif
463
464END SUBROUTINE gridinfo_export
465
466
467!> Export an arrayof_gridinfo object to a file.
468!! It receives an \a arrayof_gridinfo object which will be exported to
469!! the given file. The driver for writing to file is chosen according
470!! to the gaid associated to the first gridinfo element, and it must
471!! be the same for all the elements.
472SUBROUTINE gridinfo_export_to_file(this, filename, categoryappend)
473TYPE(arrayof_gridinfo) :: this !< array of gridinfo objects which will be written to file
474CHARACTER(len=*),INTENT(in) :: filename !< name of file to open and import, in the form [driver:]pathname
475CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
476
477INTEGER :: i, category
478CHARACTER(len=512) :: a_name
479TYPE(grid_file_id) :: output_file
480TYPE(grid_id) :: valid_grid_id
481
482IF (PRESENT(categoryappend)) THEN
483 CALL l4f_launcher(a_name,a_name_append= &
484 trim(subcategory)//"."//trim(categoryappend))
485ELSE
486 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
487ENDIF
488category=l4f_category_get(a_name)
489
490#ifdef DEBUG
491CALL l4f_category_log(category,l4f_debug, &
492 "exporting to file "//trim(filename)//" "//t2c(this%arraysize)//" fields")
493#endif
494
495valid_grid_id = grid_id_new()
496DO i = 1, this%arraysize ! find a valid grid_id in this
497 IF (c_e(this%array(i)%gaid)) THEN
498 valid_grid_id = this%array(i)%gaid
499 EXIT
500 ENDIF
501ENDDO
502
503IF (c_e(valid_grid_id)) THEN ! a valid grid_id has been found
504! open file
505 output_file = grid_file_id_new(filename, 'w', from_grid_id=valid_grid_id)
506 IF (c_e(output_file)) THEN
507!$OMP PARALLEL DEFAULT(SHARED)
508!$OMP MASTER
509 DO i = 1, this%arraysize
510!$OMP TASK FIRSTPRIVATE(i)
511 CALL export(this%array(i)) ! export information to gaid
512!$OMP CRITICAL
513 CALL export(this%array(i)%gaid, output_file) ! export gaid to file
514!$OMP END CRITICAL
515!$OMP END TASK
516 ENDDO
517!$OMP END MASTER
518!$OMP END PARALLEL
519! close file
520 CALL delete(output_file)
521 ELSE
522 CALL l4f_category_log(category,l4f_error,"opening file "//trim(filename))
523 CALL raise_error()
524 ENDIF
525ELSE ! no valid grid_id has been found
526 CALL l4f_category_log(category,l4f_error, &
527 "gridinfo object of size "//t2c(this%arraysize))
528 CALL l4f_category_log(category,l4f_error, &
529 "no valid grid id found when exporting to file "//trim(filename))
530 CALL raise_error()
531ENDIF
532
533!chiudo il logger
534CALL l4f_category_delete(category)
535
536END SUBROUTINE gridinfo_export_to_file
537
538
539!> Decode and return the data array from a grid_id object associated
540!! to a gridinfo object. This method returns a 2-d array of proper
541!! size extracted from the grid_id object associated to a gridinfo
542!! object. This can work if the gridinfo object has been correctly
543!! initialised and associated to a grid from an on-disk dataset
544!! (e.g. grib_api or gdal file). The result is an array of size \a
545!! this%griddim%dim%nx X \a this%griddim%dim%ny so it must have been
546!! properly allocated by the caller.
547FUNCTION gridinfo_decode_data(this) RESULT(field)
548TYPE(gridinfo_def),INTENT(in) :: this !< gridinfo object
549REAL :: field(this%griddim%dim%nx, this%griddim%dim%ny) ! array of decoded values
550
551CALL grid_id_decode_data(this%gaid, field)
552
553END FUNCTION gridinfo_decode_data
554
555
556!> Encode a data array into a grid_id object associated to a gridinfo object.
557!! This method encodes a 2-d array of proper size into the grid_id
558!! object associated to a gridinfo object. This can work if the
559!! gridinfo object has been correctly initialised and associated to a
560!! grid_id from an on-disk (template) dataset (grib_api or gdal
561!! file). The shape of the array must be conformal to the size of the
562!! grid previously set in the gridinfo object descriptors.
563SUBROUTINE gridinfo_encode_data(this, field)
564TYPE(gridinfo_def),INTENT(inout) :: this !< gridinfo object
565REAL,intent(in) :: field(:,:) !< data array to be encoded
566
567IF (SIZE(field,1) /= this%griddim%dim%nx &
568 .OR. SIZE(field,2) /= this%griddim%dim%ny) THEN
569 CALL l4f_category_log(this%category,l4f_error, &
570 'gridinfo_encode: field and gridinfo object non conformal, field: ' &
571 //trim(to_char(SIZE(field,1)))//'X'//trim(to_char(SIZE(field,2)))//', nx,ny:' &
572 //trim(to_char(this%griddim%dim%nx))//'X'//trim(to_char(this%griddim%dim%ny)))
573 CALL raise_error()
574 RETURN
575ENDIF
576
577CALL grid_id_encode_data(this%gaid, field)
578
579END SUBROUTINE gridinfo_encode_data
580
581
582! =========================================
583! grib_api driver specific code
584! could this be moved to a separate module?
585! =========================================
586#ifdef HAVE_LIBGRIBAPI
587SUBROUTINE gridinfo_import_gribapi(this, gaid)
588TYPE(gridinfo_def),INTENT(inout) :: this ! gridinfo object
589INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
590
591call time_import_gribapi(this%time, gaid)
592call timerange_import_gribapi(this%timerange,gaid)
593call level_import_gribapi(this%level, gaid)
594call var_import_gribapi(this%var, gaid)
595
596call normalize_gridinfo(this)
597
598END SUBROUTINE gridinfo_import_gribapi
599
600
601! grib_api driver
602SUBROUTINE gridinfo_export_gribapi(this, gaid)
603TYPE(gridinfo_def),INTENT(inout) :: this ! gridinfo object
604INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to export
605
606TYPE(conv_func) :: c_func
607REAL,ALLOCATABLE :: tmparr(:,:)
608
609! convert variable and values to the correct edition if required
610CALL volgrid6d_var_normalize(this%var, c_func, grid_id_new(grib_api_id=gaid))
611IF (this%var == volgrid6d_var_miss) THEN
612 CALL l4f_log(l4f_error, &
613 'A suitable variable has not been found in table when converting template')
614 CALL raise_error()
615ENDIF
616IF (c_func /= conv_func_miss) THEN ! convert values as well
617 tmparr = decode_gridinfo(this) ! f2003 implicit allocation
618 CALL compute(c_func, tmparr)
619 CALL encode_gridinfo(this, tmparr)
620ENDIF
621
622CALL unnormalize_gridinfo(this)
623
624CALL time_export_gribapi(this%time, gaid, this%timerange)
625CALL timerange_export_gribapi(this%timerange, gaid, this%time)
626CALL level_export_gribapi(this%level, gaid)
627CALL var_export_gribapi(this%var, gaid)
628
629END SUBROUTINE gridinfo_export_gribapi
630
631
632SUBROUTINE time_import_gribapi(this,gaid)
633TYPE(datetime),INTENT(out) :: this ! datetime object
634INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
635
636INTEGER :: EditionNumber, ttimeincr, tprocdata, centre, p2g, p2, unit, status
637CHARACTER(len=9) :: date
638CHARACTER(len=10) :: time
639
640CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
641
642IF (editionnumber == 1 .OR. editionnumber == 2) THEN
643
644 CALL grib_get(gaid,'dataDate',date )
645 CALL grib_get(gaid,'dataTime',time(:5) )
646
647 CALL init(this,simpledate=date(:8)//time(:4))
648
649 IF (editionnumber == 2) THEN
650
651 CALL grib_get(gaid,'typeOfProcessedData',tprocdata,status)
652 CALL grib_get(gaid,'typeOfTimeIncrement',ttimeincr,status)
653 IF (ttimeincr == 255) ttimeincr = 2 ! fix some MeteosWiss data
654! if analysis-like statistically processed data is encountered, the
655! reference time must be shifted to the end of the processing period
656 IF (status == grib_success .AND. ttimeincr == 1) THEN
657! old libsim convention, to be removed sometime in the future
658 CALL grib_get(gaid,'lengthOfTimeRange',p2g)
659 CALL grib_get(gaid,'indicatorOfUnitForTimeRange',unit)
660 CALL g2_interval_to_second(unit, p2g, p2)
661 this = this + timedelta_new(sec=p2)
662 ELSE IF (status == grib_success .AND. ttimeincr == 2 .AND. tprocdata == 0) THEN
663! generally accepted grib2 convention, DWD exception for cosmo
664! "accumulated" analysis is such that reftime points to the end of the
665! interval, so no time shift in that case
666 CALL grib_get(gaid,'lengthOfTimeRange',p2g)
667 CALL grib_get(gaid,'indicatorOfUnitForTimeRange',unit)
668 CALL g2_interval_to_second(unit, p2g, p2)
669 CALL grib_get(gaid,'centre',centre)
670 IF (centre /= 78) THEN
671 this = this + timedelta_new(sec=p2)
672 ENDIF
673 ELSE IF ((status == grib_success .AND. ttimeincr == 2) .OR. &
674 status /= grib_success) THEN ! usual case
675! do nothing
676 ELSE ! valid but unsupported typeOfTimeIncrement
677 CALL l4f_log(l4f_error,'typeOfTimeIncrement '//t2c(ttimeincr)// &
678 ' not supported')
679 CALL raise_error()
680 ENDIF
681 ENDIF
682
683ELSE
684 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
685 CALL raise_error()
686
687ENDIF
688
689END SUBROUTINE time_import_gribapi
691
692SUBROUTINE time_export_gribapi(this, gaid, timerange)
693TYPE(datetime),INTENT(in) :: this ! datetime object
694INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to export
695TYPE(vol7d_timerange) :: timerange ! timerange, used for grib2 coding of statistically processed analysed data
696
697INTEGER :: EditionNumber, centre
698CHARACTER(len=8) :: env_var
699LOGICAL :: g2cosmo_behavior
700
701CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
702
703IF (editionnumber == 1) THEN
704
705 CALL code_referencetime(this)
706
707ELSE IF (editionnumber == 2 )THEN
708
709 IF (timerange%p1 >= timerange%p2) THEN ! forecast-like
710 CALL code_referencetime(this)
711 ELSE IF (timerange%p1 == 0) THEN ! analysis-like
712! ready for coding with general convention
713 CALL getenv('LIBSIM_G2COSMO_BEHAVIOR', env_var)
714 g2cosmo_behavior = len_trim(env_var) > 0
715 CALL grib_get(gaid,'centre',centre)
716 IF (g2cosmo_behavior .AND. centre == 78) THEN ! DWD analysis exception
717 CALL code_referencetime(this)
718 ELSE ! cosmo or old simc convention
719 CALL code_referencetime(this-timedelta_new(sec=timerange%p2))
720 ENDIF
721 ELSE ! bad timerange
722 CALL l4f_log( l4f_error, 'Timerange with 0>p1>p2 cannot be exported in grib2')
723 CALL raise_error()
724 ENDIF
725
726ELSE
727
728 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
729 CALL raise_error()
730
731ENDIF
732
733CONTAINS
734
735SUBROUTINE code_referencetime(reftime)
736TYPE(datetime),INTENT(in) :: reftime
737
738INTEGER :: date,time
739CHARACTER(len=17) :: date_time
740
741! datetime is AAAAMMGGhhmmssmsc
742CALL getval(reftime, simpledate=date_time)
743READ(date_time(:8),'(I8)')date
744READ(date_time(9:12),'(I4)')time
745CALL grib_set(gaid,'dataDate',date)
746CALL grib_set(gaid,'dataTime',time)
747
748END SUBROUTINE code_referencetime
749
750END SUBROUTINE time_export_gribapi
751
752
753SUBROUTINE level_import_gribapi(this, gaid)
754TYPE(vol7d_level),INTENT(out) :: this ! vol7d_level object
755INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
756
757INTEGER :: EditionNumber,level1,l1,level2,l2
758INTEGER :: ltype,ltype1,scalef1,scalev1,ltype2,scalef2,scalev2
759
760call grib_get(gaid,'GRIBEditionNumber',editionnumber)
761
762if (editionnumber == 1)then
763
764 call grib_get(gaid,'indicatorOfTypeOfLevel',ltype)
765 call grib_get(gaid,'topLevel',l1)
766 call grib_get(gaid,'bottomLevel',l2)
767
768 call level_g1_to_g2(ltype,l1,l2,ltype1,scalef1,scalev1,ltype2,scalef2,scalev2)
769
770else if (editionnumber == 2)then
771
772 call grib_get(gaid,'typeOfFirstFixedSurface',ltype1)
773 call grib_get(gaid,'scaleFactorOfFirstFixedSurface',scalef1)
774 call grib_get(gaid,'scaledValueOfFirstFixedSurface',scalev1)
775 IF (scalef1 == -1 .OR. scalev1 == -1) THEN
776 scalef1 = imiss; scalev1 = imiss
777 ENDIF
778
779 call grib_get(gaid,'typeOfSecondFixedSurface',ltype2)
780 call grib_get(gaid,'scaleFactorOfSecondFixedSurface',scalef2)
781 call grib_get(gaid,'scaledValueOfSecondFixedSurface',scalev2)
782 IF (scalef2 == -1 .OR. scalev2 == -1) THEN
783 scalef2 = imiss; scalev2 = imiss
784 ENDIF
785
786else
787
788 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
789 CALL raise_error()
790
791end if
792
793! Convert missing levels and units m -> mm
794call level_g2_to_dballe(ltype1,scalef1,scalev1,ltype2,scalef2,scalev2, &
795 level1,l1,level2,l2)
796
797call init (this,level1,l1,level2,l2)
798
799END SUBROUTINE level_import_gribapi
800
801
802SUBROUTINE level_export_gribapi(this, gaid)
803TYPE(vol7d_level),INTENT(in) :: this ! vol7d_level object
804INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to export
805
806INTEGER :: EditionNumber, ltype1, scalef1, scalev1, ltype2, scalef2, scalev2, &
807 ltype, l1, l2
808
809CALL level_dballe_to_g2(this%level1, this%l1, this%level2, this%l2, &
810 ltype1, scalef1, scalev1, ltype2, scalef2, scalev2)
811
812call grib_get(gaid,'GRIBEditionNumber',editionnumber)
813
814if (editionnumber == 1)then
815
816 CALL level_g2_to_g1(ltype1,scalef1,scalev1,ltype2,scalef2,scalev2,ltype,l1,l2)
817
818 call grib_set(gaid,'indicatorOfTypeOfLevel',ltype)
819! it is important to set topLevel after, otherwise, in case of single levels
820! bottomLevel=0 overwrites topLevel (aliases in grib_api)
821 call grib_set(gaid,'bottomLevel',l2)
822 call grib_set(gaid,'topLevel',l1)
823
824else if (editionnumber == 2)then
825
826 CALL grib_set(gaid,'typeOfFirstFixedSurface',ltype1)
827 IF (.NOT.c_e(scalef1) .OR. .NOT.c_e(scalev1)) THEN ! code missing values correctly
828 CALL grib_set_missing(gaid,'scaleFactorOfFirstFixedSurface')
829 CALL grib_set_missing(gaid,'scaledValueOfFirstFixedSurface')
830 ELSE
831 CALL grib_set(gaid,'scaleFactorOfFirstFixedSurface',scalef1)
832 CALL grib_set(gaid,'scaledValueOfFirstFixedSurface',scalev1)
833 ENDIF
834
835 CALL grib_set(gaid,'typeOfSecondFixedSurface',ltype2)
836 IF (ltype2 == 255 .OR. .NOT.c_e(scalef2) .OR. .NOT.c_e(scalev2)) THEN ! code missing values correctly
837 CALL grib_set_missing(gaid,'scaleFactorOfSecondFixedSurface')
838 CALL grib_set_missing(gaid,'scaledValueOfSecondFixedSurface')
839 ELSE
840 CALL grib_set(gaid,'scaleFactorOfSecondFixedSurface',scalef2)
841 CALL grib_set(gaid,'scaledValueOfSecondFixedSurface',scalev2)
842 ENDIF
843
844else
845
846 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
847 CALL raise_error()
848
849end if
850
851END SUBROUTINE level_export_gribapi
852
853
854SUBROUTINE timerange_import_gribapi(this, gaid)
855TYPE(vol7d_timerange),INTENT(out) :: this ! vol7d_timerange object
856INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
857
858INTEGER :: EditionNumber, tri, unit, p1g, p2g, p1, p2, statproc, &
859 ttimeincr, tprocdata, status
860
861call grib_get(gaid,'GRIBEditionNumber',editionnumber)
862
863IF (editionnumber == 1) THEN
864
865 CALL grib_get(gaid,'timeRangeIndicator',tri)
866 CALL grib_get(gaid,'P1',p1g)
867 CALL grib_get(gaid,'P2',p2g)
868 CALL grib_get(gaid,'indicatorOfUnitOfTimeRange',unit)
869 CALL timerange_g1_to_v7d(tri, p1g, p2g, unit, statproc, p1, p2)
870
871ELSE IF (editionnumber == 2) THEN
872
873 CALL grib_get(gaid,'forecastTime',p1g)
874 CALL grib_get(gaid,'indicatorOfUnitOfTimeRange',unit)
875 CALL g2_interval_to_second(unit, p1g, p1)
876 call grib_get(gaid,'typeOfStatisticalProcessing',statproc,status)
877
878 IF (status == grib_success .AND. statproc >= 0 .AND. statproc < 254) THEN ! statistically processed
879 CALL grib_get(gaid,'lengthOfTimeRange',p2g)
880 CALL grib_get(gaid,'indicatorOfUnitForTimeRange',unit)
881 CALL g2_interval_to_second(unit, p2g, p2)
882
883! for forecast-like timeranges p1 has to be shifted to the end of interval
884 CALL grib_get(gaid,'typeOfProcessedData',tprocdata,status)
885 CALL grib_get(gaid,'typeOfTimeIncrement',ttimeincr)
886 IF (ttimeincr == 2 .AND. tprocdata /= 0) THEN
887 p1 = p1 + p2
888 ELSE
889 IF (p1 > 0) THEN
890 CALL l4f_log(l4f_warn,'Found p1>0 in grib2 analysis data, strange things may happen')
891 ENDIF
892 ENDIF
893 ELSE ! point in time
894 statproc = 254
895 p2 = 0
896
897 ENDIF
898
899ELSE
900
901 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
902 CALL raise_error()
903
904ENDIF
905
906CALL init(this, statproc, p1, p2)
907
908END SUBROUTINE timerange_import_gribapi
909
910
911SUBROUTINE timerange_export_gribapi(this, gaid, reftime)
912TYPE(vol7d_timerange),INTENT(in) :: this ! vol7d_timerange object
913INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to export
914TYPE(datetime) :: reftime ! reference time of data, used for coding correct end of statistical processing period in grib2
915
916INTEGER :: EditionNumber, centre, tri, currentunit, unit, p1_g1, p2_g1, p1, p2, pdtn
917CHARACTER(len=8) :: env_var
918LOGICAL :: g2cosmo_behavior
919
920CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
921
922IF (editionnumber == 1 ) THEN
923! Convert vol7d_timerange members to grib1 with reasonable time unit
924 CALL grib_get(gaid,'indicatorOfUnitOfTimeRange',currentunit)
925 CALL timerange_v7d_to_g1(this%timerange, this%p1, this%p2, &
926 tri, p1_g1, p2_g1, unit)
927! Set the native keys
928 CALL grib_set(gaid,'timeRangeIndicator',tri)
929 CALL grib_set(gaid,'P1',p1_g1)
930 CALL grib_set(gaid,'P2',p2_g1)
931 CALL grib_set(gaid,'indicatorOfUnitOfTimeRange',unit)
932
933ELSE IF (editionnumber == 2) THEN
934 CALL grib_get(gaid,'productDefinitionTemplateNumber', pdtn)
935
936 IF (this%timerange == 254) THEN ! point in time -> template 4.0
937 IF (pdtn < 0 .OR. pdtn > 7) &
938 CALL grib_set(gaid,'productDefinitionTemplateNumber', 0)
939! Set reasonable time unit
940 CALL timerange_v7d_to_g2(this%p1,p1,unit)
941! Set the native keys
942 CALL grib_set(gaid,'indicatorOfUnitOfTimeRange',unit)
943 CALL grib_set(gaid,'forecastTime',p1)
944
945 ELSE IF (this%timerange >= 0 .AND. this%timerange < 254) THEN
946! statistically processed -> template 4.8
947 IF (pdtn < 8 .OR. pdtn > 14) &
948 CALL grib_set(gaid,'productDefinitionTemplateNumber', 8)
949
950 IF (this%p1 >= this%p2) THEN ! forecast-like
951! Set reasonable time unit
952 CALL timerange_v7d_to_g2(this%p1-this%p2,p1,unit)
953 CALL grib_set(gaid,'indicatorOfUnitOfTimeRange',unit)
954 CALL grib_set(gaid,'forecastTime',p1)
955 CALL code_endoftimeinterval(reftime+timedelta_new(sec=this%p1))
956! Successive times processed have same start time of forecast,
957! forecast time is incremented
958 CALL grib_set(gaid,'typeOfStatisticalProcessing',this%timerange)
959! typeOfTimeIncrement to be replaced with a check that typeOfProcessedData /= 0
960 CALL grib_set(gaid,'typeOfTimeIncrement',2)
961 CALL timerange_v7d_to_g2(this%p2,p2,unit)
962 CALL grib_set(gaid,'indicatorOfUnitForTimeRange',unit)
963 CALL grib_set(gaid,'lengthOfTimeRange',p2)
964
965 ELSE IF (this%p1 == 0) THEN ! analysis-like
966! Set reasonable time unit
967 CALL timerange_v7d_to_g2(this%p2,p2,unit)
968 CALL grib_set(gaid,'indicatorOfUnitOfTimeRange',unit)
969 CALL grib_set(gaid,'forecastTime',0)
970 CALL code_endoftimeinterval(reftime)
971! Successive times processed have same forecast time, start time of
972! forecast is incremented
973 CALL grib_set(gaid,'typeOfStatisticalProcessing',this%timerange)
974! typeOfTimeIncrement to be replaced with typeOfProcessedData
975 CALL getenv('LIBSIM_G2COSMO_BEHAVIOR', env_var)
976 g2cosmo_behavior = len_trim(env_var) > 0
977 IF (g2cosmo_behavior) THEN
978 CALL grib_set(gaid,'typeOfProcessedData',0)
979 ELSE
980 CALL grib_set(gaid,'typeOfTimeIncrement',1)
981 ENDIF
982 CALL grib_set(gaid,'indicatorOfUnitForTimeRange',unit)
983 CALL grib_set(gaid,'lengthOfTimeRange',p2)
984
985! warn about local use
986 IF (this%timerange >= 192) THEN
987 CALL l4f_log(l4f_warn, &
988 'coding in grib2 a nonstandard typeOfStatisticalProcessing '// &
989 t2c(this%timerange))
990 ENDIF
991 ELSE ! bad timerange
992 CALL l4f_log(l4f_error, &
993 'Timerange with 0>p1>p2 cannot be exported in grib2')
994 CALL raise_fatal_error()
995 ENDIF
996 ELSE
997 CALL l4f_log(l4f_error, &
998 'typeOfStatisticalProcessing not supported: '//trim(to_char(this%timerange)))
999 CALL raise_fatal_error()
1000 ENDIF
1001
1002ELSE
1003 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
1004 CALL raise_fatal_error()
1005ENDIF
1006
1007CONTAINS
1008
1009! Explicitely compute and code in grib2 template 4.8 the end of
1010! overalltimeinterval which is not done automatically by grib_api
1011SUBROUTINE code_endoftimeinterval(endtime)
1012TYPE(datetime),INTENT(in) :: endtime
1013
1014INTEGER :: year, month, day, hour, minute, msec
1015
1016CALL getval(endtime, year=year, month=month, day=day, &
1017 hour=hour, minute=minute, msec=msec)
1018 CALL grib_set(gaid,'yearOfEndOfOverallTimeInterval',year)
1019 CALL grib_set(gaid,'monthOfEndOfOverallTimeInterval',month)
1020 CALL grib_set(gaid,'dayOfEndOfOverallTimeInterval',day)
1021 CALL grib_set(gaid,'hourOfEndOfOverallTimeInterval',hour)
1022 CALL grib_set(gaid,'minuteOfEndOfOverallTimeInterval',minute)
1023 CALL grib_set(gaid,'secondOfEndOfOverallTimeInterval',msec/1000)
1024
1025END SUBROUTINE code_endoftimeinterval
1026
1027END SUBROUTINE timerange_export_gribapi
1028
1029
1030SUBROUTINE var_import_gribapi(this, gaid)
1031TYPE(volgrid6d_var),INTENT(out) :: this ! volgrid6d_var object
1032INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
1033
1034INTEGER :: EditionNumber, centre, discipline, category, number
1035
1036call grib_get(gaid,'GRIBEditionNumber',editionnumber)
1037
1038if (editionnumber == 1) then
1039
1040 call grib_get(gaid,'centre',centre)
1041 call grib_get(gaid,'gribTablesVersionNo',category)
1042 call grib_get(gaid,'indicatorOfParameter',number)
1043
1044 call init(this, centre, category, number)
1045
1046else if (editionnumber == 2) then
1047
1048 call grib_get(gaid,'centre',centre)
1049 call grib_get(gaid,'discipline',discipline)
1050 call grib_get(gaid,'parameterCategory',category)
1051 call grib_get(gaid,'parameterNumber',number)
1052
1053 call init(this, centre, category, number, discipline)
1054
1055else
1056
1057 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
1058 CALL raise_error()
1059
1060endif
1061
1062END SUBROUTINE var_import_gribapi
1063
1064
1065SUBROUTINE var_export_gribapi(this, gaid)
1066TYPE(volgrid6d_var),INTENT(in) :: this ! volgrid6d_var object
1067INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to export
1068
1069INTEGER ::EditionNumber
1070
1071call grib_get(gaid,'GRIBEditionNumber',editionnumber)
1072
1073if (editionnumber == 1) then
1074
1075 IF (this%centre /= 255) & ! if centre missing (coming from bufr), keep template
1076 CALL grib_set(gaid,'centre',this%centre)
1077 CALL grib_set(gaid,'gribTablesVersionNo',this%category)
1078 CALL grib_set(gaid,'indicatorOfParameter',this%number)
1080else if (editionnumber == 2) then
1081
1082! this must be changed to 65535!!!!
1083 IF (this%centre /= 255) & ! if centre missing (coming from bufr), keep template
1084 CALL grib_set(gaid,'centre',this%centre)
1085 CALL grib_set(gaid,'discipline',this%discipline)
1086 CALL grib_set(gaid,'parameterCategory',this%category)
1087 CALL grib_set(gaid,'parameterNumber',this%number)
1088
1089else
1090
1091 CALL l4f_log(l4f_error,'GribEditionNumber '//t2c(editionnumber)//' not supported')
1092 CALL raise_error()
1093
1094end if
1096END SUBROUTINE var_export_gribapi
1097
1098
1099SUBROUTINE level_g2_to_dballe(ltype1,scalef1,scalev1,ltype2,scalef2,scalev2, lt1,l1,lt2,l2)
1100integer,intent(in) :: ltype1,scalef1,scalev1,ltype2,scalef2,scalev2
1101integer,intent(out) ::lt1,l1,lt2,l2
1102
1103
1104CALL g2_to_dballe(ltype1, scalef1, scalev1, lt1, l1)
1105CALL g2_to_dballe(ltype2, scalef2, scalev2, lt2, l2)
1106
1107CONTAINS
1108
1109SUBROUTINE g2_to_dballe(ltype, scalef, scalev, lt, l)
1110integer,intent(in) :: ltype,scalef,scalev
1111integer,intent(out) :: lt,l
1112
1113doubleprecision :: sl
1114
1115
1116IF (ltype == 255 .OR. ltype == -1) THEN
1117 lt = imiss
1118 l = imiss
1119ELSE IF (ltype <= 10 .OR. ltype == 101 .OR. (ltype >= 162 .AND. ltype <= 184)) THEN
1120 lt = ltype
1121 l = imiss
1122ELSE
1123 lt = ltype
1124 IF (c_e(scalef) .AND. c_e(scalev)) THEN
1125 sl = scalev*(10.d0**(-scalef))
1126! apply unit conversions
1127 IF (any(ltype == height_level)) THEN
1128 l = nint(sl*1000.d0)
1129 ELSE IF (any(ltype == thermo_level)) THEN
1130 l = nint(sl*10.d0)
1131 ELSE IF (any(ltype == sigma_level)) THEN
1132 l = nint(sl*10000.d0)
1133 ELSE
1134 l = nint(sl)
1135 ENDIF
1136 ELSE
1137 l = imiss
1138 ENDIF
1139ENDIF
1140
1141END SUBROUTINE g2_to_dballe
1142
1143END SUBROUTINE level_g2_to_dballe
1144
1145
1146SUBROUTINE level_dballe_to_g2(lt1,l1,lt2,l2, ltype1,scalef1,scalev1,ltype2,scalef2,scalev2)
1147integer,intent(in) :: lt1,l1,lt2,l2
1148integer,intent(out) :: ltype1,scalef1,scalev1,ltype2,scalef2,scalev2
1149
1150CALL dballe_to_g2(lt1, l1, ltype1, scalef1, scalev1)
1151CALL dballe_to_g2(lt2, l2, ltype2, scalef2, scalev2)
1152
1153CONTAINS
1154
1155SUBROUTINE dballe_to_g2(lt, l, ltype, scalef, scalev)
1156INTEGER,INTENT(in) :: lt,l
1157INTEGER,INTENT(out) :: ltype,scalef,scalev
1158
1159
1160IF (lt == imiss) THEN
1161 ltype = 255
1162 scalev = 0
1163 scalef = 0
1164ELSE IF (lt <= 10 .OR. lt == 101 .OR. (lt >= 162 .AND. lt <= 184)) THEN
1165 ltype = lt
1166 scalev = 0
1167 scalef = 0
1168ELSE IF (lt == 256 .AND. l == imiss) THEN ! special case for cloud level -> surface
1169 ltype = 1
1170 scalev = 0
1171 scalef = 0
1172ELSE
1173 ltype = lt
1174 scalev = l
1175 IF (any(ltype == height_level)) THEN
1176 scalef = 3
1177 ELSE IF (any(ltype == thermo_level)) THEN
1178 scalef = 1
1179 ELSE IF (any(ltype == sigma_level)) THEN
1180 scalef = 4
1181 ELSE
1182 scalef = 0
1183 ENDIF
1184ENDIF
1185
1186!Caso generale reale
1187!IF (ANY(ltype == height_level)) THEN
1188! sl=l/1000.D0
1189!ELSE
1190! sl=l
1191!END IF
1192!IF (ABS(sl) < PRECISION) THEN
1193! scalef = 0
1194! scalev = 0
1195!ELSE
1196! magn = LOG10(ABS(sl))
1197! DO scalef = magn, MAX(magn, 20)
1198! IF (ABS((sl*10.D0**(scalef) - ANINT(sl*10.D0**(scalef))) / &
1199! sl*10.D0**(scalef)) < PRECISION) EXIT
1200! ENDDO
1201! sl = scalev*(10.D0**(-scalef))
1202!ENDIF
1203
1204END SUBROUTINE dballe_to_g2
1205
1206END SUBROUTINE level_dballe_to_g2
1207
1208
1209SUBROUTINE level_g1_to_g2(ltype,l1,l2,ltype1,scalef1,scalev1,ltype2,scalef2,scalev2)
1210integer,intent(in) :: ltype,l1,l2
1211integer,intent(out) :: ltype1,scalef1,scalev1,ltype2,scalef2,scalev2
1212
1213ltype1=255
1214scalef1=0
1215scalev1=0
1216ltype2=255
1217scalef2=0
1218scalev2=0
1219
1220if (ltype > 0 .and. ltype <= 9)then
1221 ltype1=ltype
1222else if (ltype == 20) then
1223 ltype1=20
1224 scalev1=l1
1225 scalef1=2
1226else if (ltype == 100) then
1227 ltype1=100
1228 scalev1=l1*100
1229else if (ltype == 101) then
1230 ltype1=100
1231 scalev1=l1*1000
1232 ltype2=100
1233 scalev2=l2*1000
1234else if (ltype == 102) then
1235 ltype1=101
1236else if (ltype == 103) then
1237 ltype1=102
1238 scalev1=l1
1239else if (ltype == 104) then
1240 ltype1=102
1241 scalev1=l1*100
1242 ltype2=102
1243 scalev2=l2*100
1244else if (ltype == 105) then
1245 ltype1=103
1246 scalev1=l1
1247else if (ltype == 106) then
1248 ltype1=103
1249 scalev1=l1*100
1250 ltype2=103
1251 scalev2=l2*100
1252else if (ltype == 107) then
1253 ltype1=104
1254 scalef1=4
1255 scalev1=l1
1256else if (ltype == 108) then
1257 ltype1=104
1258 scalef1=2
1259 scalev1=l1
1260 ltype2=104
1261 scalef2=2
1262 scalev2=l2
1263else if (ltype == 109) then
1264 ltype1=105
1265 scalev1=l1
1266else if (ltype == 110) then
1267 ltype1=105
1268 scalev1=l1
1269 ltype2=105
1270 scalev2=l2
1271else if (ltype == 111) then
1272 ltype1=106
1273 scalef1=2
1274 scalev1=l1
1275else if (ltype == 112) then
1276 ltype1=106
1277 scalef1=2
1278 scalev1=l1
1279 ltype2=106
1280 scalef2=2
1281 scalev2=l2
1282else if (ltype == 113) then
1283 ltype1=107
1284 scalev1=l1
1285else if (ltype == 114) then
1286 ltype1=107
1287 scalev1=475+l1
1288 ltype2=107
1289 scalev2=475+l2
1290else if (ltype == 115) then
1291 ltype1=108
1292 scalev1=l1*100
1293else if (ltype == 116) then
1294 ltype1=108
1295 scalev1=l1*100
1296 ltype2=108
1297 scalev2=l2*100
1298else if (ltype == 117) then
1299 ltype1=109
1300 scalef1=9
1301 scalev1=l1
1302 if ( btest(l1,15) ) then
1303 scalev1=-1*mod(l1,32768)
1304 endif
1305else if (ltype == 119) then
1306 ltype1=111
1307 scalef1=4
1308 scalev1=l1
1309else if (ltype == 120) then
1310 ltype1=111
1311 scalef1=2
1312 scalev1=l1
1313 ltype2=111
1314 scalef2=2
1315 scalev2=l2
1316else if (ltype == 121) then
1317 ltype1=100
1318 scalev1=(1100+l1)*100
1319 ltype2=100
1320 scalev2=(1100+l2)*100
1321else if (ltype == 125) then
1322 ltype1=103
1323 scalef1=2
1324 scalev1=l1
1325else if (ltype == 128) then
1326 ltype1=104
1327 scalef1=3
1328 scalev1=1100+l1
1329 ltype2=104
1330 scalef2=3
1331 scalev2=1100+l2
1332else if (ltype == 141) then
1333 ltype1=100
1334 scalev1=l1*100
1335 ltype2=100
1336 scalev2=(1100+l2)*100
1337else if (ltype == 160) then
1338 ltype1=160
1339 scalev1=l1
1340else if (ltype == 200) then ! entire atmosphere
1341 ltype1=1
1342 ltype2=8
1343else if (ltype == 210) then ! entire ocean
1344 ltype1=1
1345 ltype2=9
1346else
1347
1348 call l4f_log(l4f_error,'level_g1_to_g2: GRIB1 level '//trim(to_char(ltype)) &
1349 //' cannot be converted to GRIB2.')
1350 call raise_error()
1351endif
1352
1353END SUBROUTINE level_g1_to_g2
1354
1355
1356SUBROUTINE level_g2_to_g1(ltype1,scalef1,scalev1,ltype2,scalef2,scalev2,ltype,l1,l2)
1357integer,intent(in) :: ltype1,scalef1,scalev1,ltype2,scalef2,scalev2
1358integer,intent(out) :: ltype,l1,l2
1359
1360if (ltype1 > 0 .and. ltype1 <= 9 .and. ltype2 == 255) then ! simple
1361 ltype = ltype1
1362 l1 = 0
1363 l2 = 0
1364else if (ltype1 == 20 .and. ltype2 == 255) then ! isothermal
1365 ltype = 2
1366 l1 = rescale2(scalef1-2,scalev1)!*100
1367 l2 = 0
1368else if (ltype1 == 100 .and. ltype2 == 255) then ! isobaric
1369 ltype = 100
1370 l1 = rescale2(scalef1+2,scalev1)!/100
1371 l2 = 0
1372else if (ltype1 == 100 .and. ltype2 == 100) then
1373 ltype = 101
1374 l1 = rescale1(scalef1+3,scalev1)!/1000
1375 l2 = rescale1(scalef2+3,scalev2)!/1000
1376else if (ltype1 == 101 .and. ltype2 == 255) then
1377 ltype = 102
1378 l1 = 0
1379 l2 = 0
1380else if (ltype1 == 102 .and. ltype2 == 255) then ! altitude over sea level
1381 ltype = 103
1382 l1 = rescale2(scalef1,scalev1)
1383 l2 = 0
1384else if (ltype1 == 102 .and. ltype2 == 102) then
1385 ltype = 104
1386 l1 = rescale1(scalef1+2,scalev1)!/100
1387 l2 = rescale1(scalef2+2,scalev2)!/100
1388else if (ltype1 == 103 .and. ltype2 == 255) then ! height over ground
1389 ltype = 105
1390 l1 = rescale2(scalef1,scalev1)
1391 l2 = 0
1392else if (ltype1 == 103 .and. ltype2 == 103) then
1393 ltype = 106
1394 l1 = rescale1(scalef1+2,scalev1)!/100
1395 l2 = rescale1(scalef2+2,scalev2)!/100
1396else if (ltype1 == 104 .and. ltype2 == 255) then ! sigma
1397 ltype = 107
1398 l1 = rescale2(scalef1,scalev1-4)!*10000
1399 l2 = 0
1400else if (ltype1 == 104 .and. ltype2 == 104) then
1401 ltype = 108
1402 l1 = rescale1(scalef1-2,scalev1)!*100
1403 l2 = rescale1(scalef2-2,scalev2)!*100
1404else if (ltype1 == 105 .and. ltype2 == 255) then ! hybrid
1405 ltype = 109
1406 l1 = rescale2(scalef1,scalev1)
1407 l2 = 0
1408else if (ltype1 == 105 .and. ltype2 == 105) then
1409 ltype = 110
1410 l1 = rescale1(scalef1,scalev1)
1411 l2 = rescale1(scalef2,scalev2)
1412else if (ltype1 == 106 .and. ltype2 == 255) then ! depth
1413 ltype = 111
1414 l1 = rescale2(scalef1-2,scalev1)!*100
1415 l2 = 0
1416else if (ltype1 == 106 .and. ltype2 == 106) then
1417 ltype = 112
1418 l1 = rescale1(scalef1-2,scalev1)!*100
1419 l2 = rescale1(scalef2-2,scalev2)!*100
1420else if (ltype1 == 107 .and. ltype2 == 255) then ! isentropic
1421 ltype = 113
1422 l1 = rescale2(scalef1,scalev1)
1423 l2 = 0
1424else if (ltype1 == 107 .and. ltype2 == 107) then
1425 ltype = 114
1426 l1 = rescale1(scalef1,scalev1)
1427 l2 = rescale1(scalef2,scalev2)
1428else if (ltype1 == 108 .and. ltype2 == 255) then ! pressure diff to ground
1429 ltype = 115
1430 l1 = rescale2(scalef1+2,scalev1)!/100
1431 l2 = 0
1432else if (ltype1 == 108 .and. ltype2 == 108) then
1433 ltype = 116
1434 l1 = rescale1(scalef1+2,scalev1)!/100
1435 l2 = rescale1(scalef2+2,scalev2)!/100
1436else if (ltype1 == 109 .and. ltype2 == 255) then ! potential vorticity surf
1437 ltype = 117
1438 l1 = rescale2(scalef1+9,scalev1)!/10**9
1439 l2 = 0
1440else if (ltype1 == 111 .and. ltype2 == 255) then ! eta level
1441 ltype = 119
1442 l1 = rescale2(scalef1-2,scalev1)!*100
1443 l2 = 0
1444else if (ltype1 == 111 .and. ltype2 == 111) then
1445 ltype = 120
1446 l1 = rescale1(scalef1-4,scalev1)!*10000
1447 l2 = rescale1(scalef2-4,scalev2)!*10000
1448else if (ltype1 == 160 .and. ltype2 == 255) then ! depth below sea
1449 ltype = 160
1450 l1 = rescale2(scalef1,scalev1)
1451 l2 = 0
1452else if (ltype1 == 1 .and. ltype2 == 8) then ! entire atmosphere
1453 ltype = 200
1454else if (ltype1 == 1 .and. ltype2 == 9) then ! entire ocean
1455 ltype = 201
1456else ! mi sono rotto per ora
1457
1458 ltype = 255
1459 l1 = 0
1460 l2 = 0
1461 call l4f_log(l4f_error,'level_g2_to_g1: GRIB2 levels '//trim(to_char(ltype1))//' ' &
1462 //trim(to_char(ltype2))//' cannot be converted to GRIB1.')
1463 call raise_error()
1464endif
1465
1466CONTAINS
1467
1468FUNCTION rescale1(scalef, scalev) RESULT(rescale)
1469INTEGER,INTENT(in) :: scalef, scalev
1470INTEGER :: rescale
1471
1472rescale = min(255, nint(scalev*10.0d0**(-scalef)))
1473
1474END FUNCTION rescale1
1475
1476FUNCTION rescale2(scalef, scalev) RESULT(rescale)
1477INTEGER,INTENT(in) :: scalef, scalev
1478INTEGER :: rescale
1479
1480rescale = min(65535, nint(scalev*10.0d0**(-scalef)))
1481
1482END FUNCTION rescale2
1483
1484END SUBROUTINE level_g2_to_g1
1485
1486! Convert timerange from grib1 to grib2-like way:
1487!
1488! Tri 2 (point in time) gives (hopefully temporarily) statproc 205.
1489!
1490! Tri 13 (COSMO-nudging) gives p1 (forecast time) 0 and a temporary
1491! 257 statproc.
1492!
1493! Further processing and correction of the values returned is
1494! performed in normalize_gridinfo.
1495SUBROUTINE timerange_g1_to_v7d(tri, p1_g1, p2_g1, unit, statproc, p1, p2)
1496INTEGER,INTENT(in) :: tri, p1_g1, p2_g1, unit
1497INTEGER,INTENT(out) :: statproc, p1, p2
1498
1499IF (tri == 0 .OR. tri == 1) THEN ! point in time
1500 statproc = 254
1501 CALL g1_interval_to_second(unit, p1_g1, p1)
1502 p2 = 0
1503ELSE IF (tri == 10) THEN ! point in time
1504 statproc = 254
1505 CALL g1_interval_to_second(unit, p1_g1*256+p2_g1, p1)
1506 p2 = 0
1507ELSE IF (tri == 2) THEN ! somewhere between p1 and p2, not good for grib2 standard
1508 statproc = 205
1509 CALL g1_interval_to_second(unit, p2_g1, p1)
1510 CALL g1_interval_to_second(unit, p2_g1-p1_g1, p2)
1511ELSE IF (tri == 3) THEN ! average
1512 statproc = 0
1513 CALL g1_interval_to_second(unit, p2_g1, p1)
1514 CALL g1_interval_to_second(unit, p2_g1-p1_g1, p2)
1515ELSE IF (tri == 4) THEN ! accumulation
1516 statproc = 1
1517 CALL g1_interval_to_second(unit, p2_g1, p1)
1518 CALL g1_interval_to_second(unit, p2_g1-p1_g1, p2)
1519ELSE IF (tri == 5) THEN ! difference
1520 statproc = 4
1521 CALL g1_interval_to_second(unit, p2_g1, p1)
1522 CALL g1_interval_to_second(unit, p2_g1-p1_g1, p2)
1523ELSE IF (tri == 13) THEN ! COSMO-nudging, use a temporary value then normalize
1524 statproc = 257 ! check if 257 is legal!
1525 p1 = 0 ! analysis regardless of p2_g1
1526 CALL g1_interval_to_second(unit, p2_g1-p1_g1, p2)
1527ELSE
1528 call l4f_log(l4f_error,'timerange_g1_to_g2: GRIB1 timerange '//trim(to_char(tri)) &
1529 //' cannot be converted to GRIB2.')
1530 CALL raise_error()
1531ENDIF
1532
1533if (statproc == 254 .and. p2 /= 0 ) then
1534 call l4f_log(l4f_warn,"inconsistence in timerange:254,"//trim(to_char(p1))//","//trim(to_char(p2)))
1535end if
1536
1537END SUBROUTINE timerange_g1_to_v7d
1538
1539
1540!0 Minute
1541!1 Hour
1542!2 Day
1543!3 Month
1544!4 Year
1545!5 Decade (10 years)
1546!6 Normal (30 years)
1547!7 Century(100 years)
1548!8-9 Reserved
1549!10 3 hours
1550!11 6 hours
1551!12 12 hours
1552! in COSMO, grib1:
1553!13 = 15 minuti
1554!14 = 30 minuti
1555! in grib2:
1556!13 Second
1557
1558SUBROUTINE g1_interval_to_second(unit, valuein, valueout)
1559INTEGER,INTENT(in) :: unit, valuein
1560INTEGER,INTENT(out) :: valueout
1561
1562INTEGER,PARAMETER :: unitlist(0:14)=(/ 60,3600,86400,2592000, &
1563 31536000,315360000,946080000,imiss,imiss,imiss,10800,21600,43200,900,1800/)
1564
1565valueout = imiss
1566IF (unit >= lbound(unitlist,1) .AND. unit <= ubound(unitlist,1)) THEN
1567 IF (c_e(unitlist(unit))) THEN
1568 valueout = valuein*unitlist(unit)
1569 ENDIF
1570ENDIF
1571
1572END SUBROUTINE g1_interval_to_second
1573
1574
1575SUBROUTINE g2_interval_to_second(unit, valuein, valueout)
1576INTEGER,INTENT(in) :: unit, valuein
1577INTEGER,INTENT(out) :: valueout
1578
1579INTEGER,PARAMETER :: unitlist(0:13)=(/ 60,3600,86400,2592000, &
1580 31536000,315360000,946080000,imiss,imiss,imiss,10800,21600,43200,1/)
1581
1582valueout = imiss
1583IF (unit >= lbound(unitlist,1) .AND. unit <= ubound(unitlist,1)) THEN
1584 IF (c_e(unitlist(unit))) THEN
1585 valueout = valuein*unitlist(unit)
1586 ENDIF
1587ENDIF
1588
1589END SUBROUTINE g2_interval_to_second
1590
1591
1592! Convert timerange from grib2-like to grib1 way:
1593!
1594! Statproc 205 (point in time, non standard, not good in grib2) is
1595! correctly converted to tri 2.
1596!
1597! Statproc 257 (COSMO nudging-like, non standard, not good in grib2)
1598! should not appear here, but is anyway converted to tri 13 (real
1599! COSMO-nudging).
1600!
1601! In case p1_g1<0 (i.e. statistically processed analysis data, with
1602! p1=0 and p2>0), COSMO-nudging tri 13 is (mis-)used.
1603SUBROUTINE timerange_v7d_to_g1(statproc, p1, p2, tri, p1_g1, p2_g1, unit)
1604INTEGER,INTENT(in) :: statproc, p1, p2
1605INTEGER,INTENT(out) :: tri, p1_g1, p2_g1, unit
1606
1607INTEGER :: ptmp, pdl
1608
1609pdl = p1 - p2
1610IF (statproc == 254) pdl = p1 ! avoid unexpected situations (necessary?)
1611
1612CALL timerange_choose_unit_g1(p1, pdl, p2_g1, p1_g1, unit)
1613IF (statproc == 0) THEN ! average
1614 tri = 3
1615ELSE IF (statproc == 1) THEN ! accumulation
1616 tri = 4
1617ELSE IF (statproc == 4) THEN ! difference
1618 tri = 5
1619ELSE IF (statproc == 205) THEN ! point in time interval, not good for grib2 standard
1620 tri = 2
1621ELSE IF (statproc == 257) THEN ! COSMO-nudging (statistical processing in the past)
1622! this should never happen (at least from COSMO grib1), since 257 is
1623! converted to something else in normalize_gridinfo; now a negative
1624! p1_g1 is set, it will be corrected in the next section
1625 tri = 13
1626! CALL second_to_gribtr(p1, p2_g1, unit, 1)
1627! CALL second_to_gribtr(p1-p2, p1_g1, unit, 1)
1628ELSE IF (statproc == 254) THEN ! point in time
1629 tri = 0
1630 p2_g1 = 0
1631ELSE
1632 CALL l4f_log(l4f_error,'timerange_v7d_to_g1: GRIB2 statisticalprocessing ' &
1633 //trim(to_char(statproc))//' cannot be converted to GRIB1.')
1634 CALL raise_error()
1635ENDIF
1636
1637IF (p1_g1 > 255 .OR. p2_g1 > 255) THEN
1638 ptmp = max(p1_g1,p2_g1)
1639 p2_g1 = mod(ptmp,256)
1640 p1_g1 = ptmp/256
1641 IF (tri /= 0) THEN ! if not instantaneous give warning
1642 CALL l4f_log(l4f_warn,'timerange_v7d_to_g1: timerange too long for grib1 ' &
1643 //trim(to_char(ptmp))//', forcing time range indicator to 10.')
1644 ENDIF
1645 tri = 10
1646ENDIF
1647
1648
1649! p1 < 0 is not allowed, use COSMO trick
1650IF (p1_g1 < 0) THEN
1651 ptmp = p1_g1
1652 p1_g1 = 0
1653 p2_g1 = p2_g1 - ptmp
1654 tri = 13
1655ENDIF
1656
1657END SUBROUTINE timerange_v7d_to_g1
1658
1659
1660SUBROUTINE timerange_v7d_to_g2(valuein, valueout, unit)
1661INTEGER,INTENT(in) :: valuein
1662INTEGER,INTENT(out) :: valueout, unit
1663
1664IF (valuein == imiss) THEN
1665 valueout = imiss
1666 unit = 1
1667ELSE IF (mod(valuein,3600) == 0) THEN ! prefer hours
1668 valueout = valuein/3600
1669 unit = 1
1670ELSE IF (mod(valuein,60) == 0) THEN ! then minutes
1671 valueout = valuein/60
1672 unit = 0
1673ELSE ! seconds
1674 valueout = valuein
1675 unit = 13
1676ENDIF
1677
1678END SUBROUTINE timerange_v7d_to_g2
1679
1680
1681! These units are tested for applicability:
1682! 0 Minute
1683! 1 Hour
1684! 2 Day
1685! 10 3 hours
1686! 11 6 hours
1687! 12 12 hours
1688SUBROUTINE timerange_choose_unit_g1(valuein1, valuein2, valueout1, valueout2, unit)
1689INTEGER,INTENT(in) :: valuein1, valuein2
1690INTEGER,INTENT(out) :: valueout1, valueout2, unit
1691
1692INTEGER :: i
1693TYPE unitchecker
1694 INTEGER :: unit
1695 INTEGER :: sectounit
1696END TYPE unitchecker
1697
1698TYPE(unitchecker),PARAMETER :: hunit(5) = (/ &
1699 unitchecker(1, 3600), unitchecker(10, 10800), unitchecker(11, 21600), &
1700 unitchecker(12, 43200), unitchecker(2, 86400) /)
1701TYPE(unitchecker),PARAMETER :: munit(3) = (/ & ! 13 14 COSMO extensions
1702 unitchecker(0, 60), unitchecker(13, 900), unitchecker(14, 1800) /)
1703
1704unit = imiss
1705IF (.NOT.c_e(valuein1) .OR. .NOT.c_e(valuein2)) THEN
1706 valueout1 = imiss
1707 valueout2 = imiss
1708 unit = 1
1709ELSE IF (mod(valuein1,3600) == 0 .AND. mod(valuein2,3600) == 0) THEN ! prefer hours
1710 DO i = 1, SIZE(hunit)
1711 IF (mod(valuein1, hunit(i)%sectounit) == 0 &
1712 .AND. mod(valuein2, hunit(i)%sectounit) == 0 &
1713 .AND. valuein1/hunit(i)%sectounit < 255 &
1714 .AND. valuein2/hunit(i)%sectounit < 255) THEN
1715 valueout1 = valuein1/hunit(i)%sectounit
1716 valueout2 = valuein2/hunit(i)%sectounit
1717 unit = hunit(i)%unit
1718 EXIT
1719 ENDIF
1720 ENDDO
1721 IF (.NOT.c_e(unit)) THEN
1722! last chance, disable overflow check and start from longest unit,
1723 DO i = SIZE(hunit), 1, -1
1724 IF (mod(valuein1, hunit(i)%sectounit) == 0 &
1725 .AND. mod(valuein2, hunit(i)%sectounit) == 0) THEN
1726 valueout1 = valuein1/hunit(i)%sectounit
1727 valueout2 = valuein2/hunit(i)%sectounit
1728 unit = hunit(i)%unit
1729 EXIT
1730 ENDIF
1731 ENDDO
1732 ENDIF
1733ELSE IF (mod(valuein1,60) == 0. .AND. mod(valuein2,60) == 0) THEN ! then minutes
1734 DO i = 1, SIZE(munit)
1735 IF (mod(valuein1, munit(i)%sectounit) == 0 &
1736 .AND. mod(valuein2, munit(i)%sectounit) == 0 &
1737 .AND. valuein1/munit(i)%sectounit < 255 &
1738 .AND. valuein2/munit(i)%sectounit < 255) THEN
1739 valueout1 = valuein1/munit(i)%sectounit
1740 valueout2 = valuein2/munit(i)%sectounit
1741 unit = munit(i)%unit
1742 EXIT
1743 ENDIF
1744 ENDDO
1745 IF (.NOT.c_e(unit)) THEN
1746! last chance, disable overflow check and start from longest unit,
1747 DO i = SIZE(munit), 1, -1
1748 IF (mod(valuein1, munit(i)%sectounit) == 0 &
1749 .AND. mod(valuein2, munit(i)%sectounit) == 0) THEN
1750 valueout1 = valuein1/munit(i)%sectounit
1751 valueout2 = valuein2/munit(i)%sectounit
1752 unit = munit(i)%unit
1753 EXIT
1754 ENDIF
1755 ENDDO
1756 ENDIF
1757ENDIF
1758
1759IF (.NOT.c_e(unit)) THEN
1760 CALL l4f_log(l4f_error,'timerange_second_to_g1: cannot find a grib1 timerange unit for coding ' &
1761 //t2c(valuein1)//','//t2c(valuein2)//'s intervals' )
1762 CALL raise_error()
1763ENDIF
1764
1765END SUBROUTINE timerange_choose_unit_g1
1766
1767
1768! Standardize variables and timerange in DB-all.e thinking:
1769!
1770! Timerange 205 (point in time interval) is converted to maximum or
1771! minimum if parameter is recognized as such and parameter is
1772! corrected as well, otherwise 205 is kept (with possible error
1773! conditions later).
1774!
1775! Timerange 257 (COSMO nudging) is converted to point in time if
1776! interval length is 0, or to a proper timerange if parameter is
1777! recognized as a COSMO statistically processed parameter (and in case
1778! of maximum or minimum the parameter is corrected as well); if
1779! parameter is not recognized, it is converted to instantaneous with a
1780! warning.
1781SUBROUTINE normalize_gridinfo(this)
1782TYPE(gridinfo_def),intent(inout) :: this
1783
1784IF (this%timerange%timerange == 254) THEN ! instantaneous
1785
1786!tmin
1787 IF (this%var == volgrid6d_var_new(255,2,16,255)) THEN
1788 this%var%number=11
1789 RETURN
1790 ENDIF
1791
1792!tmax
1793 IF (this%var == volgrid6d_var_new(255,2,15,255)) THEN
1794 this%var%number=11
1795 RETURN
1796 ENDIF
1797
1798ELSE IF (this%timerange%timerange == 205) THEN ! point in time interval
1799
1800!tmin
1801 IF (this%var == volgrid6d_var_new(255,2,16,255)) THEN
1802 this%var%number=11
1803 this%timerange%timerange=3
1804 RETURN
1805 ENDIF
1806
1807!tmax
1808 IF (this%var == volgrid6d_var_new(255,2,15,255)) THEN
1809 this%var%number=11
1810 this%timerange%timerange=2
1811 RETURN
1812 ENDIF
1813
1814! it is accepted to keep 187 since it is wind gust, not max wind
1815 IF (this%var%discipline == 255 .AND. &
1816 any(this%var%centre == cosmo_centre)) THEN ! grib1 & COSMO
1817
1818 IF (this%var%category == 201) THEN ! table 201
1819
1820 IF (this%var%number == 187) THEN ! wind gust
1821! this%var%category=2
1822! this%var%number=32
1823 this%timerange%timerange=2
1824 ENDIF
1825 ENDIF
1826 ENDIF
1827
1828ELSE IF (this%timerange%timerange == 257) THEN ! COSMO-nudging
1829
1830 IF (this%timerange%p2 == 0) THEN ! point in time
1831
1832 this%timerange%timerange=254
1833
1834 ELSE ! guess timerange according to parameter
1835
1836 IF (this%var%discipline == 255 .AND. &
1837 any(this%var%centre == cosmo_centre)) THEN ! grib1 & COSMO
1838
1839 IF (this%var%category >= 1 .AND. this%var%category <= 3) THEN ! WMO table 2
1840
1841 if (this%var%number == 11) then ! T
1842 this%timerange%timerange=0 ! average
1843
1844 else if (this%var%number == 15) then ! T
1845 this%timerange%timerange=2 ! maximum
1846 this%var%number=11 ! reset also parameter
1847
1848 else if (this%var%number == 16) then ! T
1849 this%timerange%timerange=3 ! minimum
1850 this%var%number=11 ! reset also parameter
1851
1852 else if (this%var%number == 17) then ! TD
1853 this%timerange%timerange=0 ! average
1854
1855 else if (this%var%number == 33) then ! U
1856 this%timerange%timerange=0 ! average
1857
1858 else if (this%var%number == 34) then ! V
1859 this%timerange%timerange=0 ! average
1860
1861 else if (this%var%number == 57) then ! evaporation
1862 this%timerange%timerange=1 ! accumulated
1863
1864 else if (this%var%number == 61) then ! TOT_PREC
1865 this%timerange%timerange=1 ! accumulated
1866
1867 else if (this%var%number == 78) then ! SNOW_CON
1868 this%timerange%timerange=1 ! accumulated
1869
1870 else if (this%var%number == 79) then ! SNOW_GSP
1871 this%timerange%timerange=1 ! accumulated
1872
1873 else if (this%var%number == 90) then ! RUNOFF
1874 this%timerange%timerange=1 ! accumulated
1875
1876 else if (this%var%number == 111) then ! fluxes
1877 this%timerange%timerange=0 ! average
1878 else if (this%var%number == 112) then
1879 this%timerange%timerange=0 ! average
1880 else if (this%var%number == 113) then
1881 this%timerange%timerange=0 ! average
1882 else if (this%var%number == 114) then
1883 this%timerange%timerange=0 ! average
1884 else if (this%var%number == 121) then
1885 this%timerange%timerange=0 ! average
1886 else if (this%var%number == 122) then
1887 this%timerange%timerange=0 ! average
1888 else if (this%var%number == 124) then
1889 this%timerange%timerange=0 ! average
1890 else if (this%var%number == 125) then
1891 this%timerange%timerange=0 ! average
1892 else if (this%var%number == 126) then
1893 this%timerange%timerange=0 ! average
1894 else if (this%var%number == 127) then
1895 this%timerange%timerange=0 ! average
1896
1897 endif
1898
1899 ELSE IF (this%var%category == 201) THEN ! table 201
1900
1901 if (this%var%number == 5) then ! photosynthetic flux
1902 this%timerange%timerange=0 ! average
1903
1904 else if (this%var%number == 20) then ! SUN_DUR
1905 this%timerange%timerange=1 ! accumulated
1906
1907 else if (this%var%number == 22) then ! radiation fluxes
1908 this%timerange%timerange=0 ! average
1909 else if (this%var%number == 23) then
1910 this%timerange%timerange=0 ! average
1911 else if (this%var%number == 24) then
1912 this%timerange%timerange=0 ! average
1913 else if (this%var%number == 25) then
1914 this%timerange%timerange=0 ! average
1915 else if (this%var%number == 26) then
1916 this%timerange%timerange=0 ! average
1917 else if (this%var%number == 27) then
1918 this%timerange%timerange=0 ! average
1919
1920 else if (this%var%number == 42) then ! water divergence
1921 this%timerange%timerange=1 ! accumulated
1922
1923 else if (this%var%number == 102) then ! RAIN_GSP
1924 this%timerange%timerange=1 ! accumulated
1925
1926 else if (this%var%number == 113) then ! RAIN_CON
1927 this%timerange%timerange=1 ! accumulated
1928
1929 else if (this%var%number == 132) then ! GRAU_GSP
1930 this%timerange%timerange=1 ! accumulated
1931
1932 else if (this%var%number == 135) then ! HAIL_GSP
1933 this%timerange%timerange=1 ! accumulated
1934
1935 else if (this%var%number == 187) then ! UVMAX
1936! this%var%category=2 ! reset also parameter
1937! this%var%number=32
1938 this%timerange%timerange=2 ! maximum
1939
1940 else if (this%var%number == 218) then ! maximum 10m dynamical gust
1941 this%timerange%timerange=2 ! maximum
1942
1943 else if (this%var%number == 219) then ! maximum 10m convective gust
1944 this%timerange%timerange=2 ! maximum
1945
1946 endif
1947
1948 ELSE IF (this%var%category == 202) THEN ! table 202
1949
1950 if (this%var%number == 231) then ! sso parameters
1951 this%timerange%timerange=0 ! average
1952 else if (this%var%number == 232) then
1953 this%timerange%timerange=0 ! average
1954 else if (this%var%number == 233) then
1955 this%timerange%timerange=0 ! average
1956 endif
1957
1958 ELSE ! parameter not recognized, set instantaneous?
1959
1960 call l4f_category_log(this%category,l4f_warn, &
1961 'normalize_gridinfo: found COSMO non instantaneous analysis 13,0,'//&
1962 trim(to_char(this%timerange%p2)))
1963 call l4f_category_log(this%category,l4f_warn, &
1964 'associated to an apparently instantaneous parameter '//&
1965 trim(to_char(this%var%centre))//','//trim(to_char(this%var%category))//','//&
1966 trim(to_char(this%var%number))//','//trim(to_char(this%var%discipline)))
1967 call l4f_category_log(this%category,l4f_warn, 'forcing to instantaneous')
1968
1969 this%timerange%p2 = 0
1970 this%timerange%timerange = 254
1971
1972 ENDIF ! category
1973 ENDIF ! grib1 & COSMO
1974 ENDIF ! p2
1975ENDIF ! timerange
1976
1977IF (this%var%discipline == 255 .AND. &
1978 any(this%var%centre == ecmwf_centre)) THEN ! grib1 & ECMWF
1979! useful references:
1980! http://www.ecmwf.int/publications/manuals/libraries/tables/tables_index.html
1981! http://www.ecmwf.int/products/data/technical/soil/discret_soil_lay.html
1982
1983 IF (this%var%category == 128) THEN ! table 128
1984
1985 IF ((this%var%number == 142 .OR. & ! large scale precipitation
1986 this%var%number == 143 .OR. & ! convective precipitation
1987 this%var%number == 144 .OR. & ! total snow
1988 this%var%number == 228 .OR. & ! total precipitation
1989 this%var%number == 145 .OR. & ! boundary layer dissipation
1990 this%var%number == 146 .OR. & ! surface sensible heat flux
1991 this%var%number == 147 .OR. & ! surface latent heat flux
1992 this%var%number == 169) .AND. & ! surface solar radiation downwards
1993 this%timerange%timerange == 254) THEN
1994 this%timerange%timerange = 1 ! accumulated
1995 this%timerange%p2 = this%timerange%p1 ! length of period = forecast time
1996
1997 ELSE IF ((this%var%number == 165 .OR. & ! 10m U
1998 this%var%number == 166) .AND. & ! 10m V
1999 this%level%level1 == 1) THEN
2000 this%level%level1 = 103
2001 this%level%l1 = 10000 ! 10m
2002
2003 ELSE IF ((this%var%number == 167 .OR. & ! 2m T
2004 this%var%number == 168) .AND. & ! 2m Td
2005 this%level%level1 == 1) THEN
2006 this%level%level1 = 103
2007 this%level%l1 = 2000 ! 2m
2008
2009 ELSE IF (this%var%number == 39 .OR. this%var%number == 139 .OR. this%var%number == 140) THEN ! SWVL1,STL1,SWL1
2010 this%level%level1 = 106 ! below surface
2011 this%level%l1 = 0
2012 this%level%l2 = 70 ! 7cm
2013
2014 ELSE IF (this%var%number == 40 .OR. this%var%number == 170) THEN ! SWVL2,STL2 (STL2 wrong before 2000)
2015 this%level%level1 = 106 ! below surface
2016 this%level%l1 = 70
2017 this%level%l2 = 280
2018
2019 ELSE IF (this%var%number == 171) THEN ! SWL2
2020 this%level%level1 = 106 ! below surface
2021 this%level%l1 = 70
2022 this%level%l2 = 210
2023
2024 ELSE IF (this%var%number == 41 .OR. this%var%number == 183) THEN ! SWVL3,STL3 (STL3 wrong before 2000)
2025 this%level%level1 = 106 ! below surface
2026 this%level%l1 = 280
2027 this%level%l2 = 1000
2028
2029 ELSE IF (this%var%number == 184) THEN ! SWL3
2030 this%level%level1 = 106 ! below surface
2031 this%level%l1 = 210
2032 this%level%l2 = 1000
2033
2034 ELSE IF (this%var%number == 42 .OR. this%var%number == 236 .OR. this%var%number == 237) THEN ! SWVL4,STL4,SWL4
2035 this%level%level1 = 106 ! below surface
2036 this%level%l1 = 1000
2037 this%level%l2 = 2890
2038
2039 ELSE IF (this%var%number == 121 .AND. &
2040 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! MX2T6
2041 this%timerange%timerange = 2 ! max
2042 this%timerange%p2 = 21600 ! length of period = 6 hours
2043 this%var%number=167 ! set to T2m, it could be 130 T as well
2044 this%level%level1 = 103
2045 this%level%l1 = 2000 ! 2m
2046
2047 ELSE IF (this%var%number == 122 .AND. &
2048 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! MN2T6
2049 this%timerange%timerange = 3 ! min
2050 this%timerange%p2 = 21600 ! length of period = 6 hours
2051 this%var%number=1
2052 this%var%number=167 ! set to T2m, it could be 130 T as well
2053 this%level%level1 = 103
2054 this%level%l1 = 2000 ! 2m
2055
2056 ELSE IF (this%var%number == 123 .AND. &
2057 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! 10FG6
2058 this%timerange%timerange = 2 ! max
2059 this%timerange%p2 = 21600 ! length of period = 6 hours
2060 this%level%level1 = 103
2061 this%level%l1 = 10000 ! 10m
2062
2063! set cloud cover to bufr style
2064 ELSE IF (this%var%number == 186) THEN ! low cloud cover
2065 this%var%number = 248
2066 this%level = vol7d_level_new(level1=256, level2=258, l2=1)
2067 ELSE IF (this%var%number == 187) THEN ! medium cloud cover
2068 this%var%number = 248
2069 this%level = vol7d_level_new(level1=256, level2=258, l2=2)
2070 ELSE IF (this%var%number == 188) THEN ! high cloud cover
2071 this%var%number = 248
2072 this%level = vol7d_level_new(level1=256, level2=258, l2=3)
2073
2074 ENDIF
2075 ELSE IF (this%var%category == 228) THEN ! table 228
2076
2077 IF (this%var%number == 24) THEN
2078 this%level%level1 = 4 ! Level of 0C Isotherm
2079 this%level%l1 = 0
2080 this%level%level2 = 255
2081 this%level%l2 = 0
2082
2083 ELSE IF (this%var%number == 26 .AND. &
2084 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! MX2T3
2085 this%timerange%timerange = 2 ! max
2086 this%timerange%p2 = 10800 ! length of period = 3 hours
2087 this%var%category = 128 ! reset to table version 128
2088 this%var%number=167 ! set to T2m, it could be 130 T as well
2089 this%level%level1 = 103
2090 this%level%l1 = 2000 ! 2m
2091
2092 ELSE IF (this%var%number == 27 .AND. &
2093 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! MN2T3
2094 this%timerange%timerange = 3 ! min
2095 this%timerange%p2 = 10800 ! length of period = 3 hours
2096 this%var%category = 128 ! reset to table version 128
2097 this%var%number=167 ! set to T2m, it could be 130 T as well
2098 this%level%level1 = 103
2099 this%level%l1 = 2000 ! 2m
2100
2101 ELSE IF (this%var%number == 28 .AND. &
2102 (this%timerange%timerange == 254 .OR. this%timerange%timerange == 205)) THEN ! 10FG3
2103 this%timerange%timerange = 2 ! max
2104 this%timerange%p2 = 10800 ! length of period = 3 hours
2105 this%level%level1 = 103
2106 this%level%l1 = 10000 ! 10m
2107
2108 ENDIF
2109
2110 ENDIF ! table 128
2111ENDIF ! grib1 & ECMWF
2112
2113IF (this%var%discipline == 255 .AND. &
2114 this%var%category >= 1 .AND. this%var%category <= 3) THEN ! grib1 WMO table 2
2115
2116! set cloud cover to bufr style
2117 IF (this%var%number == 73) THEN ! low cloud cover
2118 this%var%number = 71
2119 this%level = vol7d_level_new(level1=256, level2=258, l2=1)
2120 ELSE IF (this%var%number == 74) THEN ! medium cloud cover
2121 this%var%number = 71
2122 this%level = vol7d_level_new(level1=256, level2=258, l2=2)
2123 ELSE IF (this%var%number == 75) THEN ! high cloud cover
2124 this%var%number = 71
2125 this%level = vol7d_level_new(level1=256, level2=258, l2=3)
2126
2127 ENDIF
2128
2129ENDIF
2130
2131
2132END SUBROUTINE normalize_gridinfo
2133
2134
2135! Destandardize variables and timerange from DB-all.e thinking:
2136!
2137! Parameters having maximum or minimum statistical processing and
2138! having a corresponding extreme parameter in grib table, are
2139! temporarily converted to timerange 205 and to the corresponding
2140! extreme parameter; if parameter is not recognized, the max or min
2141! statistical processing is kept (with possible error conditions
2142! later).
2143SUBROUTINE unnormalize_gridinfo(this)
2144TYPE(gridinfo_def),intent(inout) :: this
2145
2146IF (this%timerange%timerange == 3) THEN ! min
2147
2148 IF (this%var == volgrid6d_var_new(255,2,11,255)) THEN ! tmin
2149 this%var%number=16
2150 this%timerange%timerange=205
2151
2152 ELSE IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2153 IF (this%var == volgrid6d_var_new(this%var%centre,128,167,255)) THEN ! tmin
2154 this%var%number=122
2155 this%timerange%timerange=205
2156
2157 ENDIF
2158 ENDIF
2159ELSE IF (this%timerange%timerange == 2) THEN ! max
2160
2161 IF (this%var == volgrid6d_var_new(255,2,11,255)) THEN ! tmax
2162 this%var%number=15
2163 this%timerange%timerange=205
2164
2165 ELSE IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2166 IF (this%var == volgrid6d_var_new(this%var%centre,128,167,255)) THEN ! tmax
2167 this%var%number=121
2168 this%timerange%timerange=205
2169
2170 ELSE IF(this%var == volgrid6d_var_new(this%var%centre,128,123,255)) THEN ! uvmax
2171 this%timerange%timerange=205
2172
2173 ELSE IF(this%var == volgrid6d_var_new(this%var%centre,228,28,255)) THEN ! uvmax
2174 this%timerange%timerange=205
2175
2176 ENDIF
2177 ELSE IF (any(this%var%centre == cosmo_centre)) THEN ! grib1 & COSMO
2178
2179! wind
2180! it is accepted to keep 187 since it is wind gust, not max wind
2181! IF (this%var == volgrid6d_var_new(255,2,32,255)) THEN
2182! this%var%category=201
2183! this%var%number=187
2184! this%timerange%timerange=205
2185! ENDIF
2186 IF (this%var == volgrid6d_var_new(this%var%centre,201,187,255)) THEN
2187 this%timerange%timerange=205
2188 ENDIF
2189 ENDIF
2190ENDIF
2191
2192! reset cloud cover to grib1 style
2193IF (this%var%discipline == 255 .AND. this%var%category == 2) THEN ! grib1 table 2
2194 IF (this%var%number == 71 .AND. &
2195 this%level%level1 == 256 .AND. this%level%level2 == 258) THEN ! l/m/h cloud cover
2196 IF (this%level%l2 == 1) THEN ! l
2197 this%var%number = 73
2198 ELSE IF (this%level%l2 == 2) THEN ! m
2199 this%var%number = 74
2200 ELSE IF (this%level%l2 == 3) THEN ! h
2201 this%var%number = 75
2202 ENDIF
2203 this%level = vol7d_level_new(level1=1) ! reset to surface
2204 ENDIF
2205ENDIF
2206
2207IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2208! reset cloud cover to grib1 style
2209 IF (this%var%discipline == 255 .AND. this%var%category == 128) THEN ! grib1 table 128
2210 IF ((this%var%number == 248 .OR. this%var%number == 164) .AND. &
2211 this%level%level1 == 256 .AND. this%level%level2 == 258) THEN ! l/m/h cloud cover
2212 IF (this%level%l2 == 1) THEN ! l
2213 this%var%number = 186
2214 ELSE IF (this%level%l2 == 2) THEN ! m
2215 this%var%number = 187
2216 ELSE IF (this%level%l2 == 3) THEN ! h
2217 this%var%number = 188
2218 ENDIF
2219 this%level = vol7d_level_new(level1=1) ! reset to surface
2220 ENDIF
2221 ENDIF
2222ENDIF
2223
2224END SUBROUTINE unnormalize_gridinfo
2225#endif
2226
2227
2228! =========================================
2229! gdal driver specific code
2230! could this be moved to a separate module?
2231! =========================================
2232#ifdef HAVE_LIBGDAL
2233SUBROUTINE gridinfo_import_gdal(this, gdalid)
2234TYPE(gridinfo_def),INTENT(inout) :: this ! gridinfo object
2235TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
2236
2237TYPE(gdaldataseth) :: hds
2238
2239
2240!call time_import_gdal(this%time, gaid)
2241this%time = datetime_new(year=2010, month=1, day=1) ! conventional year
2242
2243!call timerange_import_gdal(this%timerange,gaid)
2244this%timerange = vol7d_timerange_new(254, 0, 0) ! instantaneous data
2245
2246!call level_import_gdal(this%level, gaid)
2247hds = gdalgetbanddataset(gdalid) ! go back to dataset
2248IF (gdalgetrastercount(hds) == 1) THEN ! single level dataset
2249 this%level = vol7d_level_new(1, 0) ! surface
2250ELSE
2251 this%level = vol7d_level_new(105, gdalgetbandnumber(gdalid)) ! hybrid level
2252ENDIF
2253
2254!call var_import_gdal(this%var, gaid)
2255this%var = volgrid6d_var_new(centre=255, category=2, number=8) ! topography height
2256
2257END SUBROUTINE gridinfo_import_gdal
2258#endif
2259
2260
2261END MODULE gridinfo_class
2262
2263
2264!>\example example_vg6d_2.f90
2265!!\brief Programma esempio semplice per la lettura di file grib.
2266!!
2267!! Programma che legge i grib contenuti in un file e li organizza in un vettore di oggetti gridinfo
2268
2269
2270!>\example example_vg6d_4.f90
2271!!\brief Programma esempio semplice per la elaborazione di file grib.
2272!!
2273!! Programma che legge uno ad uno i grid contenuti in un file e li
2274!! elabora producendo un file di output contenente ancora grib
Restituiscono il valore dell'oggetto nella forma desiderata.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Copy an object, creating a fully new instance.
Quick method to append an element to the array.
Clone the object, creating a new independent instance of the object exactly equal to the starting one...
Decode and return the data array from a grid_id object associated to a gridinfo object.
Destructor, it releases every information associated with the object.
Display on standard output a description of the gridinfo object provided.
Encode a data array into a grid_id object associated to a gridinfo object.
Export gridinfo descriptors information into a grid_id object.
Import information from a file or grid_id object into the gridinfo descriptors.
Constructor, it creates a new instance of the object.
Method for inserting elements of the array at a desired position.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Method for removing elements of the array at a desired position.
Emit log message for a category with specific priority.
Apply the conversion function this to values.
Classi per la gestione delle coordinate temporali.
Module for describing geographically referenced regular grids.
This module defines an abstract interface to different drivers for access to files containing gridded...
Class for managing information about a single gridded georeferenced field, typically imported from an...
classe per la gestione del logging
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Classe per la gestione dei livelli verticali in osservazioni meteo e affini.
Classe per la gestione degli intervalli temporali di osservazioni meteo e affini.
Class for managing physical variables in a grib 1/2 fashion.
Class for expressing an absolute time value.
This object completely describes a grid on a geographic projection.
Derived type associated to a file-like object containing many blocks/messages/records/bands of gridde...
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Derived type defining a dynamically extensible array of TYPE(gridinfo_def) elements.
Object describing a single gridded message/band.
Definisce il livello verticale di un'osservazione.
Definisce l'intervallo temporale di un'osservazione meteo.
Class defining a real conversion function between units.
Definition of a physical variable in grib coding style.

Generated with Doxygen.