prep_grib_grid.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE prep_grib_grid(HGRIB,KLUOUT,HINMODEL,HGRIDTYPE,HINTERP_TYPE,TPTIME_GRIB)
!     ##########################################################################
!
!!****  *PREP_GRIB_GRID* - reads GRIB grid.
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!
!!      V. Masson (from read_all_data_grib_case)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   06/2003
!!      S. Faroux 01/2011 : to use library GRIB_API instead of GRIBEX (from
!!                          read_all_data_grib_case)
!-------------------------------------------------------------------------------
!
!*      0. DECLARATIONS
!          ------------
!
USE modd_type_date_surf
!
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc
!
USE mode_read_grib
!
USE modd_horibl, ONLY : lgloblon, lglobs, lglobn, xilo1h, xilo2h, no, ninloh, &
                        xla, xola, xolo, np, xloph
!
USE modd_prep,       ONLY : xlat_out, xlon_out, linterp
!
USE modd_grid_rotlatlon
USE modd_grid_gauss, ONLY : xila1, xilo1, xila2, xilo2, ninla, ninlo, nilen, lrotpole, xcoef, xlap, xlop, &
                            xlat, xlon
USE modd_grid_arome, ONLY : xx, xy, nx, ny, xlat0, xlon0, xlator, xlonor, xrpk, xbeta, xzx, xzy, nix
USE modd_grid_grib,  ONLY : nni, cgrib_file
USE modd_surf_par,   ONLY : xundef, nundef
USE modd_csts,       ONLY : xpi
!
USE mode_gridtype_conf_proj
USE modi_horibl_surf_init
USE modi_horibl_surf_coef
USE modi_arpege_stretch_a
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!* 0.1. Declaration of arguments
!       ------------------------
!
 CHARACTER(LEN=*),  INTENT(IN)   :: HGRIB     ! Grib file name
INTEGER,          INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6), INTENT(OUT)   :: HINMODEL  ! Grib originating model
 CHARACTER(LEN=10), INTENT(OUT)  :: HGRIDTYPE ! Grid type
 CHARACTER(LEN=6),  INTENT(OUT)   :: HINTERP_TYPE ! Grid type
type(date_time)                :: tptime_grib    ! current date and time

!
!* 0.2 Declaration of local variables
!      ------------------------------
!
 CHARACTER(LEN=50)                  :: HGRID         ! type of grid
! General purpose variables
INTEGER(KIND=kindOfInt), DIMENSION(:), ALLOCATABLE :: ININLO_GRIB
INTEGER(KIND=kindOfInt)                            :: IMISSING
INTEGER(KIND=kindOfInt)                            :: IUNIT
INTEGER(KIND=kindOfInt)                            :: IGRIB
INTEGER(KIND=kindOfInt)                            :: IRET          ! Return code from subroutines
!
! Variable involved in the task of reading the grib file
INTEGER                            :: ICENTER       ! number of center
INTEGER                            :: ISCAN, JSCAN
INTEGER                            :: ILENX ! nb points in X
INTEGER                            :: ILENY ! nb points in Y
INTEGER                            :: ITIME, IYEAR, IMONTH, IDAY
INTEGER                            :: IUNITTIME,IP1
INTEGER                            :: INO, IINLA     ! output number of points
REAL :: ZTIME
!
! Grib Grid definition variables
INTEGER                            :: JLOOP1        ! Dummy counter
!JUAN
!JUAN
INTEGER :: INFOMPI, J
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!---------------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PREP_GRIB_GRID_1',0,zhook_handle)
!
IF (nrank==npio) THEN
!
WRITE (kluout,'(A)') ' -- Grib reader started'
!
! open grib file
 CALL grib_open_file(iunit,hgrib,'R',iret)
IF (iret /= 0) THEN
  CALL abor1_sfx('PREP_GRIB_GRID: Error opening the grib file '//hgrib)
END IF
!
 CALL grib_new_from_file(iunit,igrib,iret)
! needed for HIRLAM ROTLATLON (infos in the second record)
 CALL grib_new_from_file(iunit,igrib,iret) 
!
IF (iret /= 0) THEN
  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading the grib file')
END IF
!
! close the grib file
 CALL grib_close_file(iunit)
!
!
!---------------------------------------------------------------------------------------
!* 2.  Fix originating center
!---------------------------------------------------------------------------------------
!
 CALL grib_get(igrib,'centre',icenter,iret)
IF (iret /= 0) THEN
  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading center')
END IF
!
 CALL grib_get(igrib,'typeOfGrid',hgrid,iret)
IF (iret /= 0) THEN
  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading type of grid')
END IF
!
SELECT CASE (icenter)

  CASE (96)
    SELECT CASE (hgrid)

      CASE('rotated_ll')
        WRITE (kluout,'(A)') ' | Grib file from HIRLAM - regular latlon grid '  
        hinmodel='HIRLAM'
        hgridtype='ROTLATLON '

      CASE DEFAULT
        WRITE (kluout,'(A)') ' | Grib file from HARMONY'  
        hinmodel='ALADIN'
        hgridtype='AROME     '

    END SELECT

  CASE (82)
    WRITE (kluout,'(A)') ' | Grib file from HIRLAM'  
    hinmodel='HIRLAM'
    hgridtype='ROTLATLON '

  CASE (98)
    WRITE (kluout,'(A)') ' | Grib file from European Center for Medium-range Weather Forecast'
    hinmodel = 'ECMWF '
    hgridtype= 'GAUSS     '

  CASE (85)
    SELECT CASE (hgrid)  

      CASE('regular_gg')
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
        WRITE (kluout,'(A)') 'but same grid as ECMWF model (unstretched)'
        hinmodel = 'ARPEGE'
        hgridtype= 'GAUSS     '

      CASE('reduced_gg')  
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
        WRITE (kluout,'(A)') 'but reduced grid'
        hinmodel = 'ARPEGE'
        hgridtype= 'GAUSS     '

      CASE('regular_ll')  
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Mocage model'
        hinmodel = 'MOCAGE'
        hgridtype= 'LATLON    '

      CASE('unknown_PLPresent')
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
        hinmodel = 'ARPEGE'
        hgridtype= 'ROTGAUSS  '

      CASE('reduced_stretched_rotated_gg')
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
        WRITE (kluout,'(A)') 'but reduced grid'
        hinmodel = 'ARPEGE'
        hgridtype= 'ROTGAUSS  '

      CASE('lambert')
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Aladin france model'
        hinmodel = 'ALADIN'
        hgridtype= 'AROME     '

      CASE('mercator')
        WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
        hinmodel = 'ALADIN'
        hgridtype= 'MERCATOR  '

    END SELECT

  CASE DEFAULT
    CALL abor1_sfx('PREP_GRIB_GRID: GRIB FILE FORMAT NOT SUPPORTED')

END SELECT
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_bcast(igrib,kind(igrib)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(hinmodel,len(hinmodel),mpi_character,npio,ncomm,infompi)
  CALL mpi_bcast(hgridtype,len(hgridtype),mpi_character,npio,ncomm,infompi)
#endif 
ENDIF
!
!---------------------------------------------------------------------------------------
!* 3. Number of points
!---------------------------------------------------------------------------------------
!
nx    = nundef
ny    = nundef
ninla = nundef
nilen = nundef
IF (ALLOCATED(ninlo)) DEALLOCATE(ninlo)
!
SELECT CASE (hgridtype)

     CASE ('AROME    ','MERCATOR  ')
     ! 3.1 lambert conformal projection(aladin files)
     ! or mercator projection(aladin reunion files)
       IF (nrank==npio) THEN
         CALL grib_get(igrib,'Nj',ny,iret)
         CALL grib_get(igrib,'Ni',nx,iret)
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI               
         CALL mpi_bcast(ny,kind(ny)/4,mpi_integer,npio,ncomm,infompi)
         CALL mpi_bcast(nx,kind(nx)/4,mpi_integer,npio,ncomm,infompi)
#endif         
       ENDIF
       nni= nx * ny
     !
     !
     CASE ('GAUSS     ','ROTGAUSS  ','LATLON    ')
     ! 3.2 usual or gaussian lat,lon grid(ecmwf files)
     !
   
       IF (nrank==npio) CALL grib_get(igrib,'Nj',ninla,iret)
       IF (nproc>1) THEN
#ifdef SFX_MPI               
         CALL mpi_bcast(ninla,kind(ninla)/4,mpi_integer,npio,ncomm,infompi)
#endif         
       ENDIF
       ALLOCATE (ninlo(ninla))
       nilen = 0
       IF (nrank==npio) THEN
         ALLOCATE (ininlo_grib(ninla))
        CALL grib_is_missing(igrib,'pl',imissing,iret)
        IF (iret /= 0 .OR. imissing==1)  THEN !  regular
          CALL grib_get(igrib,'Ni',ininlo_grib(1),iret)  
          ininlo_grib(2:ninla)=ininlo_grib(1)
          nilen=ninla*ininlo_grib(1)
        ELSE !  quasi-regular
          CALL grib_get(igrib,'pl',ininlo_grib)
          DO jloop1=1 ,ninla
            nilen = nilen + ininlo_grib(jloop1)
          ENDDO
        ENDIF
        ninlo = ininlo_grib !JUAN 
        DEALLOCATE(ininlo_grib)
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI               
         CALL mpi_bcast(nilen,kind(nilen)/4,mpi_integer,npio,ncomm,infompi)
         CALL mpi_bcast(ninlo,SIZE(ninlo)*kind(ninlo)/4,mpi_integer,npio,ncomm,infompi)
#endif         
       ENDIF
       nni = nilen
     CASE ('ROTLATLON ')
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'Nj',nry,iret)
       CALL grib_get(igrib,'Ni',nrx,iret)
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI               
         CALL mpi_bcast(nry,kind(nry)/4,mpi_integer,npio,ncomm,infompi)
         CALL mpi_bcast(nrx,kind(nrx)/4,mpi_integer,npio,ncomm,infompi)
#endif         
       ENDIF
       nni = nrx * nry

     CASE DEFAULT
       CALL abor1_sfx('PREP_GRIB_GRID: GRID PROJECTION NOT SUPPORTED')
     !
END SELECT
!
!---------------------------------------------------------------------------------------
!* 4.  Updates grid information
!---------------------------------------------------------------------------------------
!
xx    = xundef
xy    = xundef
xila1 = xundef
xilo1 = xundef
xila2 = xundef
xilo2 = xundef
lrotpole = .false.
xcoef = xundef
xlap  = xundef
xlop  = xundef
   
SELECT CASE (hgridtype)

     CASE ('AROME     ')  
     ! 4.1 lambert conformal projection(aladin files)
     !
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'xDirectionGridLength',ilenx)
       CALL grib_get(igrib,'yDirectionGridLength',ileny)
       xy = (ny-1)*ileny
       xx = (nx-1)*ilenx
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xy,kind(xy)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xx,kind(xx)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

       IF (nrank==npio) THEN
       CALL grib_get(igrib,'Latin1InDegrees',xlat0)
       CALL grib_get(igrib,'LoVInDegrees',xlon0)   
       IF (xlon0 > 180.) xlon0 = xlon0 - 360.
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xlat0,kind(xlat0)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlon0,kind(xlon0)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

       IF (nrank==npio) THEN
       CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xlator)
       CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xlonor)
       IF (xlonor > 180.) xlonor = xlonor - 360.
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xlator,kind(xlator)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlonor,kind(xlonor)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

       xrpk  = sin(xlat0/180.*xpi)
       xbeta = 0.
     !
     CASE ('GAUSS     ','LATLON    ')
       hgridtype = 'GAUSS     '        
     ! 4.2 usual or gaussian lat,lon grid(ecmwf files)
     !     no projection - just stores the grid definition
     !
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xila1)
       CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xilo1)
       CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xila2)
       CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xilo2)
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xila1,kind(xila1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xilo1,kind(xilo1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xila2,kind(xila2)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xilo2,kind(xilo2)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

       lrotpole = .false.
     !
     CASE ('ROTLATLON ')
     !
     ! 4.2.5 rotated lat/lon grid(hirlam)
     !
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'iScansNegatively',iscan)
       CALL grib_get(igrib,'jScansPositively',jscan)

       IF (iscan+jscan == 0 ) THEN !lon (i) positive, lat (j) negative
         CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila2)  
         CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo1)
         CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila1)         
         CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo2)
       ELSEIF (iscan+jscan == 2) THEN ! lon (i) negative, lat (j) positive
         CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila1)   
         CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo2)
         CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila2)                 
         CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo1)
       ELSEIF (iscan == 1) THEN ! lon (i) negative, lat (j) negative)
         CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila2)   
         CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo2)
         CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila1)                 
         CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo1)
       ELSEIF (jscan == 1) THEN ! lon (i) positive, lat (j) positive
         CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila1)   
         CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo1)
         CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila2)                 
         CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo2)         
       ENDIF

       CALL grib_get(igrib,'latitudeOfSouthernPoleInDegrees',xrlap)
       CALL grib_get(igrib,'longitudeOfSouthernPoleInDegrees',xrlop)                       

       CALL grib_get(igrib,'iDirectionIncrementInDegrees',xrdx)                 
       CALL grib_get(igrib,'jDirectionIncrementInDegrees',xrdy)
       WRITE(kluout,*)'XRILA1,XRILO1',xrila1,xrilo1
       WRITE(kluout,*)'XRILA2,XRILO2',xrila2,xrilo2
       WRITE(kluout,*)'XRLAP,XRLOP',xrlap,xrlop
       WRITE(kluout,*)'XRDX,XRDY',xrdx,xrdy
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xrila1,kind(xrila1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrilo1,kind(xrilo1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrila2,kind(xrila2)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrilo2,kind(xrilo2)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrlap,kind(xrlap)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrlop,kind(xrlop)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrdx,kind(xrdx)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xrdy,kind(xrdy)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

     !
     CASE ('ROTGAUSS  ')
     ! 4.3 stretched lat,lon grid(arpege files)
     !
       hgridtype = 'GAUSS     '
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xila1)
       CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xilo1)
       CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xila2)
       CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xilo2)
     
       CALL grib_get(igrib,'stretchingFactor',xcoef)
       CALL grib_get(igrib,'latitudeOfStretchingPoleInDegrees',xlap)
       CALL grib_get(igrib,'longitudeOfStretchingPoleInDegrees',xlop)
      ENDIF
      lrotpole = .true.
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xila1,kind(xila1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xilo1,kind(xilo1)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xila2,kind(xila2)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xilo2,kind(xilo2)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xcoef,kind(xcoef)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlap,kind(xlap)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlop,kind(xlop)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF

     !
     
     CASE ('MERCATOR  ')  
     ! 4.4 mercator  projection(aladin reunion files)
     !
       hgridtype = 'AROME     '     
       IF (nrank==npio) THEN
       CALL grib_get(igrib,'Dj',ileny)     
       CALL grib_get(igrib,'Di',ilenx)
       xy = (ny-1)*ileny
       xx = (nx-1)*ilenx

       CALL grib_get(igrib,'LaDInDegrees',xlat0)       

       CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xlator)
       CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xlonor)
       IF (xlonor > 180.) xlonor = xlonor - 360.
       ENDIF
       IF (nproc>1) THEN
#ifdef SFX_MPI                
         CALL mpi_bcast(xy,kind(xy)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xx,kind(xx)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlat0,kind(xlat0)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlator,kind(xlator)/4,mpi_real,npio,ncomm,infompi)
         CALL mpi_bcast(xlonor,kind(xlonor)/4,mpi_real,npio,ncomm,infompi)
#endif
       ENDIF
       xlon0 = 0.
       xrpk  = 0.
       xbeta = 0.

     CASE DEFAULT
       IF (nrank==npio) WRITE (kluout,'(A)') 'No such projection implemented in prep_grib_grid ',hgrid
       CALL abor1_sfx('PREP_GRIB_GRID: UNKNOWN PROJECTION')
     !
END SELECT
!---------------------------------------------------------------------------------------
!* 2.4 Read date
!---------------------------------------------------------------------------------------
!
WRITE (kluout,'(A)') ' | Reading date'
!
IF (nrank==npio) THEN
 CALL grib_get(igrib,'year',iyear,iret)
 CALL grib_get(igrib,'month',imonth,iret)
 CALL grib_get(igrib,'day',iday,iret)
 CALL grib_get(igrib,'time',itime,iret)
ztime=int(itime/100)*3600+(itime-int(itime/100)*100)*60
!  
 CALL grib_get(igrib,'P1',ip1,iret)
IF ( ip1>0 ) THEN
  CALL grib_get(igrib,'unitOfTimeRange',iunittime,iret)      
  SELECT CASE (iunittime)       ! Time unit indicator
    CASE (1)                    !hour
      ztime   = ztime + ip1*3600.
    CASE (0)                    !minute
      ztime   = ztime + ip1*60.
  END SELECT
ENDIF
ENDIF
!
IF (nproc>1) THEN
#ifdef SFX_MPI 
  CALL mpi_bcast(iyear,kind(iyear)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(imonth,kind(imonth)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(iday,kind(iday)/4,mpi_integer,npio,ncomm,infompi)
  CALL mpi_bcast(ztime,kind(ztime)/4,mpi_real,npio,ncomm,infompi)
#endif
ENDIF

tptime_grib%TDATE%YEAR = iyear
tptime_grib%TDATE%MONTH = imonth
tptime_grib%TDATE%DAY = iday
tptime_grib%TIME = ztime
!
!---------------------------------------------------------------------------------------
!
IF (nrank==npio) THEN
 CALL grib_release(igrib,iret)
IF (iret /= 0) THEN
  CALL abor1_sfx('PREP_GRIB_GRID: Error in releasing the grib message memory')
END IF
ENDIF
!
IF (lhook) CALL dr_hook('PREP_GRIB_GRID_1',1,zhook_handle)
IF (lhook) CALL dr_hook('PREP_GRIB_GRID_2',0,zhook_handle)
!
IF (ALLOCATED(xlat_out)) THEN

   ino = SIZE(xlat_out)

  IF (hgridtype=='GAUSS     ') THEN
    IF (ALLOCATED(xlat)) DEALLOCATE(xlat)
    IF (ALLOCATED(xlon)) DEALLOCATE(xlon)
    ALLOCATE(xlat(ino))
    ALLOCATE(xlon(ino))          
    IF (lrotpole) THEN
!* transformation of output latitudes, longitudes into rotated coordinates
      CALL arpege_stretch_a(ino,xlap,xlop,xcoef,xlat_out,xlon_out,xlat,xlon)
    ELSE
      xlat = xlat_out 
      xlon = xlon_out 
    END IF
  ELSEIF (hgridtype=='AROME     ') THEN
    IF (ALLOCATED(xzx)) DEALLOCATE(xzx)
    IF (ALLOCATED(xzy)) DEALLOCATE(xzy)
    IF (ALLOCATED(nix)) DEALLOCATE(nix)
    ALLOCATE(xzx(ino))
    ALLOCATE(xzy(ino))
    ALLOCATE(nix(ny))
    nix=nx          
    CALL xy_conf_proj(xlat0,xlon0,xrpk,xbeta,xlator,xlonor,xzx,xzy,xlat_out,xlon_out)
  ELSEIF (hgridtype=='ROTLATLON ') THEN
  ENDIF
!
  IF (ALLOCATED(no)) DEALLOCATE(no)
  IF (ALLOCATED(xla)) DEALLOCATE(xla)
  IF (ALLOCATED(xola)) DEALLOCATE(xola)
  IF (ALLOCATED(xolo)) DEALLOCATE(xolo)
  IF (ALLOCATED(ninloh)) DEALLOCATE(ninloh)

  ALLOCATE(no(ino,4))
  ALLOCATE(xola(ino),xolo(ino))
  ALLOCATE(xla(ino,4))
!
  
  !
  IF (hgridtype=='GAUSS     ') THEN
    iinla = ninla     
    ALLOCATE(ninloh(iinla+4)) 
    CALL horibl_surf_init(xila1,xilo1,xila2,xilo2,ninla,ninlo,ino,xlon,xlat, &
                          linterp,lgloblon,lglobn,lglobs,no, &
                          ninloh,xola,xolo,xilo1h,xilo2h,xla)
  ELSEIF (hgridtype=='AROME     ') THEN
    iinla = ny
    ALLOCATE(ninloh(iinla+4))
    CALL horibl_surf_init(0.,0.,xy,xx,ny,nix,ino,xzx,xzy, &
                          linterp,lgloblon,lglobn,lglobs,no, &
                          ninloh,xola,xolo,xilo1h,xilo2h,xla)        
  ENDIF              
!
  IF (ALLOCATED(np)) DEALLOCATE(np)
  IF (ALLOCATED(xloph)) DEALLOCATE(xloph)
  ALLOCATE(np(ino,12))
  ALLOCATE(xloph(ino,12))

  IF (lglobs) iinla = iinla + 2
  IF (lglobn) iinla = iinla + 2
  
  IF (hgridtype=='GAUSS     '.OR.hgridtype=='AROME     ') THEN
    CALL horibl_surf_coef(ino,linterp,lgloblon,xilo1h,xilo2h,xolo,&
                          no,ninloh(1:iinla),np,xloph)
  ENDIF
!
ENDIF
!
hinterp_type = "HORIBL"
!
IF (lhook) CALL dr_hook('PREP_GRIB_GRID_2',1,zhook_handle)
!
END SUBROUTINE prep_grib_grid