PROGRAM process_NESDIS_imssnow
!
!   PRGMMR: Ming Hu          ORG: GSD        DATE: 2009-04-15
!
! ABSTRACT: 
!     This routine read in NESDIS NESDIS SNOW/ICE data from a grib file and  
!     map them into RR mass grid
!
! 
! PROGRAM HISTORY LOG:
!
!   variable list
!
! USAGE:
!   INPUT FILES:  imssnow
!
!   OUTPUT FILES:  RRimssnow
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 + EXTENSIONS
!   MACHINE:  wJET
!
!$$$
!
!_____________________________________________________________________

  use mpi
  use grib2_read_mod, only : read_grib2_head_dim,read_grib2_head_time
  use grib2_read_mod, only : read_grib2_sngle

  implicit none
! MPI variables
  integer :: npe, mype, mypeLocal,ierror
!
!
!  for grib2
!
  real :: rlatmin,rlonmin
  real*8  :: rdx,rdy
  integer :: nxobs,nyobs
  logical :: fileexist
  integer :: ntot

  integer :: kgds(200)
  character(len=100) :: snowfile
! RR grid
  integer :: nlonRR,nlatRR
!  parameter (nlonRR=1799,nlatRR=1059)
  real,allocatable :: xlonRR(:,:)    !
  real,allocatable :: ylatRR(:,:)    !
  real,allocatable :: xlandRR(:,:)   !
  real,allocatable :: luseRR(:,:)   !
  real,allocatable :: xlandIMS(:,:)   !

  real, allocatable :: snowiceRR(:,:)    ! snow/ice in RR 
!
  integer :: i,j
!
! sea ice and snow, and mask
  REAL, allocatable :: ICEC(:,:)
  REAL, allocatable :: SNOWICEC(:,:)
  Logical*1, allocatable :: lmask(:,:)
  REAL, allocatable :: maskims(:,:)
!  integer, allocatable :: imaskims(:,:)
  integer :: iyear, imonth, iday, ihr, imm

!**********************************************************************
!
!            END OF DECLARATIONS....start of program
! MPI setup
  call MPI_INIT(ierror)
  call MPI_COMM_SIZE(mpi_comm_world,npe,ierror)
  call MPI_COMM_RANK(mpi_comm_world,mype,ierror)

!
  if(mype==0) then

     snowfile='imssnow2'
     inquire(file=trim(snowfile),exist=fileexist)
     if(fileexist) then
        call read_grib2_head_dim(snowfile,nxobs,nyobs,rlonmin,rlatmin,rdx,rdy,kgds)
        write(*,*) 'grib2 file grid =',nxobs,nyobs,rlonmin,rlatmin,rdx,rdy
        write(*,*) 'kgds=',kgds(1:11)
!
! if grib1 file
!    nxobs=6144
!    nyobs=6144
!
        allocate(ICEC(nxobs,nyobs))
        allocate(SNOWICEC(nxobs,nyobs))
        allocate(lmask(nxobs,nyobs))
        allocate(maskims(nxobs,nyobs))
        ICEC=0
        SNOWICEC=0

        call read_grib2_head_time(snowfile,iyear,imonth,iday,ihr,imm)
        write(*,'(a,5I5)') 'date: iyear,imonth,iday,ihr,imm=',&
                      iyear,imonth,iday,ihr,imm
! read field
        ntot=nxobs*nyobs
        call read_grib2_sngle(snowfile,ntot,ICEC,SNOWICEC,lmask)
        write(*,*) 'read in ice =',maxval(icec),minval(icec)
        write(*,*) 'read in snow =',maxval(SNOWICEC),minval(SNOWICEC)
!
! read from grib1 file
!     call read_issnow(icec,snowicec,lmask,kgds,iyear,imonth,iday,ihr)
!     write(*,*) 'NESDIS SNOW and ICE in ', iyear,imonth,iday,ihr
!
     
!
!  Map NESDIS SNOW/ICE to RR grid: Based on the following RUC routine:
! /whome/rucdev/code/13km/hybpre_code/snohires.f
!
        DO J = 1, nxobs
           DO I = 1, nyobs
              IF (icec(I,J) > 0. .OR. snowicec(I,J) > 0.) THEN
                 snowicec(I,J) = 1.
              ELSE
                 snowicec(I,J) = 0.
              ENDIF

!tgs Logical lmask --> .true. for land, .false. for water
             IF (lmask(I,J)) THEN
                maskims(i,j) = 1.    ! land
             ELSE
                maskims(i,j) = 0.    ! water
             ENDIF
           ENDDO
        ENDDO
        deallocate(icec)
        deallocate(lmask)
!
!
        call GET_DIM_ATT_geo('./geo_em.d01.nc',nlonRR,nlatRR)
        write(*,*) 'grid dimension =',nlonRR,nlatRR
        allocate(xlonRR(nlonRR,nlatRR))
        allocate(ylatRR(nlonRR,nlatRR))
        allocate(xlandRR(nlonRR,nlatRR))
        allocate(luseRR(nlonRR,nlatRR))
        allocate(xlandIMS(nlonRR,nlatRR))
        call GET_RR_GRID(xlonRR,ylatRR,nlonRR,nlatRR,xlandRR,luseRR)

! map to RR grid
!
        allocate(snowiceRR(nlonRR,nlatRR))
        snowiceRR=0  ! pecentage
        call map2RR(snowicec,maskims,kgds,xlandRR,nlonRR,nlatRR,xlonRR,ylatRR,snowiceRR,xlandIMS)
        deallocate(snowicec)
        deallocate(maskims)
        deallocate(xlonRR,ylatRR)
!
!
!  trim snow cover field based on NESDIS snow cover data
!
        call update_SNOWICE_netcdf_mass(snowiceRR, xlandRR, luseRR, nlonRR, nlatRR,xlandIMS)
!
     else
        write(*,*) 'file is not exist: ',trim(snowfile)
        write(*,*) 'stop update snow/ice'
     endif

     write(6,*) "=== RAPHRRR PREPROCCESS SUCCESS ==="
  endif ! mype==0

  call MPI_FINALIZE(ierror)
!

END PROGRAM process_NESDIS_imssnow

SUBROUTINE map2RR(f,maskims,kgds_src,xlandRR,nlonRR,nlatRR,xlonRR,ylatRR,snowiceRR,xlandIMS)
!
!   PRGMMR: Ming Hu          ORG: GSD        DATE: 2009-04-15
!
! ABSTRACT:
!     This routine map NESDIS SNOW/ICE data to RR grid
!
! PROGRAM HISTORY LOG:
!
!   variable list
!
! USAGE:
!   INPUT FILES:  imssnow
!
!   OUTPUT FILES:  RRimssnow
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 + EXTENSIONS
!   MACHINE:  wJET
!
!$$$
!
!_____________________________________________________________________

  use gdswzd_mod, only : gdswzd
  implicit none

  INTEGER JPDS(200),JGDS(200),KPDS_SRC(200),KGDS_SRC(200)
  REAL, intent(in):: F(6144,6144)
  REAL, intent(in):: maskims(6144,6144)
!  grid
  integer, intent(in) :: nlonRR,nlatRR
  real, intent(in):: xlandRR(nlonRR,nlatRR)
  real, intent(in):: xlonRR(nlonRR,nlatRR)    !
  real, intent(in):: ylatRR(nlonRR,nlatRR)    !
!
  real, intent(out):: snowiceRR(nlonRR,nlatRR)    !

!
  real :: xlandIMS(nlonRR,nlatRR)
  REAL :: DTR
  REAL :: XPNMC8,YPNMC8,ENNMC8,ALNMC8,ORIENT8
  REAL :: XPNMCAF,YPNMCAF,ENNMCAF,ALNMCAF,ORIENTAF

  real :: YYLAT(1),XLONG(1), RM, RAD, XX(1),YY(1),X, Y
  integer :: IS,IP1,JS,JP1, IPOINT, JPOINT
!
  integer :: iland,KOUNT
  real    :: XRATIO,YRATIO,AREA11,AREA21,AREA12,AREA22,AREA
  integer :: i,j,k,ifound,LL,JPE,JPB,IPE,IPB,NK,MK
!
!
  integer                   :: nret,nout
  real                      :: dum
  real, parameter           :: undefined_value = -1.0

!
  print *,'NESDIS grid info: kgds_src', kgds_src(1:11)

    nout=0
    DO j=1,nlatRR
    DO i=1,nlonRR
!         call gdswiz(kgds_src,-1,1,undefined_value,X,Y, &
!                     xlonRR(i,j),ylatRR(i,j),nret,0,dum,dum)
         YYLAT(1)=ylatRR(i,j)
         XLONG(1)=xlonRR(i,j)
         call gdswzd(kgds_src,-1,1,undefined_value,XX,YY, &
                     XLONG,YYLAT,nret)
         X=XX(1)
         Y=YY(1)
         if (nret /= 1) then
            nout=nout+1
!            snowiceRR(i,j) = 0.
           cycle
         else
! check if X and Y are defined
	   if(X == undefined_value .or. Y == undefined_value) then
!             snowiceRR(i,j) = 0.0 
              cycle
           endif
	   if(X < 1.00001 .or. Y < 1.00001) then
              cycle
           endif
           IS  = NINT(X)
           IP1 = IS + 1
           JS  = NINT(Y)
           JP1 = JS + 1

           ipoint = nint(X)
           jpoint = nint(Y)

!           if (f(ipoint,jpoint) == 1) then
!             snowiceRR(i,j) = 1.0
!           else
!             snowiceRR(i,j) = 0.0
!           end if

         end if


!  AND ONLY SEA POINTS ARE INTERPOLATED TO SEA POINTS (FOR ICE)
!  (NESDIS/IMS LAND MASK: SEA=0,LAND=1, WHILE THE RR MASK in XLAND IS: SEA=2,
!  LAND=1).
!

      ILAND = 1
      IF( int(xlandRR(i,j)) == 0 ) ILAND=0

!tgs - buggy      IF( int(maskims(i,j)) == ILAND ) THEN
      IF( int(maskims(is,js)) == ILAND ) THEN
        snowiceRR(i,j) = F(IPOINT,JPOINT)
      ELSE
!
!  NEAREST NEIGHBOR NOT SAME SFC TYPE, SO USE ALL 4 SURROUNDING POINTS
!
        KOUNT = 0
!
        XRATIO = X - REAL(IS)
        YRATIO = Y - REAL(JS)
!
        AREA11 = (1.0E0 - XRATIO) * (1.0E0 - YRATIO)
        AREA21 = XRATIO * (1.0E0 - YRATIO)
        AREA12 = (1.0E0 - XRATIO) * YRATIO
        AREA22 = XRATIO * YRATIO
!
!
        IF( int(maskims(IS, JS)) .EQ. ILAND) THEN
           KOUNT  = KOUNT + 1
           AREA   = AREA11
           IPOINT = IS
           JPOINT = JS
        END IF
!
        IF( int(maskims(IS, JP1)) .EQ. ILAND ) THEN
           KOUNT = KOUNT +1
           IF (KOUNT .EQ. 1) THEN
              IPOINT = IS
              JPOINT = JP1
           ELSEIF (AREA12 .GT. AREA) THEN
              AREA   = AREA12
              IPOINT = IS
              JPOINT = JP1
           END IF
        END IF
!
        IF( int(maskims(IP1, JS)) .EQ. ILAND ) THEN
           KOUNT = KOUNT + 1
           IF (KOUNT .EQ. 1) THEN
              AREA   = AREA21
              IPOINT = IP1
              JPOINT = JS
           ELSEIF (AREA21 .GT. AREA) THEN
              AREA   = AREA21
              IPOINT = IP1
              JPOINT = JS
           END IF
        END IF
!
!
        IF( int(maskims(IP1, JP1)) .EQ. ILAND ) THEN
           KOUNT = KOUNT + 1
           IF (KOUNT .EQ. 1) THEN
              AREA   = AREA22
              IPOINT = IP1
              JPOINT = JP1
           ELSEIF (AREA22 .GT. AREA) THEN
              AREA   = AREA22
              IPOINT = IP1
              JPOINT = JP1
           END IF
        END IF
!
!     DETERMINE SNO/ICE USING NEAREST NEIGHBOR WITH SAME SFC TYPE 
!
        IF(KOUNT .GT. 0) THEN
            snowiceRR(i,j) = F(IPOINT,JPOINT)
        ELSE
!
!         NO IMMEDIATELY SURROUNDING POINTS IN THE 6144 X 6144 FIELD OF
!         SNOW/ICE HAVE THE SAME LAND-SEA TYPE AS THE model POINT.  THE

!         model POINT MAY BE SMALL ISLAND OR LAKE OR SMALL BAY OR PENNIN.
!         (INVARIABLY A SMALL LAKE IN ETA GRID)
!         SO EXPAND SEARCH RADIUS AND TAKE FIRST SFC TYPE MATCH
!
            IPOINT = NINT(X)
            JPOINT = NINT(Y)
!
!  Define the frame (no. of grid points) over which to search for
!    a matching land/water type from IMS data for the model gridpoint.
            ifound=0
            DO LL=1,16
              JPE = MIN (6144, JPOINT+LL)
              JPB = MAX (1 , JPOINT-LL)
              IPE = MIN (6144, IPOINT+LL)
              IPB = MAX (1 , IPOINT-LL)
!
              DO NK=IPB,IPE
              DO MK=JPB,JPE
                 IF ( int(maskims(nk,mk)) == ILAND .and. ifound ==0 ) THEN
                    snowiceRR(i,j) = F(NK,MK)
                    ifound=1
                 ENDIF
              ENDDO  ! MK
              ENDDO  ! NK
            ENDDO  ! LL
!
!  NO LAND/SEA MASK MATCHES FOUND, SO 
!     A) NORTH OF 55N, WE ASSIGN SNOW/ICE IRRESPECTIVE OF SFC TYPE
!     B) SOUTH OF 55N, WE KEEP A PRIORI ZERO DEFAULT
!   (THE "B" OPTION BEST FOR WARMER LATS OF U.S., WHERE THIS CONDITION 
!   IS VIRTUALLY ALWAYS A SMALL ETA LAKE WITH NO COUNTERPART WATER 
!   NEARBY IN THE NESDIS/IMS GRID, E.G., SALT LAKE, WHERE WE MUST
!   AVOID GETTING SEA-ICE OWING TO SURROUNDING SNOW COVER)
!
!            IF (YYLAT .GE. 55.0 .and. ifound==0) THEN
            IF ( ifound==0 ) THEN
               snowiceRR(i,j) = F(IPOINT,JPOINT)
            ENDIF
        ENDIF   !  KOUNT .GT. 0
!
      ENDIF
!
!tgs - save NESDIS 4-km land/water mask
           xlandIMS(i,j)=maskims(IPOINT,JPOINT)
!

    ENDDO  ! nlon
    ENDDO  ! nlat
                                                                                                                       
    print*,"- WARNING: NUM of MODEL POINT OUTSIDE NESDIS GRID.",nout

end subroutine map2RR

SUBROUTINE read_issnow(icec,snowc,lb,KGDS,iyear,imonth,iday,ihr)
!
!   PRGMMR: Ming Hu          ORG: GSD        DATE: 2009-04-15
!
! ABSTRACT:
!     This routine read in NESDIS NESDIS SNOW/ICE data from a grib file
!
! PROGRAM HISTORY LOG:
!
!   variable list
!
! USAGE:
!   INPUT 
!   OUTPUT 
!      icec : ice cover
!      snowc: snow cover
!      iyear: Year
!      imonth: month
!      iday:   day
!      ihr:  hour
!   INPUT FILES:  imssnow
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 + EXTENSIONS
!   MACHINE:  wJET
!
!$$$
!
!_____________________________________________________________________


  implicit none

  INTEGER JF, MBUF
  PARAMETER (JF=6144*6144)
!  PARAMETER (JF=1024*1024)
  INTEGER JPDS(200),JGDS(200),KPDS(200),KGDS(200)
  LOGICAL*1 LB(JF)   !water/land mask
  REAL, intent(out)::  ICEC(JF)
  REAL, intent(out)::  SNOWC(JF)
!  PARAMETER(MBUF=256*1024)
  PARAMETER(MBUF=256*6144)
  CHARACTER CBUF(MBUF)

  INTEGER :: LUGB,LUGI, J, IRET, K, KF
  character(200):: FNAME1

  INTEGER :: I, iyear, imonth, iday, ihr

!-----------------------------------------------------------------------
  LUGB = 11
  LUGI = 0
  J = 0

!  FNAME1='/public/data/grids/ncep/snow/grib/latest.SNOW_IMS'
  FNAME1='./imssnow'
  call baopen(LUGB,trim(fname1),IRET)
  if ( IRET .ne. 0) then
      write(6,*) 'bad baopen!!! ', IRET
      STOP
  endif

  JPDS=-1
  JGDS=-1
  JPDS(5) = 91  ! ice cover
  call GETGB(LUGB,LUGI,JF,J,JPDS,JGDS,    &
             KF,K,KPDS,KGDS,LB,ICEC,IRET)
  if (IRET .ne. 0) then
    write(6,*) 'bad getgb ', IRET
    STOP
  endif

  JPDS(5) = 238  ! ice cover
  call GETGB(LUGB,LUGI,JF,J,JPDS,JGDS,    &
             KF,K,KPDS,KGDS,LB,SNOWC,IRET)
  if (IRET .ne. 0) then
    write(6,*) 'bad getgb ', IRET
    STOP
  endif

  iyear=2000 + KPDS(8)
  imonth=KPDS(9)
  iday=KPDS(10)
  ihr=KPDS(11)
!       POLAR STEREOGRAPHIC GRIDS
!      KGDS(i)
!          (2)   - N(I) NR POINTS ALONG LAT CIRCLE
!          (3)   - N(J) NR POINTS ALONG LON CIRCLE
!          (4)   - LA(1) LATITUDE OF ORIGIN
!          (5)   - LO(1) LONGITUDE OF ORIGIN
!          (6)   - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
!          (7)   - LOV GRID ORIENTATION
!          (8)   - DX - X DIRECTION INCREMENT
!          (9)   - DY - Y DIRECTION INCREMENT
!          (10)  - PROJECTION CENTER FLAG
!          (11)  - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28)
!  do j=1,6144,100
!    write(*,*) (lb(((j-1)*6144 + i)+0.2),i=1,6144,50)
!  enddo

   KGDS=0
   KGDS(1)=5
   KGDS(2)  = 6144
   KGDS(3)  = 6144
   KGDS(4)  = -21492
   KGDS(5)  = 235000
   KGDS(6)  = 72 
   KGDS(7)  = 280000
   KGDS(8)  = 4000
   KGDS(9)  = 4000
   KGDS(10) = 0
   KGDS(11) = 64
!        
!
end subroutine 

SUBROUTINE GET_RR_GRID(xlon,ylat,nlon,nlat,xland,luse)
!
!   PRGMMR: Ming Hu          ORG: GSD        DATE: 2009-04-15
!
! ABSTRACT:
!     This routine read in Rapid Refresh grid and land mask
!
! PROGRAM HISTORY LOG:
!
!   variable list
!
! USAGE:
!   INPUT:
!
!   OUTPUT:
!      xlon:  longitude in each grid
!      ylat:  latitude in each grid
!      xland: land mask
!
!   INPUT FILES:  imssnow
!
!   OUTPUT FILES:  RRimssnow
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 + EXTENSIONS
!   MACHINE:  wJET
!
!$$$
!
!_____________________________________________________________________


!
  implicit none
!
  INCLUDE 'netcdf.inc'

!  grid
  integer, intent(in) :: nlon,nlat
  real, intent(out):: xlon(nlon,nlat)    !
  real, intent(out):: ylat(nlon,nlat)    !
  real, intent(out):: xland(nlon,nlat)    !
  real, intent(out):: luse(nlon,nlat)    !
!
  integer :: tnlon,tnlat
!
  integer :: NCID
  CHARACTER*180   geofile
  CHARACTER*180   workPath

!
! set geogrid fle name
!
  workPath='./'
  write(geofile,'(a,a)') trim(workPath), 'geo_em.d01.nc'

  write(*,*) 'geofile', trim(geofile)
  call GET_DIM_ATT_geo(geofile,TNLON,TNLAT)
  if( (TNLON.ne.NLON) .or. (TNLAT.ne.NLAT)) then
    write(6,*) ' unmatched dimension'
    write(*,*) 'NLON,NLAT',NLON,NLAT, 'TNLON,TNLAT',TNLON,TNLAT
    stop 123
  endif 
!
!  get GSI horizontal grid in latitude and longitude
!

  call OPEN_geo(geofile, NCID)
  call GET_geo_sngl_geo(NCID,Nlon,Nlat,ylat,xlon,xland,luse)
  call CLOSE_geo(NCID)


END SUBROUTINE GET_RR_GRID