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