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

Generated with Doxygen.