SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/SURFEX/prep_grib

Go to the documentation of this file.
00001 !     #########
00002       SUBROUTINE PREP_GRIB_GRID(HGRIB,KLUOUT,HINMODEL,HGRIDTYPE,TPTIME_GRIB)
00003 !     ##########################################################################
00004 !
00005 !!****  *PREP_GRIB_GRID* - reads GRIB grid.
00006 !!
00007 !!    PURPOSE
00008 !!    -------
00009 !!
00010 !!**  METHOD
00011 !!    ------
00012 !!
00013 !!    EXTERNAL
00014 !!    --------
00015 !!
00016 !!    IMPLICIT ARGUMENTS
00017 !!    ------------------
00018 !!
00019 !!
00020 !!    REFERENCE
00021 !!    ---------
00022 !!
00023 !!
00024 !!    AUTHOR
00025 !!    ------
00026 !!
00027 !!      V. Masson (from read_all_data_grib_case)
00028 !!
00029 !!    MODIFICATIONS
00030 !!    -------------
00031 !!      Original   06/2003
00032 !!      S. Faroux 01/2011 : to use library GRIB_API instead of GRIBEX (from
00033 !!                          read_all_data_grib_case)
00034 !-------------------------------------------------------------------------------
00035 !
00036 !*      0. DECLARATIONS
00037 !          ------------
00038 !
00039 USE MODD_TYPE_DATE_SURF
00040 !
00041 USE MODE_READ_GRIB
00042 !
00043 USE MODD_GRID_ROTLATLON
00044 USE MODD_GRID_GAUSS, ONLY : XILA1, XILO1, XILA2, XILO2, NINLA, NINLO, NILEN, LROTPOLE, XCOEF, XLAP, XLOP
00045 USE MODD_GRID_AROME, ONLY : XX, XY, NX, NY, XLAT0, XLON0, XLATOR, XLONOR, XRPK, XBETA
00046 USE MODD_GRID_GRIB,  ONLY : NNI
00047 USE MODD_SURF_PAR,   ONLY : XUNDEF, NUNDEF
00048 USE MODD_CSTS,       ONLY : XPI
00049 !
00050 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00051 USE PARKIND1  ,ONLY : JPRB
00052 !
00053 USE MODI_ABOR1_SFX
00054 !
00055 IMPLICIT NONE
00056 !
00057 !* 0.1. Declaration of arguments
00058 !       ------------------------
00059 !
00060  CHARACTER(LEN=*),  INTENT(IN)   :: HGRIB     ! Grib file name
00061 INTEGER,          INTENT(IN)    :: KLUOUT    ! logical unit of output listing
00062  CHARACTER(LEN=6), INTENT(OUT)   :: HINMODEL  ! Grib originating model
00063  CHARACTER(LEN=10), INTENT(OUT)  :: HGRIDTYPE ! Grid type
00064 TYPE (DATE_TIME)                :: TPTIME_GRIB    ! current date and time
00065 
00066 !
00067 !* 0.2 Declaration of local variables
00068 !      ------------------------------
00069 ! General purpose variables
00070 INTEGER(KIND=kindOfInt)                            :: IRET          ! Return code from subroutines
00071 !
00072 ! Variable involved in the task of reading the grib file
00073 INTEGER(KIND=kindOfInt)                            :: IMISSING
00074 INTEGER(KIND=kindOfInt)                            :: IUNIT
00075 INTEGER(KIND=kindOfInt)                            :: IGRIB
00076 INTEGER                            :: ICENTER       ! number of center
00077  CHARACTER(LEN=20)                  :: HGRID         ! type of grid
00078 INTEGER                            :: ISCAN, JSCAN
00079 INTEGER                            :: ILENX ! nb points in X
00080 INTEGER                            :: ILENY ! nb points in Y
00081 INTEGER                            :: ITIME
00082 INTEGER                            :: IUNITTIME,IP1
00083 !
00084 ! Grib Grid definition variables
00085 INTEGER                            :: JLOOP1        ! Dummy counter
00086 !JUAN
00087 INTEGER(KIND=kindOfInt), DIMENSION(:), ALLOCATABLE :: ININLO_GRIB
00088 !JUAN
00089 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00090 !
00091 !---------------------------------------------------------------------------------------
00092 !
00093 IF (LHOOK) CALL DR_HOOK('PREP_GRIB_GRID',0,ZHOOK_HANDLE)
00094 WRITE (KLUOUT,'(A)') ' -- Grib reader started'
00095 !
00096 ! open grib file
00097  CALL GRIB_OPEN_FILE(IUNIT,HGRIB,'R',IRET)
00098 IF (IRET /= 0) THEN
00099   CALL ABOR1_SFX('PREP_GRIB_GRID: Error opening the grib file '//HGRIB)
00100 END IF
00101 !
00102  CALL GRIB_NEW_FROM_FILE(IUNIT,IGRIB,IRET)
00103 IF (IRET /= 0) THEN
00104   CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading the grib file')
00105 END IF
00106 !
00107 ! close the grib file
00108  CALL GRIB_CLOSE_FILE(IUNIT)
00109 !
00110 !---------------------------------------------------------------------------------------
00111 !* 2.  Fix originating center
00112 !---------------------------------------------------------------------------------------
00113 !
00114  CALL GRIB_GET(IGRIB,'centre',ICENTER,IRET)
00115 IF (IRET /= 0) THEN
00116   CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading center')
00117 END IF
00118 !
00119  CALL GRIB_GET(IGRIB,'typeOfGrid',HGRID,IRET)
00120 IF (IRET /= 0) THEN
00121   CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading type of grid')
00122 END IF
00123 !
00124 SELECT CASE (ICENTER)
00125 
00126   CASE (96)
00127     WRITE (KLUOUT,'(A)') ' | Grib file from HARMONY'  
00128     HINMODEL='ALADIN'
00129     HGRIDTYPE='AROME     '
00130 
00131   CASE (82)
00132     WRITE (KLUOUT,'(A)') ' | Grib file from HIRLAM'  
00133     HINMODEL='HIRLAM'
00134     HGRIDTYPE='ROTLATLON '
00135 
00136   CASE (98)
00137     WRITE (KLUOUT,'(A)') ' | Grib file from European Center for Medium-range Weather Forecast'
00138     HINMODEL = 'ECMWF '
00139     HGRIDTYPE= 'GAUSS     '
00140 
00141   CASE (85)
00142     SELECT CASE (HGRID)  
00143 
00144       CASE('regular_gg')
00145         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Arpege model'
00146         WRITE (KLUOUT,'(A)') 'but same grid as ECMWF model (unstretched)'
00147         HINMODEL = 'ARPEGE'
00148         HGRIDTYPE= 'GAUSS     '
00149 
00150       CASE('reduced_gg')  
00151         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Arpege model'
00152         WRITE (KLUOUT,'(A)') 'but reduced grid'
00153         HINMODEL = 'ARPEGE'
00154         HGRIDTYPE= 'GAUSS     '
00155 
00156       CASE('regular_ll')  
00157         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Mocage model'
00158         HINMODEL = 'MOCAGE'
00159         HGRIDTYPE= 'LATLON    '
00160 
00161       CASE('unknown_PLPresent')
00162         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Arpege model'
00163         HINMODEL = 'ARPEGE'
00164         HGRIDTYPE= 'ROTGAUSS  '
00165 
00166       CASE('lambert')
00167         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Aladin france model'
00168         HINMODEL = 'ALADIN'
00169         HGRIDTYPE= 'AROME     '
00170 
00171       CASE('mercator')
00172         WRITE (KLUOUT,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
00173         HINMODEL = 'ALADIN'
00174         HGRIDTYPE= 'MERCATOR  '
00175 
00176     END SELECT
00177 
00178   CASE DEFAULT
00179     CALL ABOR1_SFX('PREP_GRIB_GRID: GRIB FILE FORMAT NOT SUPPORTED')
00180 
00181 END SELECT
00182 
00183 !---------------------------------------------------------------------------------------
00184 !* 3. Number of points
00185 !---------------------------------------------------------------------------------------
00186 !
00187 NX    = NUNDEF
00188 NY    = NUNDEF
00189 NINLA = NUNDEF
00190 NILEN = NUNDEF
00191 IF (ALLOCATED(NINLO)) DEALLOCATE(NINLO)
00192 !
00193 SELECT CASE (HGRIDTYPE)
00194 
00195      CASE ('AROME    ','MERCATOR  ')
00196      ! 3.1 Lambert conformal projection (ALADIN files)
00197      ! or Mercator projection (ALADIN REUNION files)
00198        CALL GRIB_GET(IGRIB,'Nj',NY,IRET)
00199        CALL GRIB_GET(IGRIB,'Ni',NX,IRET)
00200        NNI= NX * NY
00201      !
00202      !
00203      CASE ('GAUSS    ','ROTGAUSS  ','LATLON    ')
00204      ! 3.2 Usual or Gaussian lat,lon grid (ECMWF files)
00205      !
00206        CALL GRIB_GET(IGRIB,'Nj',NINLA,IRET)
00207        ALLOCATE (NINLO(NINLA))
00208        ALLOCATE (ININLO_GRIB(NINLA))
00209        NILEN = 0
00210        CALL GRIB_IS_MISSING(IGRIB,'pl',IMISSING,IRET)
00211        IF (IRET /= 0 .OR. IMISSING==1)  THEN !  regular
00212          CALL GRIB_GET(IGRIB,'Ni',ININLO_GRIB(1),IRET)  
00213          ININLO_GRIB(2:NINLA)=ININLO_GRIB(1)
00214          NILEN=NINLA*ININLO_GRIB(1)
00215        ELSE !  quasi-regular
00216          CALL GRIB_GET(IGRIB,'pl',ININLO_GRIB)
00217          DO JLOOP1=1 ,NINLA
00218            NILEN = NILEN + ININLO_GRIB(JLOOP1)
00219          ENDDO
00220        ENDIF
00221        NNI = NILEN
00222        NINLO = ININLO_GRIB !JUAN 
00223      CASE ('ROTLATLON')
00224        CALL GRIB_GET(IGRIB,'Nj',NRY,IRET)
00225        CALL GRIB_GET(IGRIB,'Ni',NRX,IRET) 
00226        NNI = NRX * NRY
00227 
00228      CASE DEFAULT
00229        CALL ABOR1_SFX('PREP_GRIB_GRID: GRID PROJECTION NOT SUPPORTED')
00230      !
00231 END SELECT
00232 !
00233 !---------------------------------------------------------------------------------------
00234 !* 4.  Updates grid information
00235 !---------------------------------------------------------------------------------------
00236 !
00237 XX    = XUNDEF
00238 XY    = XUNDEF
00239 XILA1 = XUNDEF
00240 XILO1 = XUNDEF
00241 XILA2 = XUNDEF
00242 XILO2 = XUNDEF
00243 LROTPOLE = .FALSE.
00244 XCOEF = XUNDEF
00245 XLAP  = XUNDEF
00246 XLOP  = XUNDEF
00247    
00248 SELECT CASE (HGRIDTYPE)
00249 
00250      CASE ('AROME     ')  
00251      ! 4.1 Lambert conformal projection (ALADIN files)
00252      !
00253        CALL GRIB_GET(IGRIB,'xDirectionGridLength',ILENX)
00254        CALL GRIB_GET(IGRIB,'yDirectionGridLength',ILENY)  
00255        XY = (NY-1)*ILENY
00256        XX = (NX-1)*ILENX
00257 
00258        CALL GRIB_GET(IGRIB,'Latin1InDegrees',XLAT0) 
00259        CALL GRIB_GET(IGRIB,'LoVInDegrees',XLON0)      
00260        IF (XLON0 > 180.) XLON0 = XLON0 - 360.
00261 
00262        CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XLATOR)
00263        CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XLONOR)
00264        IF (XLONOR > 180.) XLONOR = XLONOR - 360.
00265 
00266        XRPK  = SIN(XLAT0/180.*XPI)
00267        XBETA = 0.
00268      !
00269      CASE ('GAUSS     ','LATLON    ')
00270        HGRIDTYPE = 'GAUSS     '        
00271      ! 4.2 Usual or Gaussian lat,lon grid (ECMWF files)
00272      !     No projection - just stores the grid definition
00273      !
00274        CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XILA1)
00275        CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XILO1)
00276        CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XILA2)
00277        CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XILO2)
00278 
00279        LROTPOLE = .FALSE.
00280      !
00281      CASE ('ROTLATLON ')
00282      !
00283      ! 4.2.5 Rotated lat/lon grid (HIRLAM)
00284      !
00285        CALL GRIB_GET(IGRIB,'iScansNegatively',ISCAN)
00286        CALL GRIB_GET(IGRIB,'jScansNegatively',JSCAN)
00287 
00288        IF (ISCAN+JSCAN == 0 ) THEN
00289          CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XRILA2)  
00290          CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XRILO1)
00291          CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XRILA1)         
00292          CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XRILO2)
00293        ELSEIF (ISCAN+JSCAN == 2) THEN
00294          CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XRILA1)   
00295          CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XRILO2)
00296          CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XRILA2)                 
00297          CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XRILO1)
00298        ELSEIF (ISCAN == 1) THEN
00299          CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XRILA2)   
00300          CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XRILO2)
00301          CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XRILA1)                 
00302          CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XRILO1)
00303        ELSEIF (JSCAN == 1) THEN
00304          CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XRILA1)   
00305          CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XRILO1)
00306          CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XRILA2)                 
00307          CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XRILO2)         
00308        ENDIF
00309 
00310        CALL GRIB_GET(IGRIB,'latitudeOfSouthernPoleInDegrees',XRLAP)                 
00311        CALL GRIB_GET(IGRIB,'longitudeOfSouthernPoleInDegrees',XRLOP)
00312 
00313        CALL GRIB_GET(IGRIB,'iDirectionIncrementInDegrees',XRDX)                 
00314        CALL GRIB_GET(IGRIB,'jDirectionIncrementInDegrees',XRDY)
00315 
00316        WRITE(KLUOUT,*)'XRILA1,XRILO1',XRILA1,XRILO1
00317        WRITE(KLUOUT,*)'XRILA2,XRILO2',XRILA2,XRILO2
00318        WRITE(KLUOUT,*)'XRLAP,XRLOP',XRLAP,XRLOP
00319        WRITE(KLUOUT,*)'XRDX,XRDY',XRDX,XRDY
00320      !
00321      CASE ('ROTGAUSS  ')
00322      ! 4.3 Stretched lat,lon grid (Arpege files)
00323      !
00324        HGRIDTYPE = 'GAUSS     '
00325        CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XILA1)
00326        CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XILO1)
00327        CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XILA2)
00328        CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XILO2)
00329      
00330        LROTPOLE = .TRUE.
00331        CALL GRIB_GET(IGRIB,'stretchingFactor',XCOEF)
00332        CALL GRIB_GET(IGRIB,'latitudeOfStretchingPoleInDegrees',XLAP)
00333        CALL GRIB_GET(IGRIB,'longitudeOfStretchingPoleInDegrees',XLOP)
00334      !
00335      
00336      CASE ('MERCATOR  ')  
00337      ! 4.4 Mercator  projection (ALADIN Reunion files)
00338      !
00339        HGRIDTYPE = 'AROME     '     
00340        CALL GRIB_GET(IGRIB,'Dj',ILENY)     
00341        CALL GRIB_GET(IGRIB,'Di',ILENX)
00342        XY = (NY-1)*ILENY
00343        XX = (NX-1)*ILENX
00344 
00345        CALL GRIB_GET(IGRIB,'LaDInDegrees',XLAT0)       
00346        XLON0 = 0.
00347 
00348        CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XLATOR)
00349        CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XLONOR)
00350        IF (XLONOR > 180.) XLONOR = XLONOR - 360.
00351 
00352        XRPK  = 0.
00353        XBETA = 0.
00354 
00355      CASE DEFAULT
00356        WRITE (KLUOUT,'(A)') 'No such projection implemented in prep_grib_grid ',HGRID
00357        CALL ABOR1_SFX('PREP_GRIB_GRID: UNKNOWN PROJECTION')
00358      !
00359 END SELECT
00360 !---------------------------------------------------------------------------------------
00361 !* 2.4 Read date
00362 !---------------------------------------------------------------------------------------
00363 !
00364 WRITE (KLUOUT,'(A)') ' | Reading date'
00365 !
00366  CALL GRIB_GET(IGRIB,'year',TPTIME_GRIB%TDATE%YEAR,IRET)
00367  CALL GRIB_GET(IGRIB,'month',TPTIME_GRIB%TDATE%MONTH,IRET)
00368  CALL GRIB_GET(IGRIB,'day',TPTIME_GRIB%TDATE%DAY,IRET)
00369  CALL GRIB_GET(IGRIB,'time',ITIME,IRET)
00370 TPTIME_GRIB%TIME=INT(ITIME/100)*3600+(ITIME-INT(ITIME/100)*100)*60
00371 !  
00372  CALL GRIB_GET(IGRIB,'P1',IP1,IRET)
00373 IF ( IP1>0 ) THEN
00374   CALL GRIB_GET(IGRIB,'unitOfTimeRange',IUNITTIME,IRET)      
00375   SELECT CASE (IUNITTIME)       ! Time unit indicator
00376     CASE (1)                    !hour
00377       TPTIME_GRIB%TIME   = TPTIME_GRIB%TIME + IP1*3600.
00378     CASE (0)                    !minute
00379       TPTIME_GRIB%TIME   = TPTIME_GRIB%TIME + IP1*60.
00380   END SELECT
00381 ENDIF
00382 !
00383 !---------------------------------------------------------------------------------------
00384 !
00385  CALL GRIB_RELEASE(IGRIB,IRET)
00386 IF (IRET /= 0) THEN
00387   CALL ABOR1_SFX('PREP_GRIB_GRID: Error in releasing the grib message memory')
00388 END IF
00389 !
00390 IF (LHOOK) CALL DR_HOOK('PREP_GRIB_GRID',1,ZHOOK_HANDLE)
00391 !
00392 END SUBROUTINE PREP_GRIB_GRID