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