SUBROUTINE stalst(nsta,id,rlat,rlon,nsum, &
     ihsta,jhsta,ivsta,jvsta,kclas)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: imjm=im*jm-jm/2
  INTEGER, PARAMETER :: imt=im*2-1
  INTEGER, PARAMETER :: jmt=jm/2+1
  INTEGER, PARAMETER :: jam=6+2*(jm-10)
  INTEGER, PARAMETER :: lp1=lm+1
  REAL   , PARAMETER :: d2r=1.745329252e-2
  REAL   , PARAMETER :: r2d=57.29577951
  INCLUDE "loops.f90"
  INCLUDE "mapot.f90"
  INCLUDE "masks.f90"
  INCLUDE "indx.f90"
  REAL :: rlat(nsta)
  REAL :: rlon(nsta)
  REAL :: tlat(nsta)
  REAL :: tlon(nsta)
  REAL :: rlatd(nsta)
  REAL :: rlond(nsta)
  INTEGER :: ihdum(nsta)
  INTEGER :: jhdum(nsta)
  INTEGER :: ivdum(nsta)
  INTEGER :: jvdum(nsta)
  INTEGER :: ihsta(nsta)
  INTEGER :: jhsta(nsta)
  INTEGER :: ivsta(nsta)
  INTEGER :: jvsta(nsta)
  INTEGER :: kclas(nsta)
  INTEGER :: kclasd(nsta)
  INTEGER :: iddum(nsta)
  INTEGER :: id(nsta)
  CHARACTER(len=1) :: an,aw
  REAL :: lam0
  INTEGER :: lunlst=15
  REAL :: a
  REAL :: alat
  REAL :: b
  REAL :: c
  REAL :: col
  REAL :: d
  REAL :: d1
  REAL :: d2
  REAL :: d3
  REAL :: d4
  REAL :: dlam
  REAL :: dlm1
  REAL :: dlm2
  REAL :: dlm3
  REAL :: dlm4
  REAL :: dmin
  REAL :: dphi
  INTEGER :: i
  INTEGER :: idum
  INTEGER :: idx
  INTEGER :: ihv
  INTEGER :: imin
  INTEGER :: ix
  INTEGER :: j
  INTEGER :: jdum
  INTEGER :: jmin
  INTEGER :: jx
  INTEGER :: kcls
  INTEGER :: kh
  INTEGER :: kmin
  INTEGER :: krows
  INTEGER :: kv
  INTEGER :: kvmin
  INTEGER :: n2
  INTEGER :: ncol
  INTEGER :: ndum
  INTEGER :: nrow
  INTEGER :: nsta
  INTEGER :: nsum
  REAL :: phi0
  REAL :: rlatx
  REAL :: rlonx
  REAL :: row
  REAL :: sbdb
  REAL :: tlat1
  REAL :: tlat2
  REAL :: tlat3
  REAL :: tlat4
  REAL :: tlats
  REAL :: tlatx
  REAL :: tlon1
  REAL :: tlon2
  REAL :: tlon3
  REAL :: tlon4
  REAL :: tlons
  REAL :: tlonx
  REAL :: wbdb
  REAL :: x
  REAL :: y
  REAL :: z


  phi0=tph0d*d2r
  lam0=-1.0*tlm0d*d2r
  dphi=dphd*d2r
  dlam=dlmd*d2r

  !*** sbdb and wbdb are 2 rows inside the domain boundary
  !*** to try to exclude points right on the boundary

  sbdb=sbd+2.0*dphd
  wbdb=wbd+2.0*dlmd

  !***  loop through all the stations to count how many inside domain

  nsum=0
25 READ(lunlst,50,END=100)idx,rlatx,an,rlonx,aw,kcls
  IF (an.EQ.'S') rlatx=-1.0*rlatx                  !chou
  IF (aw.EQ.'E') rlonx=360.0-rlonx                 !chou
  rlatx=rlatx*d2r
  rlonx=rlonx*d2r

  !***  convert geodetic to transformed coordinates of the station

  x=COS(phi0)*COS(rlatx)*COS(rlonx-lam0) &
       +SIN(phi0)*SIN(rlatx)
  y=-COS(rlatx)*SIN(rlonx-lam0)
  z=COS(phi0)*SIN(rlatx) &
       -SIN(phi0)*COS(rlatx)*COS(rlonx-lam0)
  tlatx=r2d*ATAN(z/SQRT(x*x+y*y))
  tlonx=r2d*ATAN(y/x)
  IF(ABS(tlatx).GT.ABS(sbdb).OR.ABS(tlonx).GT.ABS(wbdb))go to 25
  nsum=nsum+1
  rlatd(nsum)=rlatx
  rlond(nsum)=rlonx
  tlat(nsum)=tlatx
  tlon(nsum)=tlonx
  iddum(nsum)=idx
  kclasd(nsum)=kcls
  GOTO 25
100 CONTINUE

  DO n2=1,nsum

     !*** find the closest point regardless of land/sea mask

     row=tlat(n2)/dphd+jmt
     col=tlon(n2)/dlmd+im
     nrow=INT(row)
     ncol=INT(col)
     tlat(n2)=tlat(n2)*d2r
     tlon(n2)=tlon(n2)*d2r

     !***  identify the row and column of the lower left corner of the box
     !***  that surrounds each station inside the domain

     !***  first consider the situation where the station x is at
     !***
     !***              h      v
     !***
     !***
     !***                 x
     !***              v      h
     !***
     IF(MOD(nrow,2).EQ.0.AND.MOD(ncol,2).EQ.1.OR. &
          MOD(nrow,2).EQ.1.AND.MOD(ncol,2).EQ.0)THEN
        tlat1=(nrow-jmt+1)*dphi
        tlat2=tlat1-dphi
        tlon1=(ncol-im)*dlam
        tlon2=tlon1+dlam
        dlm1=tlon(n2)-tlon1
        dlm2=tlon(n2)-tlon2
        d1=ACOS(COS(tlat(n2))*COS(tlat1)*COS(dlm1) &
             +SIN(tlat(n2))*SIN(tlat1))
        d2=ACOS(COS(tlat(n2))*COS(tlat2)*COS(dlm2) &
             +SIN(tlat(n2))*SIN(tlat2))
        IF(d1.GT.d2)THEN
           ncol=ncol+1
        ELSE
           nrow=nrow+1
        ENDIF

        !***  now consider the situation where the point x is at
        !***
        !***              v      h
        !***
        !***
        !***                 x
        !***              h      v
        !***
     ELSE
        tlat1=(nrow-jmt)*dphi
        tlat2=tlat1+dphi
        tlon1=(ncol-im)*dlam
        tlon2=tlon1+dlam
        dlm1=tlon(n2)-tlon1
        dlm2=tlon(n2)-tlon2
        d1=ACOS(COS(tlat(n2))*COS(tlat1)*COS(dlm1) &
             +SIN(tlat(n2))*SIN(tlat1))
        d2=ACOS(COS(tlat(n2))*COS(tlat2)*COS(dlm2) &
             +SIN(tlat(n2))*SIN(tlat2))
        IF(d1.GT.d2)THEN
           nrow=nrow+1
           ncol=ncol+1
        ENDIF
     ENDIF

     !***  now we can find the k value

     krows=((nrow-1)/2)*imt
     IF(MOD(nrow,2).EQ.1)THEN
        kmin=krows+(ncol+1)/2
     ELSE
        kmin=krows+im+ncol/2
     ENDIF
     jx=(kmin-1)/imt+1
     ix=kmin-(jx-1)*imt
     IF(ix.LE.im)THEN
        idum=ix
        jdum=2*jx-1
     ELSE
        idum=ix-im
        jdum=2*jx
     ENDIF

     IF (kclasd(n2).GE.10.AND.kclasd(n2).LT.20) THEN

        !*** we want this to be a land point, if not, find the closest land point

        IF (sm(idum,jdum).GT.0.5) THEN
           dmin=99999.
           DO j=1,jm
              tlat1=(j-jmt)*dphi
              DO i=1,im
                 ihv=i*2-MOD(j,2)
                 tlon1=(ihv-im)*dlam
                 dlm1=tlon(n2)-tlon1
                 IF(sm(i,j).LT.0.5) THEN
                    d1=ACOS(COS(tlat(n2))*COS(tlat1)*COS(dlm1)+ &
                         SIN(tlat(n2))*SIN(tlat1))
                    IF(d1.LT.dmin) THEN
                       dmin=d1
                       ncol=ihv
                       nrow=j
                    ENDIF
                 ENDIF
              ENDDO
           ENDDO
        ENDIF
        krows=((nrow-1)/2)*imt
        IF(MOD(nrow,2).EQ.1)THEN
           kmin=krows+(ncol+1)/2
        ELSE
           kmin=krows+im+ncol/2
        ENDIF
     ENDIF

     IF (kclasd(n2).GE.20.AND.kclasd(n2).LT.30) THEN

        !*** we want this to be a sea point, if not, find the closest sea point

        IF (sm(idum,jdum).LT.0.5) THEN
           dmin=99999.
           DO j=1,jm
              tlat1=(j-jmt)*dphi
              DO i=1,im
                 ihv=i*2-MOD(j,2)
                 tlon1=(ihv-im)*dlam
                 dlm1=tlon(n2)-tlon1
                 IF(sm(i,j).GT.0.5) THEN
                    d1=ACOS(COS(tlat(n2))*COS(tlat1)*COS(dlm1)+ &
                         SIN(tlat(n2))*SIN(tlat1))
                    IF(d1.LT.dmin) THEN
                       dmin=d1
                       ncol=ihv
                       nrow=j
                    ENDIF
                 ENDIF
              ENDDO
           ENDDO
        ENDIF
        krows=((nrow-1)/2)*imt
        IF(MOD(nrow,2).EQ.1)THEN
           kmin=krows+(ncol+1)/2
        ELSE
           kmin=krows+im+ncol/2
        ENDIF
     ENDIF

     !*** find the geodetic lat/lon

     tlats=(nrow-jmt)*dphi
     tlons=(ncol-im)*dlam
     a=SIN(tlats)*COS(phi0)+COS(tlats)*SIN(phi0)*COS(tlons)
     alat=ASIN(a)
     rlatd(n2)=alat
     b=COS(tlats)*COS(tlons)/(COS(alat)*COS(phi0))
     d=b-TAN(alat)*TAN(phi0)
     IF (ABS(d).GT.1.0) d=ABS(d)/d
     c=ACOS(d)
     c=SIGN(c,-1.0*tlons)
     rlond(n2)=lam0+c

     !***  this is the h point, now find the closest neighboring
     !***  v point that is on the same step
     !***
     !***                 2          kmin+im-1
     !***
     !***              1  h  3  kmin-1  h   kmin
     !***                  
     !***                 4           kmin-im
     !***       

     tlat1=(nrow-jmt)*dphi
     tlon1=(ncol-im-1)*dlam
     tlat2=tlat1+dphi
     tlon2=tlon1+dlam
     tlat3=tlat1
     tlon3=tlon2+dlam
     tlat4=tlat1-dphi
     tlon4=tlon2
     dlm1=tlon(n2)-tlon1
     dlm2=tlon(n2)-tlon2
     dlm3=tlon(n2)-tlon3
     dlm4=tlon(n2)-tlon4
     d1=ACOS(COS(tlat(n2))*COS(tlat1)*COS(dlm1) &
          +SIN(tlat(n2))*SIN(tlat1))
     d2=ACOS(COS(tlat(n2))*COS(tlat2)*COS(dlm2) &
          +SIN(tlat(n2))*SIN(tlat2))
     d3=ACOS(COS(tlat(n2))*COS(tlat3)*COS(dlm3) &
          +SIN(tlat(n2))*SIN(tlat3))
     d4=ACOS(COS(tlat(n2))*COS(tlat4)*COS(dlm4) &
          +SIN(tlat(n2))*SIN(tlat4))
     dmin=99999.
     CALL k2ijh(kmin,imin,jmin)
     IF(lmh(imin,jmin).EQ.lmv(imin+ihw(jmin),jmin).AND. &
          d1.LT.dmin)THEN
        dmin=d1
        kvmin=kmin-1
     ENDIF
     IF(lmh(imin,jmin).EQ.lmv(imin,jmin+1).AND.d2.LT.dmin)THEN
        dmin=d2
        kvmin=kmin+im-1
     ENDIF
     IF(lmh(imin,jmin).EQ.lmv(imin+ihe(jmin),jmin).AND. &
          d3.LT.dmin)THEN
        dmin=d3
        kvmin=kmin
     ENDIF
     IF(lmh(imin,jmin).EQ.lmv(imin,jmin-1).AND.d4.LT.dmin)THEN
        dmin=d4
        kvmin=kmin-im
     ENDIF
     IF(dmin.EQ.99999.)THEN
        kvmin=kmin
        WRITE(6,180)iddum(n2),rlatd(n2)*r2d,rlond(n2)*r2d,kmin, &
             lmh(imin,jmin),lmv(imin+ihw(jmin),jmin),lmv(imin,jmin+1), & 
             lmv(imin+ihe(jmin),jmin),lmv(imin,jmin-1)
     ENDIF

     kh=kmin
     kv=kvmin

     jx=(kh-1)/imt+1
     ix=kh-(jx-1)*imt
     IF(ix.LE.im)THEN
        ihdum(n2)=ix
        jhdum(n2)=2*jx-1
     ELSE
        ihdum(n2)=ix-im
        jhdum(n2)=2*jx
     ENDIF

     jx=(kv-1)/imt+1
     ix=kv-(jx-1)*imt
     IF(ix.LE.im-1)THEN
        ivdum(n2)=ix
        jvdum(n2)=2*jx-1
     ELSE
        ivdum(n2)=ix-im+1
        jvdum(n2)=2*jx
     ENDIF

  END DO

  !***  remove duplicate points

  ndum=0
  DO i=1,nsum
     IF(iddum(i).NE.-99999)THEN
        DO j=1,nsum
           IF(ihdum(i).EQ.ihdum(j).AND.jhdum(i).EQ.jhdum(j).AND. &
                ivdum(i).EQ.ivdum(j).AND.jvdum(i).EQ.jvdum(j).AND. &
                i.NE.j)THEN
              WRITE(6,900)iddum(i),iddum(j),iddum(j)
              ihdum(j)=-99999
              jhdum(j)=-99999
              ivdum(j)=-99999
              jvdum(j)=-99999
              iddum(j)=-99999
           ENDIF
        ENDDO
        ndum=ndum+1
        ihsta(ndum)=ihdum(i)
        jhsta(ndum)=jhdum(i)
        ivsta(ndum)=ivdum(i)
        jvsta(ndum)=jvdum(i)
        id(ndum)=iddum(i)
        rlat(ndum)=rlatd(i)
        rlon(ndum)=rlond(i)
        kclas(ndum)=MOD(kclasd(i),10)
     ENDIF
  ENDDO
  nsum=ndum
50 FORMAT(i5,f6.2,a1,f7.2,a1,4x,i3)
180 FORMAT(i5,2f8.2,i8,5i4,' no neighbor on same step')
900 FORMAT(' stalst',i8,' and',i8,' are at the same grid pts', &
       i8,' is removed')
END SUBROUTINE stalst
