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

Generated with Doxygen.