PROGRAM TRAJECTORY
IMPLICIT NONE
!
!
! indices and counters
integer :: grid !model grid index
integer :: traj !trajectory index
integer :: time !timestep in trajectory integration
!
!
! Constants
real, parameter :: prad=3389.5e03 !radius in meters
real, parameter :: sechr=3698.967 !seconds in Mars hour
real, parameter :: deg2rad=atan(1.0)/45.
real, parameter :: rad2deg=45./atan(1.0)
real  :: wall !actual time
!
!
! Initial Release location and time
real*8, parameter :: lat0= 4.0
real*8, parameter :: lon0= 137.5
real, parameter :: agl0=2.0 !start height above ground level
real, parameter :: t0=21.3*sechr !time of release in seconds from OOOO UTC 
real, parameter :: rel_rad=30.0  !radius for random release 
!
!
! Integration Parameters
!real, parameter :: t_final=47.3*sechr !end time of simulation from 0000 UTC
real, parameter :: t_final=46.*sechr !end time of simulation from 0000 UTC
real, parameter :: tstep=sechr/180. !integration time step in seconds
real*8 :: dx,dlat,dy,dlon,dz
integer, parameter :: tsize=(int(t_final-t0)/tstep)+1
!
!
! Trajectory variable set-up
integer, parameter :: n_traj=1 !number of trajectories to compute
real*8, dimension (n_traj,tsize) :: xlat,xlon
real , dimension (n_traj,tsize) :: agl,u_traj,v_traj,w_traj,xtopo
real :: x,y,z,tindex !fractional x,y,z,t location in model grid space
real :: alt !altitude above aeroid
!
!
! MRAMS model input set-up
integer, parameter :: ngrids=5 !number of availabe MRAMS grids
character(len=80) :: sfcfile="Gale_Ls195_7grid-sfc-S-g0"
character(len=80) :: atmfile="Gale_Ls195_7grid-atm-S-g0"
character(len=80) :: gfile
integer, dimension(ngrids), parameter :: nxp=(/60,59,98,89,98/)
integer, dimension(ngrids), parameter :: nyp=(/45,62,128,110,95/)
integer, parameter :: nxpmax=98!max x dim
integer, parameter :: nypmax=128 !max y dim 
integer, parameter :: nzp=50
real, parameter :: grid_t0=18.0*sechr !starting UTC of output time 1 in s
real, parameter :: grid_dt=sechr/3. !time between output times in s
real, dimension(ngrids,nxpmax,nypmax) :: topo,glat,glon
real, dimension(ngrids,nzp,nxpmax,nypmax) :: umod,vmod,wmod,zmod !temp vars
!
!
! Initialization
!    Read in glat,glon,topo for all grids
DO grid=1,ngrids
        call filename(grid,sfcfile,gfile)
        call grads_read_sfc(gfile,glat(grid,1:nxp(grid),1:nyp(grid)),  &
                            glon(grid,1:nxp(grid),1:nyp(grid)),  &
                            topo(grid,1:nxp(grid),1:nyp(grid)),  &
                            nxp(grid),nyp(grid))
ENDDO
!
! Start Trajectory loop
!print *,'starting at t0 ',t0
!print *,'ending at tfinal ',t_final
!print *,'using time step ',tstep
!print *,'there will be ',tsize,' time steps'
DO traj=1,n_traj
!        print *,'working trajectory ',traj
        time=1 !time step 1 
        xlat(traj,time)=lat0
        xlon(traj,time)=lon0
        agl(traj,time)=agl0 
        wall=t0+(time-1)*tstep
!        print *,'starting at agl0',agl0
        write(*,'(4(A7,1x),3(A6,1X),4(A9,1x))'), &
         & "time","lat","lon","topo","agl","u","v","w","dx","dy","dz"
                x=-1.  ! starting x not yet known from lat/lon
                y=-1.  ! starting y not yet known from lat/lon
                !locate the lat/lon on highest res grid
                grid=ngrids
                DO WHILE (x < 0. .and. y < 0. .and. grid>0)
                        call searchll(nxp(grid),nyp(grid), &
                                      xlat(traj,time),xlon(traj,time), &
                                      glat(grid,1:nxp(grid),1:nyp(grid)), &
                                      glon(grid,1:nxp(grid),1:nyp(grid)), &
                                      x,y)
                         !now have x,y corresponding to xlat,xlon
                         grid=grid-1
                END DO
                grid=grid+1 !reset to current active grid
                call gettopo(nxp(grid),nyp(grid),&
                     topo(grid,1:nxp(grid),1:nyp(grid)),xtopo(traj,time),x,y)
        !
        DO WHILE (time<tsize)
                wall=t0+(time-1)*tstep
                !
                !
                ! find the time output fraction
                call findtime(t0,time,tstep,grid_t0,grid_dt,tindex)
                !
                !  Get u, v, w at point x, y, z
                !
                  call filename(grid,atmfile,gfile)
                  !this call really slows things down because have to open and
                  !read a file every timestep.
                  !maybe find a way to speed this up later.
                  !check and see if we need to read in new time and grid data or
                  !whether the stuff from the last time step is good.
                  call getzuvw(gfile,nzp,nxp(grid),nyp(grid),x,y,tindex, &
                              u_traj(traj,time),v_traj(traj,time), &
                              w_traj(traj,time),agl(traj,time),   &
                              xtopo(traj,time),topo(grid,1:nxp(grid),1:nyp(grid)))
                !ADVECT       
                dx=u_traj(traj,time)*tstep
                dlon=rad2deg*dx/(prad*cos(xlat(traj,time)*deg2rad))
                xlon(traj,time+1)=xlon(traj,time)+dlon
                dy=v_traj(traj,time)*tstep
                dlat=rad2deg*dy/prad
                xlat(traj,time+1)=xlat(traj,time)+dlat
                dz=w_traj(traj,time)*tstep
                !now that we have updated lat,lon we should find the new x and y
                !and get the topo at that location, which is topo at time+1
                x=-1.  ! starting x not yet known from lat/lon
                y=-1.  ! starting y not yet known from lat/lon
                !locate the lat/lon on highest res grid
                grid=ngrids
                DO WHILE (x < 0. .and. y < 0. .and. grid>0)
                        call searchll(nxp(grid),nyp(grid), &
                                      xlat(traj,time+1),xlon(traj,time+1), &
                                      glat(grid,1:nxp(grid),1:nyp(grid)), &
                                      glon(grid,1:nxp(grid),1:nyp(grid)), &
                                      x,y)
                         !now have x,y corresponding to xlat,xlon
                         grid=grid-1
                END DO
                grid=grid+1 !reset to current active grid
                call gettopo(nxp(grid),nyp(grid),&
                     topo(grid,1:nxp(grid),1:nyp(grid)),xtopo(traj,time+1),x,y)
                agl(traj,time+1)=agl(traj,time)+xtopo(traj,time)+dz- &
                                 xtopo(traj,time+1) 
                if(agl(traj,time+1) < 0.) then
                   agl(traj,time+1)=0.001 !don't let go below ground
                   dz=agl(traj,time+1)-agl(traj,time)-xtopo(traj,time+1)+xtopo(traj,time)
                endif
                !xtopo(traj,time+1)=xtopo(traj,time)
          write(*,'(f9.2,1x,2(f7.3,1x),e10.4,1x,f9.3,2(f7.3,1x),4(e10.4,1x))'), &
                 wall,xlat(traj,time+1),xlon(traj,time+1),xtopo(traj,time+1), &
                 agl(traj,time+1), &
                 u_traj(traj,time),v_traj(traj,time),w_traj(traj,time), &
                 dx,dy,dz
                time=time+1
        END DO !end time integration   
END DO  !end trajectory loop

END PROGRAM TRAJECTORY
!
!
!
SUBROUTINE gettopo(nxp,nyp,topo,xtopo,x,y)
IMPLICIT NONE
integer :: nxp,nyp,ix,iy
real, dimension(nxp,nyp) :: topo
real :: xtopo
real :: x,y,fracx,fracy,var5,var6
ix=int(x)
iy=int(y)
fracx=x-ix
fracy=y-iy
!interpolate these in x at bounding y values
!var1 u1 zinterp to ix,iy
!var2 u1 zinterp to ix+1,iy
!var3 u1 zinterp to ix,iy+1
!var4 u1 zinterp to ix+1,iy+1
!     topo(ix+1,iy),topo(ix,iy),topo(ix+1,iy+1),topo(ix,iy+1)

var5=fracx*topo(ix+1,iy)+(1.-fracx)*topo(ix,iy) !valid at iy
var6=fracx*topo(ix+1,iy+1)+(1.-fracx)*topo(ix,iy+1) !valid at iy+1
!now interpolate to single y value
xtopo=fracy*var6+(1.-fracy)*var5 
END SUBROUTINE gettopo
!
!*****************

SUBROUTINE getzuvw(gfile,nzp,nxp,nyp,x,y,t, &
                  u_traj,v_traj,w_traj,agl,xtopo,topo) 
!x,y is the fractional x,y grid location we want to interpolate
!agl is height above ground at x,y that we want
!xtopo is topo at x,y
!height to interpolate to is agl+xtopo=alt
!u_traj,v_traj,w_traj will be the velocities interpolated to x,y,alt
!u_traj,v_traj,w_traj will be the velocities interpolated to x,y,alt

IMPLICIT NONE
integer :: nzp,nxp,nyp,t1,nrecsize,icount,k
integer :: ix,iy,iz,it
real :: u_traj,v_traj,w_traj,alt,x,y,t,agl,xtopo
real :: fracx,fracy,fracz,fract
real :: var1,var2,var3,var4,var5
real, dimension(nxp,nyp) :: topo
real, dimension(nzp,nxp,nyp) :: u1,u2
real, dimension(nzp,nxp,nyp) :: v1,v2
real, dimension(nzp,nxp,nyp) :: w1,w2
real, dimension(nzp,nxp,nyp) :: agl1,agl2
character(len=80) :: gfile
nrecsize=nxp*nyp*4
!need to get records bounding t
t1=int(t)
open(unit=21,file=gfile,access='direct',recl=nrecsize)
!skip to the correct record that is t1 output times
icount=4*nzp*t1 + 1 
!read z_lyrmid_agl
do k=1,nzp
  read(21,rec=icount) agl1(k,:,:)
  icount=icount+1
enddo
agl1(1,:,:)=0.
do k=1,nzp
  read(21,rec=icount) u1(k,:,:)
  icount=icount+1
enddo
u1(1,:,:)=u1(2,:,:)
u1=0.
do k=1,nzp
  read(21,rec=icount) v1(k,:,:)
  icount=icount+1
enddo
v1(1,:,:)=v1(2,:,:)
v1=0.
do k=1,nzp
  read(21,rec=icount) w1(k,:,:)
  icount=icount+1
enddo
w1=1.
w1(1,:,:)=0.
do k=1,nzp
  read(21,rec=icount) agl2(k,:,:)
  icount=icount+1
enddo
agl2(1,:,:)=0.
do k=1,nzp
  read(21,rec=icount) u2(k,:,:)
  icount=icount+1
enddo
u2(1,:,:)=u2(2,:,:)
u2=0.
do k=1,nzp
  read(21,rec=icount) v2(k,:,:)
  icount=icount+1
enddo
v2(1,:,:)=v2(2,:,:)
v2=0.
do k=1,nzp
  read(21,rec=icount) w2(k,:,:)
  icount=icount+1
enddo
w2=1.
! set lowest level to w=0 rather than whatever it may be
!w2(1,:,:)=0.
!given the x,y locations and the altitude to interpolate 
!do vertical inteporlation and then horizontal interpolation
alt=xtopo+agl !altitude is the topo value + the agl
call varinterp(u1,u2,agl1,agl2,topo,nzp,nxp,nyp,x,y,alt,t,0,u_traj)
call varinterp(v1,v2,agl1,agl2,topo,nzp,nxp,nyp,x,y,alt,t,0,v_traj)
call varinterp(w1,w2,agl1,agl2,topo,nzp,nxp,nyp,x,y,alt,t,1,w_traj)
close(21)
return
END SUBROUTINE getzuvw

!
!
!
!****************
SUBROUTINE varinterp(u1,u2,agl1,agl2,topo,nzp,nxp,nyp,x,y,alt,t,wtopo,avg)
!alt is aeroid level we want
!agl1 and agl2 are model AGL levels at time1 and time2 (they are the same)
IMPLICIT NONE
INTEGER :: nzp,nxp,nyp,wtopo !wtopo is flag for cases of w and topo
real, dimension(nzp,nxp,nyp) :: u1,u2,agl1,agl2
real, dimension(nxp,nyp) :: topo
real :: x,y,alt,agl,t,fracx,fracy,fracz,fract,var1,var2,var3,var4,avg
real :: var5,var6,var7,var8
INTEGER :: ix,iy,iz,it
ix=int(x)
iy=int(y)
it=int(t)
fracx=x-ix
fracy=y-iy
fract=t-it
! remember that the first t-point is k=1 which should be below ground
! but it's now been hardwired to be z=0 when data is read in. 
!
!first do time 1
!for model columns surrounding point x,y, do a vertical interpolation
!to desired altitude level
!
if(wtopo==0) then
call zinterp(u1(1:nzp,ix,iy),topo(ix,iy),agl1(1:nzp,ix,iy),nzp,alt,var1)
call zinterp(u1(1:nzp,ix+1,iy),topo(ix+1,iy),agl1(1:nzp,ix+1,iy),nzp,alt,var2)
call zinterp(u1(1:nzp,ix,iy+1),topo(ix,iy+1),agl1(1:nzp,ix,iy+1),nzp,alt,var3)
call zinterp(u1(1:nzp,ix+1,iy+1),topo(ix+1,iy+1),agl1(1:nzp,ix+1,iy+1),nzp,alt,var4)
elseif (wtopo==1) then
call zinterpw(u1(1:nzp,ix,iy),topo(ix,iy),agl1(1:nzp,ix,iy),nzp,alt,var1)
call zinterpw(u1(1:nzp,ix+1,iy),topo(ix+1,iy),agl1(1:nzp,ix+1,iy),nzp,alt,var2)
call zinterpw(u1(1:nzp,ix,iy+1),topo(ix,iy+1),agl1(1:nzp,ix,iy+1),nzp,alt,var3)
call zinterpw(u1(1:nzp,ix+1,iy+1),topo(ix+1,iy+1),agl1(1:nzp,ix+1,iy+1),nzp,alt,var4)
else
! do topo interpolation
endif

!var1 u1 zinterp to ix,iy
!var2 u1 zinterp to ix+1,iy
!var3 u1 zinterp to ix,iy+1
!var4 u1 zinterp to ix+1,iy+1
!
!interpolate these in x at bounding y values, but check for data below ground
if(var2 < -9998. .and. var1 < -9998) then
        var5=-9999.
elseif (var2 < -9998.) then
        var5=var1
elseif (var1 < -9998.) then
       var5=var2
else
       var5=fracx*var2+(1.-fracx)*var1 !valid at iy
endif
if(var4 < -9998. .and. var3 < -9998.) then
        var6=-9999.
        elseif (var4 < -9998.) then
                var6=var3
        elseif (var3 < -9998.) then
                var6=var4
        else
                var6=fracx*var4+(1.-fracx)*var3 !valid at iy
endif
!now interpolate to single y value

if(var6 <= -9998. .and. var5 <= -9998.) then
        var7=-9999.
        print *,'no valid value...all below ground'
        stop
        elseif (var6 <= -9998.) then
                var7=var5
        elseif (var5 <= -9998.) then
                var7=var6
        else
                var7=fracy*var6+(1.-fracy)*var5 
endif
!var 7 is now the interpolated value at x,y,z for time it
!
!repeat at the next time
!
!interpolate in z
if(wtopo.eq.0) then
call zinterp(u2(1:nzp,ix,iy),topo(ix,iy),agl2(1:nzp,ix,iy),nzp,alt,var1)
call zinterp(u2(1:nzp,ix+1,iy),topo(ix+1,iy),agl2(1:nzp,ix+1,iy),nzp,alt,var2)
call zinterp(u2(1:nzp,ix,iy+1),topo(ix,iy+1),agl2(1:nzp,ix,iy+1),nzp,alt,var3)
call zinterp(u2(1:nzp,ix+1,iy+1),topo(ix+1,iy+1),agl2(1:nzp,ix+1,iy+1),nzp,alt,var4)
elseif(wtopo.eq.1) then
call zinterpw(u2(1:nzp,ix,iy),topo(ix,iy),agl2(1:nzp,ix,iy),nzp,alt,var1)
call zinterpw(u2(1:nzp,ix+1,iy),topo(ix+1,iy),agl2(1:nzp,ix+1,iy),nzp,alt,var2)
call zinterpw(u2(1:nzp,ix,iy+1),topo(ix,iy+1),agl2(1:nzp,ix,iy+1),nzp,alt,var3)
call zinterpw(u2(1:nzp,ix+1,iy+1),topo(ix+1,iy+1),agl2(1:nzp,ix+1,iy+1),nzp,alt,var4)
else
!do topo
endif
!var1 u1 zinterp to ix,iy
!var2 u1 zinterp to ix+1,iy
!var3 u1 zinterp to ix,iy+1
!var4 u1 zinterp to ix+1,iy+1
!
!interpolate these in x at bounding y values, but check for data below ground

if(var2 < -9998. .and. var1 < -9998) then
        var5=-9998.
elseif (var2 < -9998.) then
        var5=var1
elseif (var1 < -9998.) then
       var5=var2
else
       var5=fracx*var2+(1.-fracx)*var1 !valid at iy
endif
if(var4 < -9998. .and. var3 < -9998.) then
        var6=-9998.
elseif (var4 < -9998.) then
       var6=var3
elseif (var3 < -9998.) then
       var6=var4
        else
                var6=fracx*var4+(1.-fracx)*var3 !valid at iy
endif
!now interpolate to single y value

if(var6 <= -9997. .and. var5 <= -9997.) then
        var8=-9997.
        stop
        elseif (var6 <= -9997.) then
                var8=var5
        elseif (var5 <= -9997.) then
                var8=var6
        else
                var8=fracy*var6+(1.-fracy)*var5 
endif
!
!avg is now the interpolated value at x,y,z for time it+1
!should have a valid value at each time
avg=fract*var8+(1-fract)*var7
!
END SUBROUTINE varinterp
!
!
!
!*************
SUBROUTINE zinterp(var,topo,agl,nzp,alt,var1)
IMPLICIT NONE
INTEGER :: nzp,k,ksave
real, dimension(nzp) :: var !variable to be interpolated
real, dimension(nzp) :: agl !model height array 
real :: alt,var1  !aeroid altitude level to be interpolated and the asnwer
real :: malt1,malt2 ! model altitude
real :: frac !linear interpolation weighting factor
real :: topo !model topo 
ksave=-1
k=2
if(alt < topo) then !desired level is below ground
    var1=-99999. !flag as below ground
    return
endif
if(alt.eq.topo) then !don't bother doing
       var1=var(1) 
       return
endif 
do while(ksave<0 .and. k<nzp+1)
        malt1=agl(k-1)+topo
        malt2=agl(k)+topo
        if(alt>=malt1 .and. alt<malt2) then
           ksave=k !level above desired agl
        endif
        k=k+1
enddo
!do k=2,nzp
!        malt1=agl(k-1)+topo
!        malt2=agl(k)+topo
!        if(alt>=malt1 .and. alt<malt2) then
!           ksave=k !level at desired agl
!        endif
!enddo
if(ksave < 0) then 
  print *,'could not find desired level =',alt
  print *,agl(1),agl(2),malt1,malt2,topo
  stop
endif
frac=(alt-malt1)/(malt2-malt1)
var1=frac*var(ksave)+(1.-frac)*var(ksave-1)        
return
END SUBROUTINE zinterp
!
!
!
!
!*************
SUBROUTINE zinterpw(var,topo,agl,nzp,alt,var1)
IMPLICIT NONE
INTEGER :: nzp,k,ksave
real, dimension(nzp) :: var !variable to be interpolated
real, dimension(nzp) :: agl !model height array
real :: alt,var1  !agl level to be interpolated and the asnwer
real :: malt1,malt2
real :: frac !linear interpolation weighting factor
real :: topo !model topo
ksave=-1
k=2
!it's possible that the altitude we want, alt, is below ground of
!the column we are looking at 
if(alt .lt. topo) then
       var1=-9999.
       return
endif 
if(alt.eq.topo) then !don't bother doing
       var1=var(1) 
       return
endif 
  do while(ksave<0 .and. k<nzp+1)
        malt1=agl(k-1)+topo
        malt2=agl(k)+topo
        if(alt>=malt1 .and. alt<malt2) then
           ksave=k !level above desired agl
        endif
        k=k+1
   enddo
if(ksave < 0) then
  print *,'could not find desired level =',alt
endif
if(ksave>2) then !deal with w=0 b.c.
        frac=(alt-malt1)/(malt2-malt1)
        var1=frac*var(ksave)+(1.-frac)*var(ksave-1)        
else
        frac=(alt-topo)/(malt2-topo)
        var1=frac*var(ksave)+(1.-frac)*0. !force w=0
endif
return
END SUBROUTINE zinterpw
!
!**********************
SUBROUTINE findtime(t0,time,tstep,grid_t0,grid_dt,tindex)
IMPLICIT NONE
integer :: time !integration time step
real :: t0, grid_t0, grid_dt, t,delt,tindex,tstep
t=t0+(time-1)*tstep !current time in seconds
delt=t-grid_t0 !time since first model output time
tindex=(delt/grid_dt)+1. !fractional time output
return
END SUBROUTINE findtime




!
!*******************************
SUBROUTINE searchll(nxp,nyp,xlat,xlon,glat,glon,x,y)
IMPLICIT NONE
integer :: nxp,nyp
real, dimension(nxp,nyp) :: glat,glon
real*8 :: xlat,xlon !lat and lon we want
real :: x, y !x,y location in model grid of xlat,xlon
real :: dist1,dist2,weighti,weightj,wtot
integer :: isave,jsave
integer :: i,j
!find the smallest distance between target and model point
dist1=1.e31
dist2=1.e31
isave=-1
jsave=-1
do i=1,nxp
        do j=1,nyp
          dist1=(xlat-glat(i,j))*(xlat-glat(i,j))+  &
               (xlon-glon(i,j))*(xlon-glon(i,j))
          if(dist1 < dist2) then
              isave=i
              jsave=j
              dist2=dist1     
          endif
        enddo
enddo
!
! Make sure it is not on the edge of grid
! point must be surrounded by other valid points
if (isave==1 .or. isave==nxp .or. jsave==1 .or. jsave==nyp) then
      isave=-1
      jsave=-1
      x=-1.
      y=-1.
      !print *,'location not found on current grid'
      return
endif
!
! do a 9-point weighted average
wtot=0.
weighti=0.
weightj=0.
do i=isave-1,isave+1
        do j=jsave-1,jsave+1
          dist1=(xlat-glat(i,j))*(xlat-glat(i,j))+  &
               (xlon-glon(i,j))*(xlon-glon(i,j))
          dist1=max(dist1,1.e-6)
          weighti=weighti+float(i)/dist1
          weightj=weightj+float(j)/dist1
          wtot=wtot+1./dist1
        enddo
enddo
x=weighti/wtot
y=weightj/wtot
return
!
!
!
!********
END SUBROUTINE searchll
!
!
!
!*******************************
SUBROUTINE searchll2(nxp,nyp,xlat,xlon,glat,glon,x,y)
! an trial alternate search schemd
IMPLICIT NONE
integer :: nxp,nyp
real, dimension(nxp,nyp) :: glat,glon
real*8 :: xlat,xlon !lat and lon we want
real :: x, y !x,y location in model grid of xlat,xlon
real :: dist1,dist2,dist3,weighti,weightj,wtot
integer :: isave,jsave
integer :: i,j,istart,iseek,jstart,jseek
isave=-1
jsave=-1
istart=1
jstart=1

!find the first i,j point to exceed latitude
istart=nxp/2
jstart=nyp/2
iseek=0
jseek=0
dist3=-2.
dist2=-2.
dist1=1.
!let's find out which direction to go
dist1=(xlat-glat(istart+iseek,jstart))*(xlat-glat(istart,jstart))+  &
      (xlon-glon(istart+iseek,jstart))*(xlon-glon(istart,jstart))
dist2=(xlat-glat(istart+iseek+1,jstart))*(xlat-glat(istart+iseek+1,jstart))+  &
      (xlon-glon(istart+iseek+1,jstart))*(xlon-glon(istart+iseek+1,jstart))
dist3=(xlat-glat(istart+iseek-1,jstart))*(xlat-glat(istart+iseek-1,jstart))+  &
      (xlon-glon(istart+iseek-1,jstart))*(xlon-glon(istart+iseek-1,jstart))
print *,'initial: nxp, istart, dist3,dist1,dist2,iseek',nxp,istart,dist3,dist1,dist2,iseek
do while(dist2 .le. dist1)
    print *,'dist2 <= dist1, iseek',dist2,dist1,iseek
    iseek=iseek+1
    dist1=dist2
    dist2=(xlat-glat(istart+iseek+1,jstart))*(xlat-glat(istart+iseek+1,jstart))+  &
      (xlon-glon(istart+iseek+1,jstart))*(xlon-glon(istart+iseek+1,jstart))
enddo
do while(dist3 .lt. dist1)
    print *,'dist3 < dist1, iseek',dist3,dist1,iseek
    iseek=iseek-1
    dist1=dist3
    dist3=(xlat-glat(istart+iseek-1,jstart))*(xlat-glat(istart+iseek-1,jstart))+  &
      (xlon-glon(istart+iseek-1,jstart))*(xlon-glon(istart+iseek-1,jstart))
enddo
!now we know the i point = istart+iseek
!find the jpoint
dist1=(xlat-glat(istart+iseek,jstart))*(xlat-glat(istart+iseek,jstart))+  &
      (xlon-glon(istart+iseek,jstart))*(xlon-glon(istart+iseek,jstart))
dist2=(xlat-glat(istart+iseek,jstart+jseek+1))*(xlat-glat(istart+iseek+1,jstart+jseek+1))+  &
      (xlon-glon(istart+iseek,jstart+jseek+1))*(xlon-glon(istart+iseek+1,jstart+jseek+1))
dist3=(xlat-glat(istart+iseek,jstart+jseek-1))*(xlat-glat(istart+iseek,jstart+jseek-1))+  &
      (xlon-glon(istart+iseek,jstart+jseek-1))*(xlon-glon(istart+iseek,jstart+jseek-1))
print *,'initial: nyp, jstart, dist3,dist1,dist2,jseek',nyp,jstart,dist3,dist1,dist2,iseek
do while(dist2 .le. dist1)
    print *,'dist2 <= dist1, jseek',dist2,dist1,jseek
    jseek=jseek+1
    dist1=dist2
    dist2=(xlat-glat(istart+iseek,jstart+jseek+1))*(xlat-glat(istart+iseek+1,jstart+jseek+1))+  &
      (xlon-glon(istart+iseek,jstart+jseek+1))*(xlon-glon(istart+iseek+1,jstart+jseek+1))
enddo
do while(dist3 .lt. dist1)
    print *,'dist3 < dist1, jseek',dist3,dist1,jseek
    jseek=jseek-1
    dist1=dist3
    dist3=(xlat-glat(istart+iseek,jstart+jseek-1))*(xlat-glat(istart+iseek,jstart+jseek-1))+  &
      (xlon-glon(istart+iseek,jstart+jseek-1))*(xlon-glon(istart+iseek,jstart+jseek-1))
enddo
print *,'istart,iseek,jstart,jseek',istart,iseek,jstart,jseek
print *,'isave jsave ',istart+iseek,jstart+jseek

 
return
!
!
!
!********
END SUBROUTINE searchll2
!
!
!
SUBROUTINE filename(grid,sfcfile,gfile)
integer :: grid
character(len=80) :: sfcfile, gfile
write(gfile,"(A,i1,a)")trim(sfcfile),grid,".dat"
return
END SUBROUTINE filename

!
!
!******************************
SUBROUTINE grads_read_sfc(sfcfile,glat,glon,topo,nxp,nyp)
character(len=*) :: sfcfile
integer :: nxp,nyp
real, dimension(nxp,nyp) :: glat,glon,topo
integer :: nrecsize,icount
nrecsize=nxp*nyp*4
open(unit=21,file=sfcfile,access='direct',recl=nrecsize)
icount=0
         ICOUNT=ICOUNT+1  
         read(21,rec=ICOUNT) topo
         ICOUNT=ICOUNT+1  
         read(21,rec=ICOUNT) glat 
         ICOUNT=ICOUNT+1  
         read(21,rec=ICOUNT) glon 
close(21)
return
END SUBROUTINE grads_read_sfc