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

Generated with Doxygen.