       program ssmask

!----------------------------------------------------------------------
!   Sea-Salt Emissions Model Version 1.0
!
!   Written: 2/11/03 George Pouliot, NOAA
!   Revised:  4/2/03 per Michelle Mebust recommendations
!   Revised:  1/8/04 by KM Zhang, UC Davis
!            added coastline calculation
!            added freshwater locations in VISTAS 12&36km domains
!   Revised: 5/27/04 by Charles Chang, CSC
!            added freshwater locations in ConUS 32km domain
!   Revised: 8/06/04 by Prakash Bhave, NOAA
!            documented & removed extraneous code
!   Revised: 12/14/04 by Charles Chang, CSC
!            added freshwater locations in BRACE 2&8km domains
!   Revised: 12/22/04 by Prakash Bhave, NOAA
!            corrected MASK_CELL for Northern & Southern coastal cells
!            fixed bug in OPEN calculation for cells with LAND < 10%
!             OLDCODE-- open(c,r) = landfract(c,r)/100.
!            reduced surf-zone width from 100m to 50m
!   Revised: 6/23/05 by George Pouliot, NOAA
!             to read FIPCODE to automatically determine salt water
!   Revised: 7/08/05 by Prakash Bhave, NOAA
!             updated in-line comments; removed obsolete code

!   INPUTS:   LANDTOTAL file, containing a "LAND" field with the
!             percent of each grid cell covered by land and a "FIPCODE"
!             field that identifies Pacific and Atlantic Ocean grid cells
!   OUTPUTS:  OCEANfile, containing 3 fields
!     MASK = 0 for land or freshwater
!            1 for coastline
!            2 for sea water
!     SURF = fraction of grid cell in the surf zone
!     OPEN = fraction of grid cell that is open ocean
!            mask_cell can be used to identify coastline cells
!----------------------------------------------------------------------

c *** Variable declarations
      implicit none
      
      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

      REAL    :: XCELL, YCELL

      INTEGER :: NCOLS, NROWS
      integer    c, r
      integer    cc , rr, x, y, nx, ny, sx, sy
        
      REAL, ALLOCATABLE :: MASK_CELL( :,: )
      INTEGER, ALLOCATABLE :: FIPS(: ,: )
      REAL, ALLOCATABLE :: surf( :,: )
      REAL, ALLOCATABLE :: open( :,: )
      INTEGER, ALLOCATABLE :: ocean( :,: )
      real, allocatable :: landfract(:,:)
      REAL :: AREA
      real :: sz_wid

      CHARACTER(LEN=16) :: progname
      CHARACTER(LEN=16) :: met, land, szwidth, vname, outfile
      CHARACTER(LEN=16) :: outfile2    

      LOGICAL :: USE_10MWIND

      CHARACTER(LEN=80) :: MESG

      INTEGER :: istatus         

      INTEGER :: TRIMLEN
      real       envreal

      EXTERNAL TRIMLEN, envreal

c *** Initialize variables
      progname = 'SSMASK'
      land     = 'LAND_TOTALS'
      
      vname    = 'LAND'
      szwidth  = 'SZ_WID'
      outfile  = 'OUTFILE'
      outfile2 = 'OUTFILE2'


      sz_wid = envreal(szwidth,'surf zone width',50.,istatus)

      if ( .not. open3( land, FSREAD3, progname ) ) THEN
         MESG = 'Could not open file "' //
     &   land( 1: TRIMLEN(land))
     &   // '" for input'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

      IF ( .NOT. DESC3(land))THEN
         MESG = 'Could not get description info for file "' //
     &            land( 1: TRIMLEN( land) ) //'"'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      ENDIF

      ncols = NCOLS3D
      nrows = NROWS3D
      xcell = XCELL3D
      ycell = YCELL3D

      area = xcell*ycell

      allocate ( MASK_CELL(ncols, nrows) )
      allocate ( surf(ncols, nrows) )
      allocate ( open(ncols, nrows) )
      allocate ( landfract(ncols, nrows))
      allocate ( ocean(nx, ny) )      

      allocate (fips(ncols, nrows))



      if ( .not.
     &       read3(land,'OPEN',1,0,0,open)
     &       ) then
             mesg = 'Error reading '//'FIPS'//'from file '//
     &         land( 1: TRIMLEN( land ) )
             call m3exit( progname, 0, 0, MESG, 2 )
      end if


c ------------------ Window to smaller domain with same resolution ------%


       open(40,file=outfile2,status='NEW',form='FORMATTED')

c       write(*,*) nx, ny, sx, sy

       ocean = 0.
       do x = 1, nx
        do y = 1, ny
         cc = x + (sx - 1)
         rr = y + (sy - 1)
         ocean(x,y) = int(open(cc,rr))
        enddo
       enddo

c        do x = 1, ncols
c          write(40,1088) ((int(open(x,y))),y=1,nrows)
        do y = nrows,1,-1
         write(40,1088) ((int(open(x,y))),x=1,ncols) 
c        write(*,*) x,((int(open(x,y))),y=1,nrows)
        enddo

        close(40)

1088    format(500(i1))

c ----------------------------- Write output file -------------------------%
      NVARS3D = 3
      
      VNAME3D(1) = 'MASK'
      UNITS3D(1) = 'none'
      VTYPE3D(1) = M3REAL
      VDESC3D(1) = '2=open ocean, 1=coastline, 0=other'
      
      VNAME3D(2) = 'SURF'
      UNITS3D(2) = 'none'
      VTYPE3D(2) = M3REAL
      VDESC3D(2) = 'surf zone area / total area'
      
      VNAME3D(3) = 'OPEN'
      UNITS3D(3) = 'none'
      VTYPE3D(3) = M3REAL
      VDESC3D(3) = 'open ocean area / total area'

      if ( .not. open3( outfile, FSUNKN3, progname ) ) THEN
          MESG = 'Could not open file "' //
     &     outfile( 1: TRIMLEN(outfile))
     &     // '" for output'
           CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

      if ( .not.
     &        write3(outfile,'MASK',0,0,mask_cell(1,1))
     &        ) then
              mesg = 'Error writing '//'MASK'//'from file '//
     &          outfile( 1: TRIMLEN( outfile ) )
              call m3exit( progname, 0, 0, MESG, 2 )
      end if

      if ( .not.
     &        write3(outfile,'SURF',0,0,surf(1,1))
     &        ) then
              mesg = 'Error writing '//'SURF'//'from file '//
     &          outfile( 1: TRIMLEN( outfile ) )
              call m3exit( progname, 0, 0, MESG, 2 )
      end if

      if ( .not.
     &        write3(outfile,'OPEN',0,0,open(1,1))
     &        ) then
              mesg = 'Error writing '//'OPEN'//'from file '//
     &          outfile( 1: TRIMLEN( outfile ) )
              call m3exit( progname, 0, 0, MESG, 2 )
      end if

      END
