MODULE getsst_grib
!
! READ SST and LAND/SEA MASK FROM GRIB FILE
! Biju Thomas on 2016-12-15
!

USE grib_mod
USE params
IMPLICIT NONE
  CONTAINS

  SUBROUTINE getsst_grib2

    INTEGER, PARAMETER :: nvar = 2
    INTEGER, PARAMETER :: lugb=11,lugi=31
    LOGICAL :: unpack=.TRUE.
    CHARACTER(LEN=7) :: fnameg, fnamei
    INTEGER, DIMENSION(nvar) :: ig2_parm_cat, ig2_parm_num, ig2_lev_typ, &
                                ig2_lev_val, ig2_discip_num
    INTEGER, DIMENSION(200) :: jids,jpdt,jgdt
    INTEGER :: jskp,jdisc,jpdtn,jgdtn,krec,iret,ipack
    TYPE(gribfield) :: gfld
    REAL, ALLOCATABLE ::  xlon(:), ylat(:), fld(:,:)
    REAL :: lon1, lon2, lat1, lat2, dlon, dlat

    INTEGER :: nlon, nlat, nlon2
    INTEGER :: i,j, iv, kf, ifa, np
    INTEGER :: input_fhr,input_mm,input_dd,input_hh
    INTEGER :: igoret,iioret

    data ig2_parm_cat  /0, 0 /
    data ig2_parm_num  /0, 0 /
    data ig2_lev_typ   /1, 1 /
    data ig2_lev_val   /0, 0 /
    data ig2_discip_num /0, 2/

    OPEN (10, FILE='gribinfo', STATUS='OLD')
    READ(10,*)input_fhr
    READ(10,*)input_mm
    READ(10,*)input_dd
    READ(10,*)input_hh
    CLOSE(10)

    fnameg = "fort.11"
    fnamei = "fort.99"

    CALL baopenr (lugb,fnameg,igoret)
    CALL baopenr (lugi,fnamei,iioret)

    DO iv = 1, nvar
      gfld%idsect => NULL()
      gfld%local => NULL()
      gfld%list_opt => NULL()
      gfld%igdtmpl => NULL()
      gfld%ipdtmpl => NULL()
      gfld%coord_list => NULL()
      gfld%idrtmpl => NULL()
      gfld%bmap => NULL()
      gfld%fld => NULL()
!      jdisc=0  ! This indicates meteorological products
      jdisc = ig2_discip_num(iv)
      jids=-9999
      jpdtn=0  ! 0 = analysis or forecast; 1 = ens fcst
      jgdtn=0  ! lat/lon grid
      jgdt=-9999
      jpdt=-9999

! Set values of jids to ensure the data we are reading
! in comes from the correct mm, dd, hh cycle...

      jids(7) = input_mm  ! Fcst initial month to search for
      jids(8) = input_dd  ! Fcst initial day to search for
      jids(9) = input_hh  ! Fcst initial cycle to search for

      jpdt(1) = ig2_parm_cat(iv)
      jpdt(2) = ig2_parm_num(iv)
      jpdt(8) = 1  ! Value of 1 indicates fcst hr unit is hours
      jpdt(9) = input_fhr*1  ! Fcst lead time to search for
      jpdt(10) = ig2_lev_typ(iv)
      jpdt(12) = ig2_lev_val(iv) !

      jskp = 0
    
      CALL getgb2(lugb,lugi,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt  &
                    ,unpack,krec,gfld,iret)
      IF ( iret == 0) THEN
        ipack = 40
        IF ( gfld%idrtnum.eq.0 ) then      ! Simple packing
          ipack = 0
         ELSEIF ( gfld%idrtnum.eq.2 ) then  ! Complex packing
           ipack = 2
         ELSEIF ( gfld%idrtnum.eq.3 ) then  ! Complex & spatial packing
           ipack = 31
         ELSEIF ( gfld%idrtnum.eq.40.or.gfld%idrtnum.eq.15 ) then
                ! JPEG 2000 packing
           ipack = 40
         ELSEIF ( gfld%idrtnum.eq.41 ) then  ! PNG packing
           ipack = 41
         ENDIF
         kf = gfld%ngrdpts  ! Number of gridpoints returned from read
         nlon = gfld%igdtmpl(8)
         nlat = gfld%igdtmpl(9)
         nlon2 = nlon/2
         IF ( .NOT. ALLOCATED(fld) ) ALLOCATE( fld(nlon,nlat), stat = ifa)
         IF (ifa /=0) THEN
           PRINT* ,'!!! ERROR in allocating f array in grib2'
           PRINT*, 'STOP ALLOC Error(fld) ', ifa
           STOP
         ENDIF
         IF ( .NOT. ALLOCATED(xlon) ) ALLOCATE( xlon(nlon), stat = ifa)
         IF (ifa /=0) THEN
           PRINT* ,'!!! ERROR in allocating f array in grib2'
           PRINT*, 'STOP ALLOC Error(xlon) ', ifa
           STOP
         ENDIF
         IF ( .NOT. ALLOCATED(ylat) ) ALLOCATE( ylat(nlat), stat = ifa)
         IF (ifa /=0) THEN
           PRINT* ,'!!! ERROR in allocating f array in grib2'
           PRINT*, 'STOP ALLOC Error(ylat) ', ifa
           STOP
         ENDIF
         np = 1
         DO j = 1, nlat
           DO i = 1, nlon2
              fld(i+nlon2,nlat-j+1)  = gfld%fld(np)
              np = np + 1
           ENDDO
           DO i = 1+nlon2, nlon
              fld(i-nlon2,nlat-j+1)  = gfld%fld(np)
              np = np + 1
           ENDDO
         ENDDO
         IF ( iv == 1 ) THEN
           WRITE (74) fld
         ELSE IF ( iv == 2 ) THEN
           WRITE (77) fld
         END IF
       ELSE
         PRINT* ,'!!! iret /= 0   in getgb2 CALL (getsst_grib.f90)'
         PRINT* ,'!!! iret = ', iret
         STOP
       ENDIF
    ENDDO
    lat1 = FLOAT(gfld%igdtmpl(12))*1.E-06
    lat2 = FLOAT(gfld%igdtmpl(15))*1.E-06
    lon1 = FLOAT(gfld%igdtmpl(13))*1.E-06
    lon2 = FLOAT(gfld%igdtmpl(16))*1.E-06
    dlat = FLOAT(gfld%igdtmpl(17))*1.E-06
    dlon = FLOAT(gfld%igdtmpl(18))*1.E-06
    DO j = 1, nlat
      ylat(j) = lat2 + (j-1)*dlat
    ENDDO
    DO i = 1, nlon2
      xlon(i+nlon2) = lon1 + (i-1)*dlon 
    ENDDO
    DO i = 1+nlon2, nlon
      xlon(i-nlon2) = lon1 + (i-1)*dlon -360
    ENDDO

    WRITE(23,'(1X,2I5)')nlon,nlat
    WRITE(23,'(1X,10E10.4)')xlon,ylat
    
  END SUBROUTINE getsst_grib2

END MODULE getsst_grib