! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10)
      Use NetCDF
      Implicit None
      Real, Parameter :: Miss=-99.0, Lat_min=-11.0  ! Cat phan Autralia
      Integer, parameter :: Yr1=1961, Yr2=2015, Nmon=12
      Integer, Parameter :: Yr_sta=1948, Nidx=24
      Real, Parameter :: Slon1=60, Slon2=180, Slat1=-20, Slat2=50
      ! "drawrec 89 147 -11 27 "
      Real, Parameter :: lon1_sea=89.0, lon2_sea=147.0, lat1_sea=-11.0, lat2_sea=27.0
      Real, Parameter :: Miss_quv=-9.96921e+36

      !---------------------------------------------------------------------------------

      Integer :: Nyear ! Thoi ky

      Integer :: i,j,L
      ! Mo ta thong tin SL dau vao:
      !---------------------------------------------------------------------------------
      Real, allocatable :: lon(:), lat(:), lev(:), xtime(:)
      Integer :: Nlon, Nlat, Nlev, Ntime, ilon, ilat
      Logical :: Status

      Real, allocatable :: qu_all(:,:,:), qv_all(:,:,:)   ! Nlon,Nlat,Ntime

      Real, allocatable :: q_flux(:,:)                             ! Ntime
      !---------------------------------------------------------------------------------
      Character (Len=255) :: Fname
      !---------------------------------------------------------------------------------

      Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/quv.mon.mean.nc"
      Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,Ntime,qu_all,"qu")
      Deallocate (lon,lat,xtime)
      Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,Ntime,qv_all,"qv")

      Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/" &
              //"Moisture_flux.nc"
      Call Moisture_Flux (Fname) 

      CONTAINS
      !================================================================================= 
      Subroutine Moisture_Flux (Fname)
        Character (Len=*) :: Fname
        Integer :: it, i1,i2,j1,j2,ilon,ilat, i,j, imon, imon1, Nc, Nc1
        Real :: Qx(Nlon,Nlat), Qy(Nlon,Nlat), tmp(4)
        Real :: lon1,lat1,lon2,lat2
        Real :: Flux, x1,x2,y1,y2,dx,dy,dqx,dqy
        Integer :: Nx1, Nx2, Ny1, Ny2, Nx, Ny, Move_Ave, M2, Nt
        Real, Allocatable :: xlon(:),xlat(:), X1d(:), Y1d(:)
        Real :: Xtb(Nmon), Ytb(Nmon), tmp0, tmp1

        Move_Ave=13   ! Moving average
        M2 = Move_Ave/2; Nt = Ntime-Move_Ave+1
        lon1=lon1_sea; lat1=lat1_sea; lon2=lon2_sea; lat2=lat2_sea
        Do ilon=1,Nlon-1
           if (Lon(ilon)<=lon1.AND.Lon(ilon+1)>=lon1) then
              i1=ilon; exit
           endif
        Enddo
        Do ilon=1,Nlon-1
           if (Lon(ilon)<=lon2.AND.Lon(ilon+1)>=lon2) then
              i2=ilon+1; exit
           endif
        Enddo
        Do ilat=1,Nlat-1
           if (Lat(ilat)<=lat1.AND.Lat(ilat+1)>=lat1) then
              j1=ilat; exit
           endif
        Enddo
        Do ilat=1,Nlat-1
           if (Lat(ilat)<=lat2.AND.Lat(ilat+1)>=lat2) then
              j2=ilat+1; exit
           endif
        Enddo
        !-------------------------------------------
        !print*,"i1,i2,j1,j2=",i1,i2,j1,j2
        Nx1=i1; Nx2=i2; Ny1=j1; Ny2=j2
        Nx=Nx2-Nx1+1; Ny=Ny2-Ny1+1
        Allocate (xlon(Nx),xlat(Ny))
        xlon(1:Nx)=Lon(Nx1:Nx2); xlat(1:Ny)=Lat(Ny1:Ny2)

        Allocate (q_flux(Ntime,4), X1d(Ntime), Y1d(Ntime))
        !q_flux(Ntime,1) : monthly
        !q_flux(Ntime,2) : monthly moving average
        !q_flux(Ntime,3) : monthly Anomaly
        q_flux = Miss; tmp(1)=1; tmp(2)=2; tmp(3)=3; tmp(4)=4

        Do it=1,Ntime
           Qx(:,:) = qu_all(:,:,it); Qy(:,:) = qv_all(:,:,it)
           Flux=0
           Do i=1,Nx-1
              lon1=xlon(i); lon2=xlon(i+1)
              Do j=1,Ny-1
                 lat1=xlat(j); lat2=xlat(j+1)
                 dx=Distance(lon1,lat1,lon2,lat1)
                 dy=Distance(lon1,lat1,lon1,lat2)
                 dqx=Qx(i+1,j)-Qx(i,j)
                 dqy=Qy(i,j+1)-Qy(i,j)
                 Flux=Flux+(dqy/dx-dqx/dy)*dx*dy
              Enddo
             ! pause
           Enddo
           q_flux(it,1) = Flux 
        Enddo
        X1d(1:Ntime) = q_flux(1:Ntime,1)   ! Monthly timeseries
        ! Moving average:
        Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d)
        q_flux(M2+1:Ntime-M2,2) = Y1d(1:Nt)

        ! Anomaly:
        Do imon=1,Nmon
           Nc=0; tmp0=0; Nc1=0; tmp1=0
           Do it=1,Ntime
              imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon
              if (imon1==imon) then
                 if (Y1d(it)/=Miss) then 
                    Nc=Nc+1; tmp0=tmp0+Y1d(it)
                 endif
                 if (X1d(it)/=Miss) then
                    Nc1=Nc1+1; tmp1=tmp1+X1d(it)
                 endif
              endif
           Enddo
           Xtb(imon)=tmp0/real(Nc)
           Ytb(imon)=tmp1/real(Nc1)
        Enddo
        Do it=1,Ntime
           imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon
           if (Y1d(it)/=Miss) then
              q_flux(it,3)=Y1d(it)-Xtb(imon1)
           endif
           if (X1d(it)/=Miss) then
              q_flux(it,4)=X1d(it)-Ytb(imon1)
           endif
        Enddo
        ! Moving average:
        !Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d)
        !q_flux(M2+1:Ntime-M2,2) = Y1d(1:Nt)
        !Call wr_ncdf_1d (Fname, xtime, Ntime, q_flux, "qflux", &
        !     "Moisture fluxes", "kg/(m^2s)","time", Miss)
        Call wr_ncdf_2d (Fname,xtime,tmp,Ntime,4,q_flux,"qflux", &
              "Moisture fluxes","kg/(m^2s)","time","tmp",Miss)
      End Subroutine Moisture_Flux
      !---------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m,Y)
         Integer :: N,m, i,j,k   ! m must be an odd number
         Real :: X(N), Y(N), tmp ! Y(N-m+1), tmp

         k=m/2; Y=Miss
         Do j=1,N-m+1
            tmp=0
            Do i=j,j+m-1
               tmp=tmp+X(i)
            Enddo
            Y(j+k)=tmp/Real(m)
         Enddo
      End Subroutine Moving_Ave
      !---------------------------------------------------------------------------------
      Logical Function Empty_Ar1(X,N)
         Integer :: N,i
         Real :: X(N)
         Logical :: Empty
         Empty=.False.
         Do i=1,N
            if (X(i)==Miss.OR.IsNaN(X(i)).OR.ABS(X(i))>1000) then
               Empty=.True.; Exit
            endif
         Enddo
         Empty_Ar1=Empty
      End Function Empty_Ar1

      !----------------------------------------------------------------------------------!
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/Distance.f90"
      END
      !----------------------------------------------------------------------------------!