      SUBROUTINE SUBSURF(iyear,iyears,iyeare,months,monthe,
     *                   idays,idaye,lons,lone,lats,late,
     *                   depth,temp,nlevs,nlevout,nparams,ireqpars,
     *                   ip2,latindex,lonindex,itimeindex,xgrid,
     *                   imple)

C     SUBSURF EXTRACTS ONLY THAT PORTION OF A SURFACE CRUISE
C     (STATION) WHICH A USER REQUESTS

c     imple - set to 1 to implement change in temp array

      parameter (maxlevel=6000, maxcalc=200,bmiss=-1E10)
      parameter (ndays=366,idim0=360,jdim0=180)

      parameter (masknum=100,idimax=1440,jdimax=720)

      dimension depth(maxlevel),temp(maxlevel,maxcalc)
      dimension ip2(0:maxlevel)
      dimension monthend(12)
      data monthend/31,29,31,30,31,30,31,31,30,31,31,31/

      common /geog/ m1(idimax,jdimax), mtot(masknum), mex(idimax,jdimax)
      common /geopnt/ igpnt,xlon,xlat,xdist

      j1=0
      j2=0
      i1=0
      i2=0
      i3=0
      i4=0

      idim=float(idim0)/xgrid
      jdim=float(jdim0)/xgrid

c     Figure julian days to look for

      iy1=iyear
      iy2=iyear
      idatecheck=1

      if ( iyears .eq. 0 .and. iyeare .eq. 0 .and.
     *     months .eq. 0 .and. months .eq. 0 .and.
     *     idays .eq. 0 .and. idaye .eq. 0 ) idatecheck=0

      if ( idatecheck .eq. 1 ) then
       if ( iyears .gt. 0 .and. iyeare .gt. 0 ) then

        iy1=iyears
        iy2=iyeare

       endif

c     set months

       if ( months .gt. 0 .and. monthe .gt. 0  ) then
        im1=months
        im2=monthe
       else
        im1=1
        im2=12
       endif 

c     set days

       if ( idays .gt. 0 .and. idaye .gt. 0 ) then
        id1=idays
        id2=idaye
       else
        id1=1
        id2=monthend(im2)
        if ( im2 .eq. 2 .and. (iy2/4)*4 .ne. iy2 ) id2=id2-1
       endif

c     calculate beginning and ending julian days

       call juliant(iyear,iy1,im1,id1,0.0,0,
     *              xjdays,jsigs)
       call juliant(iyear,iy2,im2,id2,0.0,0,
     *              xjdaye,jsigs)
       xjdaye=xjdaye+1.
      endif

c    Set latitudes, longitudes

      iposcheck=1
      if ( lats .eq. 0 .and. late .eq. 0 .and. lons .eq. 0 .and.
     *     lone .eq. 0 .and. igpnt .eq. 0 ) iposcheck=0

      if ( iposcheck .gt. 0 .and. lons .gt. 0 ) then

       if ( lats .gt. 0 .and. late .gt. 0 ) then
        j1=lats
        j2=late
       else
        j1=1
        j2=jdim
       endif

       if ( lons .gt. 0 .and. lone .gt. 0 ) then
        if ( lons .le. lone ) then
         i1=lons
         i2=lone
         i3=0
         i4=0
        else
         i1=1
         i2=lone
         i3=lons
         i4=idim
        endif

       endif

      endif
         
       
c    Remove levels not within user requested positions/dates

      nend=nlevs
c     nlevs=0
      nlevout=0

      do 50 n=1,nend

       iyes=0

       if ( n .le. maxlevel ) then
        d1=depth(n)
        ic=n
        np=0
       else
        call overmax(maxlevel,nparams,n,0,ic,np)
        d1=temp(ic,np)
       endif

c     check position

       if ( iposcheck .eq. 1 .or. igpnt .eq. 1 ) then
       
        rlat=temp(ic,np+latindex)
        rlon=temp(ic,np+lonindex)
        if ( rlon .gt. bmiss .and. rlat .gt. bmiss ) then

         if ( iposcheck .eq. 1 .and. lons .gt. 0) then
          call grid(rlat,rlon,lat,lon,xgrid)

          if ( lat .ge. j1 .and. lat .le. j2 .and.
     *         (lon .ge. i1 .and. lon .le. i2 .or.
     *          lon .ge. i3 .and. lon .le. i4 ) ) 
     *     iyes=1
         elseif ( igpnt .eq. 1 ) then
          call distance(xlat,rlat,xlon,rlon,xd)
          if ( xd .le. xdist ) iyes=1
         elseif ( lons .lt. 0 ) then
          call grid(rlat,rlon,lat,lon,xgrid)
          if ( m1(lon,lat) .eq. lone ) iyes=1
         endif

        endif
     
       else
        iyes=1
       endif

c      check date

c      if ( iyes .eq. 1 ) write(6,*) 'yes',
c    * idatecheck,temp(ic,np+latindex),
c    * temp(ic,np+lonindex)
       if ( idatecheck .eq. 1 .and. iyes .eq. 1) then

        rjul=temp(ic,np+itimeindex)
        if ( rjul .gt. bmiss ) then

c      if a specific month/day is desired, but not a
c      specific year, reduce the julian day to the 
c      year specified in the header (iyear)

         if ( iyears .eq. 0 .and. iyeare .eq. 0 ) then
          call reducejulian(iyear,rjul,iyearadd)
         endif

        endif

c       write(6,*) 'rjul',rjul,xjdays,xjdaye
        if ( rjul .ge. xjdays .and. rjul .lt. xjdaye ) iyes=2

       else

        iyes=iyes+1

       endif

       if ( iyes .ne. 2 ) then

        if ( n .le. maxlevel ) then
         depth(n)=bmiss
        else
         temp(ic,np)=bmiss
        endif
      
        do 305 ii=1,ireqpars
305      if ( ip2(ii) .ne. 14 )
     *    temp(ic,np+ip2(ii))=bmiss

c    Keep information from within user prescribed positions/dates

       elseif ( imple .eq. 1 ) then

        nlevout=nlevout+1
        
        if ( nlevout .le. maxlevel ) then
         depth(nlevout)=d1
         ic2=nlevout
         np2=0
        else
         call overmax(maxlevel,nparams,nlevout,0,ic2,np2)
         temp(ic2,np2)=d1
        endif
      
        do 300 ii=1,ireqpars
300      if ( ip2(ii) .ne. 14 )
     *    temp(ic2,np2+ip2(ii))=temp(ic,np+ip2(ii))

       else

        nlevout=nlevout+1

       endif

50    continue

      return
      end
