c-----*----------------------------------------------------6---------7--
c     Computing the normalized seasonality suv of vector field V=ui+vj,
c     where u and v denote zonal and meridional winds, respectively.
c        suv=DT(V)=||(V1-V2)||/||(V3+V4)/2||,
c     where ||*|| is a norm on the domain S of integration.
c        S=[2*grdx*ma]X[2*grdy*mb], ma and mb are input integers.
c        Generally, ma=mb=1. For the static normalized seasonality (SNS),
c          V1=V3 is the climatological January wind vector, and 
c          V2=V4 is the climatological July wind vector.
c        However, for the dynamical normalized seasonality (DNS), 
c          they are different, i.e., 
c          V1=V3 is the climatological January wind,
c          V4 is the climatological July wind vector, and 
c          V2 is summer (e.g., June, July, August, or their average, etc.) wind vector every year.
c     input: mx,my,v1,v2,v3,v4,fai0,grdx
c       mx: number of gird-point of the fileds in x direction
c       my: number of grid-point of the fileds in y direction
c       v1(,,1): eastward velocity component of V1;
c       v1(,,2): northward velocity component of V1;
c       v2, v3 and v4: V2, V3 and V4 (like V1 mentioned above);
c       fai0: the initial latitude(e.g., 90.,87.5,-90.,etc.)
c             fai0>0 denotes the data starting from the North Pole.
c             fai0<0 denotes the data starting from the South Pole.
c       grdy: grid resolution in y direction (units in degree), e.g. 
c             grdy=2.5 for NCEP/NCAR reanalysis data.
c     output: suv(mx,my)
c       suv: normalized seasonality filed
c     work array: vw1,vw2,vw3,vw4,uv1,uv2
c     Refs.
c     Li, J.P., and Q.C. Zeng, 2000: Significance of the normalized seasonality of wind field and
c        its rationality for characterizing the monsoon. Science in China (D), 43(6), 647-653.
c     Li, J.P., and Q.C. Zeng, 2002: A unified monsoon index, Geophysical Research Letters, 29(8), 1151-1154.
c     Li, J.P., and Q.C. Zeng, 2003: A new monsoon index and the geographical distribution of 
c        the global monsoon. Adv. Atmos. Sci., 20, 299-302.
      subroutine seasonality(mx,my,v1,v2,v3,v4,fai0,grdy,suv)
      parameter(ma=1,mb=1)
      dimension v1(mx,my,2),v2(mx,my,2),v3(mx,my,2),v4(mx,my,2)
      dimension suv(mx,my)
      dimension vw1((-ma+1):(mx+ma),(-mb+1):(my+mb),2) !work array
      dimension vw2((-ma+1):(mx+ma),(-mb+1):(my+mb),2) !work array
      dimension vw3((-ma+1):(mx+ma),(-mb+1):(my+mb),2) !work array
      dimension vw4((-ma+1):(mx+ma),(-mb+1):(my+mb),2) !work array
      dimension uv1((-ma+1):(mx+ma),(-mb+1):(my+mb)) !work array
      dimension uv2((-ma+1):(mx+ma),(-mb+1):(my+mb)) !work array
      mc=-ma+1
      md=-mb+1
      do k=1,2
        call processing(mx,my,ma,mb,fai0,grdy,v1(1,1,k),vw1(mc,md,k))
        call processing(mx,my,ma,mb,fai0,grdy,v2(1,1,k),vw2(mc,md,k))
        call processing(mx,my,ma,mb,fai0,grdy,v3(1,1,k),vw3(mc,md,k))
        call processing(mx,my,ma,mb,fai0,grdy,v4(1,1,k),vw4(mc,md,k))
      enddo
      do j=-mb+1,my+mb
        do i=-ma+1,mx+ma
          a1=(vw1(i,j,1)-vw2(i,j,1))**2
          b1=(vw1(i,j,2)-vw2(i,j,2))**2
          uv1(i,j)=a1+b1
          a2=(vw3(i,j,1)+vw4(i,j,1))**2
          b2=(vw3(i,j,2)+vw4(i,j,2))**2
          uv2(i,j)=a2+b2
	  enddo
	enddo
      do j=1,my
        do i=1,mx
          su1=4.*uv1(i,j)+uv1(i-1,j)+uv1(i+1,j)+uv1(i,j-1)+uv1(i,j+1)
          su2=4.*uv2(i,j)+uv2(i-1,j)+uv2(i+1,j)+uv2(i,j-1)+uv2(i,j+1)
          suv(i,j)=2.*sqrt(su1)/sqrt(su2)
	  enddo
	enddo
      return
      end
c-----*----------------------------------------------------6---------7--
c     Preprocessing for the original field.
c     input: mx,my,ma,mb,fai0,grdx,u
c       fai0: the initial latitude(e.g., 90.,87.5,-90.,etc.)
c             fai0>0 denotes the data starting from the North Pole.
c             fai0<0 denotes the data starting from the South Pole.
c       grdy: grid resolution in y direction (units in degree).
c     output: wu
      subroutine processing(mx,my,ma,mb,fai0,grdy,u,wu)
      dimension u(mx,my),wu((-ma+1):(mx+ma),(-mb+1):(my+mb))
      dimension fai(my) !work array
      pi=3.1415926
      fai(1)=-pi/2. 
      if(fai0.ge.0.) then
        do j=1,my
          fai(j)=(fai0-grdy*(j-1))*pi/180.
        enddo
      else
        do j=1,my
          fai(j)=(fai0+grdy*(j-1))*pi/180.
        enddo
      endif
      do j=1,my
      do i=1,mx
        wu(i,j)=u(i,j)*sqrt(cos(fai(j))) 
      enddo
      enddo
      if(ma.eq.0.and.mb.eq.0) return
      if(ma.gt.0)then
        do j=1,my
          do i=-ma+1,0
            wu(i,j)=wu(mx+i,j)
          enddo
          do i=1,ma
            wu(mx+i,j)=wu(i,j)
          enddo
        enddo
      endif
	mx1=mx/2
      if(mb.gt.0)then
        do i=-ma+1,mx1-1
          do j=-mb+1,0
            wu(i,j)=wu(i+mx1,-j+1)
          enddo
          do j=1,mb
            wu(i,my+j)=wu(i+mx1,my-j+1)
          enddo
	  enddo
        do i=mx1,mx+ma
          do j=-mb+1,0
            wu(i,j)=wu(i-mx1,-j+1)
          enddo
          do j=1,mb
            wu(i,my+j)=wu(i-mx1,my-j+1)
          enddo
	  enddo
      endif
      return
      end