subroutine getsst_nemsio(cin)
      ! This program to read SST and landsea mask from NEMSIO binary format GFS surface
      ! analysis file. The purpose of this program is to mimic the getsst program,
      ! which read SST and landsea mask from sigio foramt GFS surface analysis file,
      ! and write out the SST, landsea mask, as well as lon, lat data in needed order
      ! and format.
      ! 
      ! Note: In GFS NEMSIO binary output file, the longitude ranges from 0. to 360.
      !   and the latitude ranges from 90. to -90. However the POM init needed SST
      !   field should correspond to lon and lat coordinates ranging from -180. to
      !   180. and -90. to 90. respectively. So, after reading the lon, lat
      !   coordinates and the data, they are tranformed into the needed order.
      !    
      ! Sample input file: 
      !   for11 linked to gfs.t00z.sfcanl
      ! Ouput files:
      !   fort.23: lon, lat coordinates in ascii format
      !   fort.74: SST file in binary format
      !   fort.77: landsea mask in binary format
      !
      ! Created by Bin Liu, NOAA/NCEP/EMC, 12/08/2016
      ! 
      use nemsio_module
      implicit none
      !
      type(nemsio_gfile) :: gfile
      integer :: im,jm,nframe,fieldsize,iret
      character(8) :: vname
      character(16) :: vlevtyp
      character(len=*),intent(in) :: cin
      real(4),allocatable :: var1d(:)
      real(4),allocatable :: sst(:,:), landmask(:,:)
      real(4),allocatable :: lon(:,:), lat(:,:)
      
      integer :: i,j
      !
      !--- Initialize
      call nemsio_init(iret=iret)
      !
      !--- Opean a NEMSIO file
      !call getarg(1,cin)
      !cin='for11'
      call nemsio_open(gfile,trim(cin),'READ',iret=iret)
      !
      !--- Get dimension
      call nemsio_getfilehead(gfile,iret=iret,dimx=im,dimy=jm,    
     &                        nframe=nframe)
      !
      !---Allocate array
      fieldsize=(im+2*nframe)*(jm+2*nframe)
      allocate(var1d(fieldsize))
      allocate(sst(im,jm))
      allocate(landmask(im,jm))
      allocate(lon(im,jm))
      allocate(lat(im,jm))
      !
      call nemsio_getfilehead(gfile,lat=var1d,iret=iret)
      ! This is the regular reshape from 1d to 2d array
      !do j=1,jm
      !do i=1,im
      !  lat(i,j)=var1d(i+(j-1)*im)
      !enddo
      !enddo
      ! However, POM needs the array on -180 to 180, and -90 to 90 coordinates
      do j=1,jm
      do i=1,im/2
        lat(im/2+i,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      do i=im/2+1,im
        lat(i-im/2,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      enddo
      
      call nemsio_getfilehead(gfile,lon=var1d,iret=iret)
      do j=1,jm
      do i=1,im/2
        lon(im/2+i,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      do i=im/2+1,im
        lon(i-im/2,jm-j+1)=var1d(i+(j-1)*im)-360.
      enddo
      enddo
      
      print*, 'im,jm=',im,jm
      print*, 'lon(1,1),lon(im/2,jm/2),lon(im,jm)=',   
     &         lon(1,1),lon(im/2,jm/2),lon(im,jm)
      print*, 'lat(1,1),lat(im/2,jm/2),lat(im,jm)=',   
     &         lat(1,1),lat(im/2,jm/2),lat(im,jm)
      print*, 'lon(:,jm/2)=',lon(:,jm/2)
      print*, 'lat(im/2,:)=',lat(im/2,:)
      
      !--- Get one data field out by giving field name, levtyp, and level
      call nemsio_readrecv(gfile,'tmp','sfc',1,var1d,iret=iret)
      do j=1,jm
      do i=1,im/2
        sst(im/2+i,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      do i=im/2+1,im
        sst(i-im/2,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      enddo
      
      print*, 'sst(im/2),jm/2)=',sst(im/2,jm/2)
      
      call nemsio_readrecv(gfile,'land','sfc',1,var1d,iret=iret)
      do j=1,jm
      do i=1,im/2
        landmask(im/2+i,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      do i=im/2+1,im
        landmask(i-im/2,jm-j+1)=var1d(i+(j-1)*im)
      enddo
      enddo
      
      print*, 'landmask(im/2),jm/2)=',landmask(im/2,jm/2)
      
      !write out data
      write(74) sst
      write(77) landmask
      
      write(23,124) im,jm
      write(23,123) lon(:,1),lat(1,:)
      
 124  FORMAT(1X,2I5)
 123  FORMAT(1X,10E10.4)
      
      !
      !--- Close the NEMSIO file
      call nemsio_close(gfile,iret=iret)
      !
      !--- Finalize
      call nemsio_finalize()
      !
      end subroutine getsst_nemsio