!### 2015/06/21 Yuta Ando (Mie University)

subroutine advection(u,v,tmp,xmax,ymax,zmax,lon,lat,lev,uT,vT,undef)
  implicit none
  real,parameter :: pai=4.*atan(1.)
  integer :: x,y,z,xmax,ymax,zmax

  real,intent(in) :: u(xmax,ymax,zmax)   !u-wind velcity
  real,intent(in) :: v(xmax,ymax,zmax)   !v-wind velcity
  real,intent(in) :: tmp(xmax,ymax,zmax) !temperature
  real,intent(in) :: lon(xmax),lat(ymax),lev(zmax) !lon[deg] lat[deg] lev[hPa]

  real,intent(out) :: uT(xmax,ymax,zmax),vT(xmax,ymax,zmax)

  real :: dy,dx
  real :: dTy(xmax,ymax,zmax),dTx(xmax,ymax,zmax)
  integer :: jx1,jx2,jz1,jz2

  real :: undef

  do z=1,zmax
    do y=2,ymax-1
      dy=(lat(y+1)-lat(y-1))*((2.0*pai*6378150)/360.0)
      do x=1,xmax
        jx1=x+1
        jx2=x-1
        if(x==1)jx2=xmax
        if(x==xmax)jx1=1
        dx=(lon(jx1)-lon(jx2))*cos(lat(y)*pai/180.0)*((2.0*pai*6378150)/360.0)
        dTx(x,y,z)=(tmp(jx1,y,z)-tmp(jx2,y,z))/dx
        if (tmp(jx1,y,z)==undef .or. tmp(jx2,y,z)==undef) then
          dTx(x,y,z)=undef
        endif
        dTy(x,y,z)=(tmp(x,y+1,z)-tmp(x,y-1,z))/dy
        if (tmp(x,y+1,z)==undef .or. tmp(x,y-1,z)==undef) then
          dTy(x,y,z)=undef
        endif
      enddo !x=1,xmax
    enddo !y=2,ymax-1
  enddo !z=1,zmax

  do z=1,zmax
    do y=2,ymax-1
      do x=1,xmax
        uT(x,y,z)=-u(x,y,z)*dTx(x,y,z)
        vT(x,y,z)=-v(x,y,z)*dTy(x,y,z)
        if (u(x,y,z)==undef .or. dTx(x,y,z)==undef) then
          uT(x,y,z)=undef
        endif
        if (v(x,y,z)==undef .or. dTy(x,y,z)==undef) then
          vT(x,y,z)=undef
        endif
      enddo
    enddo
  enddo

  do x=1,xmax
    uT(x,1,:)=undef
    vT(x,1,:)=undef
    uT(x,ymax,:)=undef
    vT(x,ymax,:)=undef
  enddo

  do z=1,zmax
    do y=1,ymax
      do x=1,xmax
        uT(x,y,z)=uT(x,y,z)
        vT(x,y,z)=vT(x,y,z)
      enddo
    enddo
  enddo

end subroutine