!### 2015/05/03 Yuta Ando (Mie University)

subroutine epflux(u,v,air,zonalair,xmax,ymax,zmax,lat,lev,Fy,Fp)
  implicit none
  integer :: x,y,z,xmax,ymax,zmax
  real,parameter :: a=6.37e+6 !mean radius of the earth[m],pai=4.*atan(1.)
  real,parameter:: Cp=1004.,Rd=287.,omg=7.292e-5,ps0=1.e+3

  real,intent(in) :: u(xmax,ymax,zmax)   !zonal anomaly of u-wind velocity[m/s]
  real,intent(in) :: v(xmax,ymax,zmax)   !zonal anomaly of v-wind velocity[m/s]
  real,intent(in) :: air(xmax,ymax,zmax) !zonal anomaly of temperature[K]
  real,intent(in) :: zonalair(ymax,zmax) !zonal average of temperature[K]
  real,intent(in) :: lat(ymax),lev(zmax) !lat[deg],lev[Pa]

  real,intent(out) :: Fy(ymax,zmax) !meridional component of EP flux
  real,intent(out) :: Fp(ymax,zmax) !vertical component of EP flux

  real :: f(ymax),fai(ymax),dp
  real :: sita(xmax,ymax,zmax)
  real :: zonaldsitadp(ymax,zmax),zonaluv(ymax,zmax),zonalvsita(ymax,zmax)
  real :: zonalsita(ymax,zmax)
  integer :: jz1,jz2


  fai(1:ymax)=pai*lat(1:ymax)/180.
  f(1:ymax)=2.*omg*sin(fai(1:ymax))
  do z=1,zmax
    do y=1,ymax
      do x=1,xmax
        sita(x,y,z)=air(x,y,z)*((ps0*100.0)/lev(z))**(Rd/cp)
      enddo
      zonalsita(y,z)=zonalair(y,z)*((ps0*100.0)/lev(z))**(Rd/cp)
    enddo
  enddo

  zonaldsitadp=0.
  zonaluv=0.
  zonalvsita=0.

  do z=1,zmax
    jz1=z+1
    jz2=z-1
    if(z==1)jz2=1
    if(z==zmax)jz1=zmax
    dp=-(log(lev(jz1))-log(lev(jz2)))
    do y=2,ymax-1
      do x=1,xmax
        zonaluv(y,z)=zonaluv(y,z)+u(x,y,z)*v(x,y,z)/real(xmax)
        zonalvsita(y,z)=zonalvsita(y,z)+v(x,y,z)*sita(x,y,z)/real(xmax)
      enddo
      zonaldsitadp(y,z)=(1.0/lev(z))*(zonalsita(y,jz1)-zonalsita(y,jz2))/dp
    enddo
  enddo


  do z=1,zmax
    do y=2,ymax-1
      Fy(y,z)=cos(fai(y))*(-a*cos(fai(y))*zonaluv(y,z))/(a*pai)
      Fp(y,z)=(cos(fai(y))*(f(y)*a*cos(fai(y))*zonalvsita(y,z)))*(1.0/zonaldsitadp(y,z))/(100000.0)
    enddo
  enddo


  Fy(1,:)=0.
  Fy(ymax,:)=0.
  Fp(1,:)=0.
  Fp(ymax,:)=0.

  do z=1,zmax
    do y=1,ymax
      Fy(y,z)=Fy(y,z)*sqrt(100000.0/lev(z))
      Fp(y,z)=Fp(y,z)*sqrt(100000.0/lev(z))
      if (lev(z)<=10000.0) then
        Fy(y,z)=Fy(y,z)*5.0
        Fp(y,z)=Fp(y,z)*5.0
      endif
    enddo
  enddo


end subroutine