program ex1 
      implicit none

****
*     This program reads a file containing SST data  
*     stored in ASCII format on a lat/lon tabular format 
*     into the RPN standard file format.
*     The data are obtained from the NOAA/NESDIS 14 km SST archive
****

****
*     To run the program, execute "compile_ssttofst" file then   
*     take the binary file "ssttofst" and redirect the input ASCII
*     file to it: e.g. ssttofst < input_text_file
****

****
*     Declare variables used by the RPN standard file library
****
      character*2 nomvar
      character*1 typvar, grtyp
      character*8 etiket
  
****  The ikeys are identifiers for the IPs and IGs        
      integer ikey1, ikey2, ikey3 
      integer key, dateo, deet, npas, ni, nj, nk, npak, datyp 
      integer ip1, ip2, ip3
      integer ig1, ig2, ig3, ig4

****  xlat0 and xlon0 are the location of SW corner of data
****  dlat and dlon are the grid spacing
      real xlat0, xlon0, dlat, dlon
      real xg1, xg2, xg3, xg4

****
*     Declare the name and type of the RPN standard file functions   
****
      external fstecr
      external fnom, fstouv, fclos, fstfrm

      integer fstecr
      integer fnom, fstouv, fclos, fstfrm

****
*     Declare other variables used by the program
****
      integer ier, nrecs
      integer i,j,ii,jj,iun
      integer month
****  The default data domain has 226 points east-west and 361 north-south
      real fld(226, 361), work(226, 361), WORK2
      real fld2(226, 361)
****  1D arrays to hold all the values of lats and lons the data is on
      real lats(361), lons(226)
      
****
*     Association of the RPN standard file produced by the
*     program with the FORTRAN logical unit 1.
****
      iun = 1
      ier = fnom(iun, 'Karen_14km_sst_oct14_cres.fst', 'STD+RND', 0)
      if (ier.lt.0) then
         print *, 'Fatal error while opening the file (FNOM)'
         stop
      endif

****
*     Opening of the standard file
****
      iun = 1
      ier = fstouv(iun, 'RND')
      if (ier.lt.0) then
         print *, 'Cannot open unit:', iun, 
     *            ' in random access mode (FSTOUV)'
         stop
      endif


****
*     Initialization of the standard file attributes that remain
*     constant for all fields
****
      typvar = 'C'
      nomvar = 'TM'
      etiket = 'OCT14_00'

      ip1 = 0
      ip2 = 0
      ip3 = 0

      ni = 226
      nj = 361
      nk = 1

      deet = 0
      npas = 0

      grtyp  = 'L'
      ig1 = 0
      ig2 = 0
      ig3 = 0
      ig4 = 0

      xlat0 = 15.0
      xlon0 = 220.0
      dlat = 0.125 
      dlon = 0.400 

      datyp = 1
      npak = -16

****
*        read in the sst values 
****  
      read(5,*) ((fld(i,j),i=1,226),j=1,361) 
      call land_mask(fld,ni,nj)
****
*        Set a date equal the 1st of each month in 1999
****
      dateo = month * 10000000 + 0199000

** This section in CCs is cut from Serge's routine 
* CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
*
*
*     initialization of igx for a grid type 'L'
*     where the origin is at (0,0) and
*     where the grid length is 1.0 degree in each direction
*
      xg1=0.0
      xg2=0.0
      xg3=1.0
      xg4=1.0
*
      call cxgaig('L',ig1,ig2,ig3,ig4,xg1,xg2,xg3,xg4)

*     --> Building the grid descriptor
*
      ikey1=11
      ikey2=22
      ikey3=33

* Construct the lons and lats column arrays
      call build_latlon(lons,lats,xlon0,xlat0,dlon,dlat,ni,nj)
      
 
      ier=fstecr(lons,WORK2,-32,iun,dateo,deet,npas,ni,1,nk,
     *           ikey1,ikey2,ikey3,typvar,'>>','EAST','L',
     *           ig1,ig2,ig3,ig4,datyp,.false.)

      ier=fstecr(lats,WORK2,-32,iun,dateo,deet,npas,1,nj,nk,
     *           ikey1,ikey2,ikey3,typvar,'^^','NORTH','L',
     *           ig1,ig2,ig3,ig4,datyp,.false.)

*
*     Writing the field in STANDARD file
            
      ip1=0
      ip2=0
      ip3=0

* CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

****
*        Write the TM field to a standard file record
****
      ier = fstecr(fld, WORK, npak, iun, dateo, deet, npas, ni, nj,
     *                 nk, ip1, ip2, ip3, typvar, nomvar, etiket, 'Z',
     *                 ikey1,ikey2,ikey3,0,datyp, .false.)
      
      print*, 'IG1= ',ig1, 'IG2= ',ig2, 'IG3= ',ig3, 'IG4= ',ig4


****
*     Close the standard file
****
      ier = fstfrm(1)

****
*     Unlink the unit 1 from the output.fst file 
****
      ier = fclos(1)

      stop
      end

**** SUBROUTINES

*     This routine is not really necessary for the post cressman data  
      subroutine land_mask(fld,ni,nj)
      integer ni,nj,i,j
      real fld(ni,nj)
      
      do j=1,nj
         do i=1,ni
            if (fld(i,j) .eq. 999.9) fld(i,j)=-1.
*           print*,i,j,fld(i,j)
         enddo
      enddo
      
      return
      end



      subroutine build_latlon(lons,lats,xlon0,xlat0,dlon,dlat,ni,nj)
      real xlat0, xlon0, dlat, dlon 
      integer ni, nj, i, j 
      real lats(nj), lons(ni)

      i=1
      j=1 

      lats(j)=xlat0
      lons(i)=xlon0

      do j=2,nj
         lats(j)=lats(j-1)+dlat
      enddo

      do i=2,ni
         lons(i)=lons(i-1)+dlon
      enddo
          
      return
      end