!--Supporting functions/subroutine for TC_detect program

!-- Calculate vorticity, Rx and Ry is in Degree
subroutine vorticity(Vor,U,V,Nx,Ny,Rx,Ry)
implicit none
integer :: Nx,Ny,i,j
real, dimension(Nx,Ny) :: Vor, U, V
real, dimension(Nx) :: Rx
real, dimension(Ny) :: Ry
real :: Pi, Re, Dx, Dy, Fac, Deg2Rad

Pi = 3.14159265 
Re = 6378161.  ! Earth radius in m
Deg2Rad = Pi/180.
Fac = Deg2Rad*Re


Do i=2,Nx-1
   Do j=2,Ny-1
     Dx = (Rx(i+1)-Rx(i-1))*Fac*cos(Deg2Rad*Ry(j))
     Dy = (Ry(j+1)-Ry(j-1))*Fac
     Vor(i,j) = (V(i+1,j) -  V(i-1,j))/Dx - (U(i,j+1) -  U(i,j-1))/Dy
   Enddo
Enddo

end subroutine vorticity

!---
!--- This subroutine calculate the azimuthal mean value (Fmean) of a scalar field F
!--- This subroutine use birational interpolation via birational coeficiton Coef 
!---   that is generated from RAT2D subroutine
!--- This also calculate the center value Fcenter and anomaly value Fano=Fcenter-Fmean
!---
subroutine azimean(F,Nx,Ny,Rx,Ry,Xc,Yc,R,Fmean,Fcenter,Fano,iStatus)
implicit none
   !--Input and output
   integer :: Nx,Ny
   real, dimension(Nx,Ny) :: F
   real, dimension(Nx) :: Rx
   real, dimension(Ny) :: Ry
   real :: Xc,Yc,R,Fmean,Fcenter,Fano      ! R is in Degree (or the same unit with Rx,Ry)
   !--Local variables
   real, dimension(Nx,Ny,4,4) :: Coef
   integer :: i,iErr,Npoints,iStatus
   real :: dangle,angle,Pi,X,Y,pp,Ftmp
   pp=0.1              !  birational interpolation parameter (must >=-1)
   Npoints=8           !  here we use only 8 points around the center 
   Pi = 3.14159265 
   dangle=2*Pi/Npoints
   !write(*,*) "In azimean subroutine..."
   ! Calculate birational int. coef.
   call RAT2D(F,rx,ry,nx,ny,nx,ny,PP,Coef,iErr)
   IF (iErr.EQ.1) THEN
         PRINT*,' azimean: problem with birational interpolation!'
         iStatus=0
         RETURN
   ENDIF

   Fmean=0.
   do i=1,Npoints
      angle=(i-1)*dangle
      X = Xc + R*cos(angle)
      Y = Yc + R*sin(angle)
      !write(*,*)"Azimean:",X,Y,"-->",Ftmp
      !write(*,*) X,Rx(1),Rx(Nx),Y,Ry(1),Ry(Ny) 
      
      if (X<Rx(1) .or. X>Rx(Nx) .or. Y<Ry(1) .or. Y>Ry(Ny)) then
         write(*,*)'  azimean: Out of bound!'
         iStatus=0
         return
      endif
      call RATVA(Coef,Rx,Ry,Nx,Ny,Nx,Ny,PP,X,Y,Ftmp)
      Fmean=Fmean+Ftmp
   enddo
   Fmean=Fmean/Npoints 
   call RATVA(Coef,Rx,Ry,Nx,Ny,Nx,Ny,PP,Xc,Yc,Fcenter)
   Fano = Fcenter - Fmean
   iStatus=1
end subroutine

!---------
!-- Calculate Outer core strength OCS 
!-- OCS = avarage of tangential wind from radius of 1degree to 2.5degree
subroutine calc_OCS_old(U,V,Nx,Ny,Rx,Ry,Xc,Yc,OCS,iStatus)
   integer :: Nx,Ny,i,j,ir,Npoints,iStatus
   real, dimension(Nx,Ny) :: U,V,Vt
   real, dimension(Nx) :: Rx   ! Note: Rx, Ry,Xc,Yc is in Degree 
   real, dimension(Nx) :: Ry
   real :: Xc,Yc,OCS, dx, dy
   real :: Pi, Re, Reg2Rad,Fac, Cosa,Sina,Dis, pp,Ftmp
   real, dimension(Nx,Ny,4,4) :: Coef
   Real, dimension(4) :: R = (/1.,1.5,2.,2.5/)
   pp=0.1
   Pi = 3.14159265 
   Re = 6378161.  ! Earth radius in m
   Deg2Rad = Pi/180.
   Fac = Deg2Rad*Re
   Npoints=8           !  here we use only 8 points around the center 
   dangle=2*Pi/Npoints

! calculate tangential wind
   do i=1,Nx 
      do j=1,Ny
         dy=Ry(j)-Yc
         dx=(Rx(i)-Xc)*cos(Deg2Rad*Ry(j))
         dis=sqrt(dx*dx+dy*dy)
         if (dis.gt.0) then
            cosa=dx/dis
            sina=dy/dis
            Vt(i,j)=-U(i,j)*sina + V(i,j)*cosa
         else
            Vt(i,j)=0 
         endif
      enddo
   enddo
   ! Calculate birational int. coef.
   call RAT2D(Vt,rx,ry,nx,ny,nx,ny,PP,Coef,iErr)
   
   
   OCS=0.
   do ir=1,4
      do i=1,Npoints
      angle=(i-1)*dangle
      X = Xc + R(ir)*cos(angle)
      Y = Yc + R(ir)*sin(angle)

      if (X<Rx(1) .or. X>Rx(Nx) .or. Y<Ry(1) .or. Y>Ry(Ny)) then
         write(*,*)'  azimean: Out of bound!'
         iStatus=0
         return
      endif

      call RATVA(Coef,Rx,Ry,Nx,Ny,Nx,Ny,PP,X,Y,Ftmp)
      OCS=OCS+Ftmp
      enddo
   enddo
   OCS=OCS/Npoints/4
   iStatus=1

end subroutine calc_OCS_old


!---------
!-- Calculate Outer core strength OCS 
!-- OCS = avarage of tangential wind from radius of 1degree to 2.5degree
!-- New subroutine: simpler and faster
subroutine calc_OCS(U,V,Nx,Ny,Rx,Ry,Xc,Yc,OCS)
   integer :: Nx,Ny,i,j,Npoints
   real, dimension(Nx,Ny) :: U,V
   real, dimension(Nx) :: Rx   ! Note: Rx, Ry,Xc,Yc is in Degree 
   real, dimension(Nx) :: Ry
   real :: Xc,Yc,OCS, dx, dy
   real :: Pi, Re, Reg2Rad,Fac, Cosa,Sina,Dis
   Pi = 3.14159265 
   Re = 6378161.  ! Earth radius in m
   Deg2Rad = Pi/180.
  
! calculate tangential wind
   OCS = 0.
   Npoints = 0
   do i=1,Nx 
      do j=1,Ny
         dy=Ry(j)-Yc
         dx=(Rx(i)-Xc)*cos(Deg2Rad*Ry(j))
         dis=sqrt(dx*dx+dy*dy)
         if ((dis.ge.1).and.(dis.le.2.5))then
            cosa=dx/dis
            sina=dy/dis
            OCS = OCS - U(i,j)*sina + V(i,j)*cosa
            Npoints = Npoints+1
         endif 
      enddo
   enddo
   OCS = OCS/Npoints

end subroutine calc_OCS

!-- Calculate distance (km) beetween two point (degree)
subroutine geodistance(lat1,lon1,lat2,lon2,dis)
real :: lat1,lon1,lat2,lon2,dis
real :: Pi,R_e,Fac,lat11,lon11,lat22,lon22,cos_distance,rad_distance
   Pi = 3.14159265 
   R_e=6378.161  ! Earth radius in Km
   FAC=Pi/180.
   lat11=FAC*lat1
   lon11=FAC*lon1
   lat22=FAC*lat2
   lon22=FAC*lon2
   cos_distance= sin(lat11)*sin(lat22) + cos(lat11)*cos(lat22)*cos(lon22-lon11)
   !print "Cos Distance:",cos_distance
   if ((lat1.eq.lat2).and.(lon1.eq.lon2)) then
      rad_distance=0.
   else 
      rad_distance=acos( cos_distance )
   endif
   dis=rad_distance*R_e
   return 
end subroutine geodistance


!
! 4 points moothing
!
subroutine smooth_4p(F,Nx,Ny)
integer :: Nx,Ny, i,j
real, dimension(Nx,Ny) :: F, Ftemp
   Ftemp=F
   do i=2,Nx-1
     do j=2,Ny-1
       Ftemp(i,j) = 0.25* ( F(i-1,j-1) + F(i-1,j+1) + F(i+1,j-1) + F(i+1,j+1)  )
     enddo
   enddo
   F=Ftemp
end subroutine smooth_4p