libsim  Versione7.1.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 
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
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
267 #ifdef DEBUG
268 call l4f_category_log(this%category,l4f_debug,"displaying gridinfo " )
269 #endif
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
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*,"--------------------------------------------------------------"
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)
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 END SUBROUTINE gridinfo_import
354 
355 
362 SUBROUTINE gridinfo_import_from_file(this, filename, categoryappend)
363 TYPE(arrayof_gridinfo) :: this
364 CHARACTER(len=*),INTENT(in) :: filename
365 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
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  IF (.NOT. 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 
426 SUBROUTINE gridinfo_export(this)
427 TYPE(gridinfo_def),INTENT(inout) :: this
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 
464 SUBROUTINE gridinfo_export_to_file(this, filename, categoryappend)
465 TYPE(arrayof_gridinfo) :: this
466 CHARACTER(len=*),INTENT(in) :: filename
467 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
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 
531 FUNCTION gridinfo_decode_data(this) RESULT(field)
532 TYPE(gridinfo_def),INTENT(in) :: this
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 
547 SUBROUTINE gridinfo_encode_data(this, field)
548 TYPE(gridinfo_def),INTENT(inout) :: this
549 REAL,intent(in) :: field(:,:)
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 INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
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 INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
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 IF (editionnumber == 1 .OR. 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  IF (status == grib_success .AND. 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  ELSE IF (status == grib_success .AND. ttimeincr == 2 .AND. 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  ELSE IF ((status == grib_success .AND. ttimeincr == 2) .OR. &
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
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  IF (g2cosmo_behavior .AND. 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 INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
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  IF (scalef1 == -1 .OR. 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  IF (scalef2 == -1 .OR. 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  IF (.NOT.c_e(scalef1) .OR. .NOT.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  IF (ltype2 == 255 .OR. .NOT.c_e(scalef2) .OR. .NOT.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 INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
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  IF (status == grib_success .AND. statproc >= 0 .AND. 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  IF (ttimeincr == 2 .AND. 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  IF (pdtn < 0 .OR. 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  ELSE IF (this%timerange >= 0 .AND. this%timerange < 254) THEN
929 ! statistically processed -> template 4.8
930  IF (pdtn < 8 .OR. 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 INTEGER,INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
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()
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
1085 
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 IF (ltype == 255 .OR. ltype == -1) THEN
1100  lt = imiss
1101  l = imiss
1102 ELSE IF (ltype <= 10 .OR. ltype == 101 .OR. (ltype >= 162 .AND. ltype <= 184)) THEN
1103  lt = ltype
1104  l = imiss
1105 ELSE
1106  lt = ltype
1107  IF (c_e(scalef) .AND. 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 ELSE IF (lt <= 10 .OR. lt == 101 .OR. (lt >= 162 .AND. lt <= 184)) THEN
1148  ltype = lt
1149  scalev = 0
1150  scalef = 0
1151 ELSE IF (lt == 256 .AND. 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 if (ltype > 0 .and. 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 if (ltype1 > 0 .and. ltype1 <= 9 .and. ltype2 == 255) then ! simple
1344  ltype = ltype1
1345  l1 = 0
1346  l2 = 0
1347 else if (ltype1 == 20 .and. ltype2 == 255) then ! isothermal
1348  ltype = 2
1349  l1 = rescale2(scalef1-2,scalev1)!*100
1350  l2 = 0
1351 else if (ltype1 == 100 .and. ltype2 == 255) then ! isobaric
1352  ltype = 100
1353  l1 = rescale2(scalef1+2,scalev1)!/100
1354  l2 = 0
1355 else if (ltype1 == 100 .and. ltype2 == 100) then
1356  ltype = 101
1357  l1 = rescale1(scalef1+3,scalev1)!/1000
1358  l2 = rescale1(scalef2+3,scalev2)!/1000
1359 else if (ltype1 == 101 .and. ltype2 == 255) then
1360  ltype = 102
1361  l1 = 0
1362  l2 = 0
1363 else if (ltype1 == 102 .and. ltype2 == 255) then ! altitude over sea level
1364  ltype = 103
1365  l1 = rescale2(scalef1,scalev1)
1366  l2 = 0
1367 else if (ltype1 == 102 .and. ltype2 == 102) then
1368  ltype = 104
1369  l1 = rescale1(scalef1+2,scalev1)!/100
1370  l2 = rescale1(scalef2+2,scalev2)!/100
1371 else if (ltype1 == 103 .and. ltype2 == 255) then ! height over ground
1372  ltype = 105
1373  l1 = rescale2(scalef1,scalev1)
1374  l2 = 0
1375 else if (ltype1 == 103 .and. ltype2 == 103) then
1376  ltype = 106
1377  l1 = rescale1(scalef1+2,scalev1)!/100
1378  l2 = rescale1(scalef2+2,scalev2)!/100
1379 else if (ltype1 == 104 .and. ltype2 == 255) then ! sigma
1380  ltype = 107
1381  l1 = rescale2(scalef1,scalev1-4)!*10000
1382  l2 = 0
1383 else if (ltype1 == 104 .and. ltype2 == 104) then
1384  ltype = 108
1385  l1 = rescale1(scalef1-2,scalev1)!*100
1386  l2 = rescale1(scalef2-2,scalev2)!*100
1387 else if (ltype1 == 105 .and. ltype2 == 255) then ! hybrid
1388  ltype = 109
1389  l1 = rescale2(scalef1,scalev1)
1390  l2 = 0
1391 else if (ltype1 == 105 .and. ltype2 == 105) then
1392  ltype = 110
1393  l1 = rescale1(scalef1,scalev1)
1394  l2 = rescale1(scalef2,scalev2)
1395 else if (ltype1 == 106 .and. ltype2 == 255) then ! depth
1396  ltype = 111
1397  l1 = rescale2(scalef1-2,scalev1)!*100
1398  l2 = 0
1399 else if (ltype1 == 106 .and. ltype2 == 106) then
1400  ltype = 112
1401  l1 = rescale1(scalef1-2,scalev1)!*100
1402  l2 = rescale1(scalef2-2,scalev2)!*100
1403 else if (ltype1 == 107 .and. ltype2 == 255) then ! isentropic
1404  ltype = 113
1405  l1 = rescale2(scalef1,scalev1)
1406  l2 = 0
1407 else if (ltype1 == 107 .and. ltype2 == 107) then
1408  ltype = 114
1409  l1 = rescale1(scalef1,scalev1)
1410  l2 = rescale1(scalef2,scalev2)
1411 else if (ltype1 == 108 .and. ltype2 == 255) then ! pressure diff to ground
1412  ltype = 115
1413  l1 = rescale2(scalef1+2,scalev1)!/100
1414  l2 = 0
1415 else if (ltype1 == 108 .and. ltype2 == 108) then
1416  ltype = 116
1417  l1 = rescale1(scalef1+2,scalev1)!/100
1418  l2 = rescale1(scalef2+2,scalev2)!/100
1419 else if (ltype1 == 109 .and. ltype2 == 255) then ! potential vorticity surf
1420  ltype = 117
1421  l1 = rescale2(scalef1+9,scalev1)!/10**9
1422  l2 = 0
1423 else if (ltype1 == 111 .and. ltype2 == 255) then ! eta level
1424  ltype = 119
1425  l1 = rescale2(scalef1-2,scalev1)!*100
1426  l2 = 0
1427 else if (ltype1 == 111 .and. ltype2 == 111) then
1428  ltype = 120
1429  l1 = rescale1(scalef1-4,scalev1)!*10000
1430  l2 = rescale1(scalef2-4,scalev2)!*10000
1431 else if (ltype1 == 160 .and. ltype2 == 255) then ! depth below sea
1432  ltype = 160
1433  l1 = rescale2(scalef1,scalev1)
1434  l2 = 0
1435 else if (ltype1 == 1 .and. ltype2 == 8) then ! entire atmosphere
1436  ltype = 200
1437 else if (ltype1 == 1 .and. 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 IF (tri == 0 .OR. 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 if (statproc == 254 .and. 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 IF (unit >= lbound(unitlist,1) .AND. 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 IF (unit >= lbound(unitlist,1) .AND. 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 IF (p1_g1 > 255 .OR. 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 IF (.NOT.c_e(valuein1) .OR. .NOT.c_e(valuein2)) THEN
1689  valueout1 = imiss
1690  valueout2 = imiss
1691  unit = 1
1692 ELSE IF (mod(valuein1,3600) == 0 .AND. 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  IF (.NOT.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 ELSE IF (mod(valuein1,60) == 0. .AND. 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  IF (.NOT.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 IF (.NOT.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  IF (this%var%discipline == 255 .AND. &
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  IF (this%var%discipline == 255 .AND. &
1820  any(this%var%centre == cosmo_centre)) THEN ! grib1 & COSMO
1821 
1822  IF (this%var%category == 2) THEN ! 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 IF (this%var%discipline == 255 .AND. &
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  IF ((this%var%number == 142 .OR. & ! large scale precipitation
1969  this%var%number == 143 .OR. & ! convective precipitation
1970  this%var%number == 144 .OR. & ! total snow
1971  this%var%number == 228 .OR. & ! total precipitation
1972  this%var%number == 145 .OR. & ! boundary layer dissipation
1973  this%var%number == 146 .OR. & ! surface sensible heat flux
1974  this%var%number == 147 .OR. & ! surface latent heat flux
1975  this%var%number == 169) .AND. & ! 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  ELSE IF ((this%var%number == 165 .OR. & ! 10m U
1981  this%var%number == 166) .AND. & ! 10m V
1982  this%level%level1 == 1) THEN
1983  this%level%level1 = 103
1984  this%level%l1 = 10000 ! 10m
1985 
1986  ELSE IF ((this%var%number == 167 .OR. & ! 2m T
1987  this%var%number == 168) .AND. & ! 2m Td
1988  this%level%level1 == 1) THEN
1989  this%level%level1 = 103
1990  this%level%l1 = 2000 ! 2m
1991 
1992  ELSE IF (this%var%number == 39 .OR. this%var%number == 139 .OR. 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  ELSE IF (this%var%number == 40 .OR. 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  ELSE IF (this%var%number == 41 .OR. 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  ELSE IF (this%var%number == 42 .OR. this%var%number == 236 .OR. 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  ELSE IF (this%var%number == 121 .AND. &
2023  (this%timerange%timerange == 254 .OR. 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  ELSE IF (this%var%number == 122 .AND. &
2031  (this%timerange%timerange == 254 .OR. 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  ELSE IF (this%var%number == 123 .AND. &
2040  (this%timerange%timerange == 254 .OR. 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  ELSE IF (this%var%number == 26 .AND. &
2067  (this%timerange%timerange == 254 .OR. 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  ELSE IF (this%var%number == 27 .AND. &
2076  (this%timerange%timerange == 254 .OR. 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  ELSE IF (this%var%number == 28 .AND. &
2085  (this%timerange%timerange == 254 .OR. 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 IF (this%var%discipline == 255 .AND. this%var%category == 2) THEN ! grib1 table 2
2097 
2098 ! set cloud cover to bufr style
2099  IF (this%var%number == 73) THEN ! low cloud cover
2100  this%var%number = 71
2101  this%level = vol7d_level_new(level1=256, level2=258, l2=1)
2102  ELSE IF (this%var%number == 74) THEN ! medium cloud cover
2103  this%var%number = 71
2104  this%level = vol7d_level_new(level1=256, level2=258, l2=2)
2105  ELSE IF (this%var%number == 75) THEN ! high cloud cover
2106  this%var%number = 71
2107  this%level = vol7d_level_new(level1=256, level2=258, l2=3)
2108 
2109  ENDIF
2110 
2111 ENDIF
2112 
2113 
2114 END SUBROUTINE normalize_gridinfo
2115 
2116 
2117 ! Destandardize variables and timerange from DB-all.e thinking:
2118 !
2119 ! Parameters having maximum or minimum statistical processing and
2120 ! having a corresponding extreme parameter in grib table, are
2121 ! temporarily converted to timerange 205 and to the corresponding
2122 ! extreme parameter; if parameter is not recognized, the max or min
2123 ! statistical processing is kept (with possible error conditions
2124 ! later).
2125 SUBROUTINE unnormalize_gridinfo(this)
2126 TYPE(gridinfo_def),intent(inout) :: this
2127 
2128 IF (this%timerange%timerange == 3) THEN ! min
2129 
2130  IF (this%var == volgrid6d_var_new(255,2,11,255)) THEN ! tmin
2131  this%var%number=16
2132  this%timerange%timerange=205
2133 
2134  ELSE IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2135  IF (this%var == volgrid6d_var_new(this%var%centre,128,167,255)) THEN ! tmin
2136  this%var%number=122
2137  this%timerange%timerange=205
2138 
2139  ENDIF
2140  ENDIF
2141 ELSE IF (this%timerange%timerange == 2) THEN ! max
2142 
2143  IF (this%var == volgrid6d_var_new(255,2,11,255)) THEN ! tmax
2144  this%var%number=15
2145  this%timerange%timerange=205
2146 
2147  ELSE IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2148  IF (this%var == volgrid6d_var_new(this%var%centre,128,167,255)) THEN ! tmax
2149  this%var%number=121
2150  this%timerange%timerange=205
2151 
2152  ELSE IF(this%var == volgrid6d_var_new(this%var%centre,128,123,255)) THEN ! uvmax
2153  this%timerange%timerange=205
2154 
2155  ELSE IF(this%var == volgrid6d_var_new(this%var%centre,228,28,255)) THEN ! uvmax
2156  this%timerange%timerange=205
2157 
2158  ENDIF
2159  ELSE IF (any(this%var%centre == cosmo_centre)) THEN ! grib1 & COSMO
2160 
2161 ! wind
2162 ! it is accepted to keep 187 since it is wind gust, not max wind
2163 ! IF (this%var == volgrid6d_var_new(255,2,32,255)) THEN
2164 ! this%var%category=201
2165 ! this%var%number=187
2166 ! this%timerange%timerange=205
2167 ! ENDIF
2168  IF (this%var == volgrid6d_var_new(this%var%centre,201,187,255)) THEN
2169  this%timerange%timerange=205
2170  ENDIF
2171  ENDIF
2172 ENDIF
2173 
2174 ! reset cloud cover to grib1 style
2175 IF (this%var%discipline == 255 .AND. this%var%category == 2) THEN ! grib1 table 2
2176  IF (this%var%number == 71 .AND. &
2177  this%level%level1 == 256 .AND. this%level%level2 == 258) THEN ! l/m/h cloud cover
2178  IF (this%level%l2 == 1) THEN ! l
2179  this%var%number = 73
2180  ELSE IF (this%level%l2 == 2) THEN ! m
2181  this%var%number = 74
2182  ELSE IF (this%level%l2 == 3) THEN ! h
2183  this%var%number = 75
2184  ENDIF
2185  this%level = vol7d_level_new(level1=1) ! reset to surface
2186  ENDIF
2187 ENDIF
2188 
2189 IF (any(this%var%centre == ecmwf_centre)) THEN ! ECMWF
2190 ! reset cloud cover to grib1 style
2191  IF (this%var%discipline == 255 .AND. this%var%category == 128) THEN ! grib1 table 128
2192  IF ((this%var%number == 248 .OR. this%var%number == 164) .AND. &
2193  this%level%level1 == 256 .AND. this%level%level2 == 258) THEN ! l/m/h cloud cover
2194  IF (this%level%l2 == 1) THEN ! l
2195  this%var%number = 186
2196  ELSE IF (this%level%l2 == 2) THEN ! m
2197  this%var%number = 187
2198  ELSE IF (this%level%l2 == 3) THEN ! h
2199  this%var%number = 188
2200  ENDIF
2201  this%level = vol7d_level_new(level1=1) ! reset to surface
2202  ENDIF
2203  ENDIF
2204 ENDIF
2205 
2206 END SUBROUTINE unnormalize_gridinfo
2207 #endif
2208 
2209 
2210 ! =========================================
2211 ! gdal driver specific code
2212 ! could this be moved to a separate module?
2213 ! =========================================
2214 #ifdef HAVE_LIBGDAL
2215 SUBROUTINE gridinfo_import_gdal(this, gdalid)
2216 TYPE(gridinfo_def),INTENT(inout) :: this ! gridinfo object
2217 TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
2218 
2219 TYPE(gdaldataseth) :: hds
2220 
2221 
2222 !call time_import_gdal(this%time, gaid)
2223 this%time = datetime_new(year=2010, month=1, day=1) ! conventional year
2224 
2225 !call timerange_import_gdal(this%timerange,gaid)
2226 this%timerange = vol7d_timerange_new(254, 0, 0) ! instantaneous data
2227 
2228 !call level_import_gdal(this%level, gaid)
2229 hds = gdalgetbanddataset(gdalid) ! go back to dataset
2230 IF (gdalgetrastercount(hds) == 1) THEN ! single level dataset
2231  this%level = vol7d_level_new(1, 0) ! surface
2232 ELSE
2233  this%level = vol7d_level_new(105, gdalgetbandnumber(gdalid)) ! hybrid level
2234 ENDIF
2235 
2236 !call var_import_gdal(this%var, gaid)
2237 this%var = volgrid6d_var_new(centre=255, category=2, number=8) ! topography height
2238 
2239 END SUBROUTINE gridinfo_import_gdal
2240 #endif
2241 
2242 
2243 END MODULE gridinfo_class
2244 
2245 
2250 
2251 
Classi per la gestione delle coordinate temporali.
Export gridinfo descriptors information into a grid_id object.
Derived type defining a dynamically extensible array of TYPE(gridinfo_def) elements.
Functions that return a trimmed CHARACTER representation of the input variable.
Method for removing elements of the array at a desired position.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Definition of a physical variable in grib coding style.
This module defines an abstract interface to different drivers for access to files containing gridded...
Derived type associated to a file-like object containing many blocks/messages/records/bands of gridde...
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:234
Decode and return the data array from a grid_id object associated to a gridinfo object.
Encode a data array into a grid_id object associated to a gridinfo object.
Import information from a file or grid_id object into the gridinfo descriptors.
Destructor, it releases every information associated with the object.
This object completely describes a grid on a geographic projection.
Definition: grid_class.F90:267
Classe per la gestione dei livelli verticali in osservazioni meteo e affini.
Quick method to append an element to the array.
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Operatore di resto della divisione.
Classe per la gestione degli intervalli temporali di osservazioni meteo e affini. ...
Copy an object, creating a fully new instance.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Object describing a single gridded message/band.
Class for managing physical variables in a grib 1/2 fashion.
Display on standard output a description of the gridinfo object provided.
Clone the object, creating a new independent instance of the object exactly equal to the starting one...
Constructor, it creates a new instance of the object.
classe per la gestione del logging
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Method for inserting elements of the array at a desired position.
Class defining a real conversion function between units.
Apply the conversion function this to values.
Class for managing information about a single gridded georeferenced field, typically imported from an...
Emit log message for a category with specific priority.

Generated with Doxygen.