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

Generated with Doxygen.