! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1980, Yr2=2020, Nidx=24 Integer, Parameter :: Nmode=6 !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40 Real, Parameter :: Slon1=110, Slon2=290, Slat1=-30, Slat2=30 Real, Allocatable :: Lon(:), Lat(:), xtime(:), Regs(:) Real, Allocatable :: Lon_sst(:), Lat_sst(:), Mask_sst(:,:) Real, Allocatable :: Mask(:,:), xlon(:), xlat(:), ylon(:), ylat(:) Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j Integer :: SeaPoint, Nsst, Nlon_sst, Nlat_sst Integer :: Yr_sta Real :: Miss_sst !---------------------------------------------------------------- ! For MCA: Real, Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:) Real, Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:) Real, Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:) Real, Allocatable :: A(:,:), B(:,:) ! Timeseries of X & Y Real, Allocatable :: Uk(:,:), Vk(:,:), WWk(:,:), Wk(:) Real, Allocatable :: Ak(:,:), Bk(:,:), Corre(:), Probs(:) Real, Allocatable :: Accum(:), Percent(:) ! Accumulated percentage of Singular values Real, Allocatable :: XX(:), YY(:), Dr_IDX(:,:,:) ! (Nlon,Nlat,Ntime) Real, Allocatable :: U3D(:,:,:),V3D(:,:,:) ! Uk,Vk in 3D array Real, Allocatable :: Homo(:,:),Hete(:,:) ! Homo and Hetero maps (Nlon,Nlat,Nmode) Real, Allocatable :: p_Ho(:,:),p_He(:,:) ! Pval Integer :: ilen, Nt, Nx, Ny, isys Integer :: Mov_Ave_Step ! Buoc truot trong Moving average Real :: Sumx,Corr,prob Real, Allocatable :: xx1d(:), yy1d(:) Logical :: Calc_mtx, Normalize, SVD_MatLab, IsAnomaly !-------------------------------------------------------------------- Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab, Dir_Out !Character (Len=15) :: Dat_src="Ens" ! "APH" ! "Ens" Character (Len=15) :: Dat_src="Stn" ! "Ens" Character (Len=15) :: varn="HWD" Character (Len=3) :: Normalize_Opt Character (Len=20), Allocatable :: Stname(:) !-------------------------------------------------------------------- ! LAPACK: Integer :: M, N, LDA, LDU, LDVT, LWMAX, INFO, LWORK DOUBLE PRECISION, Allocatable :: AA(:,:), UU(:,:), VT(:,:), SS(:), WORK(:) EXTERNAL :: DGESVD !-------------------------------------------------------------------- ! Mov_Ave_Step = 0 ! Smooth with moving average !-------------------------------------------------------------------- ! print*,"Set options ..." ! Normalize = .false. ! Normalize data (true/false) ! Calc_mtx = .true. ! Calculate C-matrix (true/false) ! SVD_MatLab = .True. ! Use SVD procedure from MatLab (true) or R (false) ! IsAnomaly = .false. ! Di thuong = Origin - Monthly mean Include "MCA_Options.par" if (Normalize) then Normalize_Opt="Cor"; else Normalize_Opt="Cov"; endif print*,"Normalize,Calc_mtx,SVD_MatLab,IsAnomaly: ",Normalize,Calc_mtx,SVD_MatLab,IsAnomaly !-------------------------------------------------------------------- Dir_Out="/home/tanpv/Progs/HW_Ha/NC/" Fname="Heatwave.txt" Call Read_HW_data (Fname) print*,"Reading HW data ... done !" Fname="/work/projects/G09/SEACLID/SST/sst.mon.mean_COBE.nc" !HadISST_sst.nc" Yr_sta=1891; Miss_sst=1.e+20 ! Depend on sst source Call Read_sst_data (Fname,varn) print*,"Reading SST ... done !" Nt=Ntime; Nx=Landpoint; Ny=SeaPoint Nst=Landpoint !print*,"Ntime,Landpoint,SeaPoint=",Ntime,Landpoint,SeaPoint; pause Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(Ny), Sx(Nx), Sy(Ny)) !Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nt,Nx), B(Nt,Ny)) Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt)) Allocate (WW(Ny,Ny)) Allocate (Uk(Nx,Nmode), Vk(Ny,Nmode), WWk(Nmode,Nmode), Wk(Nmode)) !Allocate (Ak(Nt,Nmode), Bk(Nt,Nmode), Corre(Nmode),Probs(Nmode)) Allocate (Ak(Nmode,Nt), Bk(Nmode,Nt), Corre(Nmode),Probs(Nmode)) Allocate (Accum(Nmode),Percent(Nmode)) Allocate (U3D(Nlon,Nlat,Nmode),V3D(Nlon_sst,Nlat_sst,Nmode)) Allocate (XX(Nt),YY(Nt)) !-------------------------------------------------------------------- ! print*,"Centering data ..." ! Normalize = .false. ! Normalize data (true/false) ! Calc_mtxi = .true. ! Calculate C-matrix (true/false) ! SVD_MatLab = .True. ! Use SVD procedure from MatLab (true) or R (false) ! Anomalous = .true. ! Di thuong = Origin - Monthly mean ! if (Normalize) then ! Normalize_Opt="Cor"; ! else ! Normalize_Opt="Cov"; ! endif !-------------------------------------------------------------------- Call Centered_Data (Normalize) ! Su dung de tinh Ak, Bk if (Calc_mtx) then ! Chi tinh Ma tran Covariance ! Covariance matrix: print*,"Calculating Covariance matrix..." C = MatMul(Xc,Transpose(Yc))/Real(Nt-1) ! Write C to nc file: Allocate (xx1d(Nx), yy1d(Ny)) Do i=1,Nx; xx1d(i)=i; Enddo Do i=1,Ny; yy1d(i)=i; Enddo Call wr_ncdf_2d (trim(Dir_Out)//"C_Matrix.nc",xx1d,yy1d, & Nx,Ny,C,"cor","-", "-", "x", "y", Miss) Deallocate (xx1d,yy1d) ! Calculating SVD: if (SVD_MatLab) then ! MatLab isys = system ("./run_mlab.sh ") Fnu=trim(Dir_Out)//"U.nc" Fnv=trim(Dir_Out)//"V.nc" Fnab=trim(Dir_Out)//"S.nc" Call Read_UVW_MatLab (Fnu,Fnv,Fnab) else ! R isys = system ("./run_r.sh ") Fnu=trim(Dir_Out)//"u.nc" Fnv=trim(Dir_Out)//"v.nc" Fnab=trim(Dir_Out)//"w.nc" Call Read_UVW (Fnu,Fnv,Fnab) endif Call Calc_3D_AB !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ else ! Xu ky KQ sau khi thuc hien SVD bang MatLab hoac R if (SVD_MatLab) then ! Use MatLab run Fnu=trim(Dir_Out)//"U.nc" Fnv=trim(Dir_Out)//"V.nc" Fnab=trim(Dir_Out)//"S.nc" Call Read_UVW_MatLab (Fnu,Fnv,Fnab) else ! Use R run Fnu=trim(Dir_Out)//"u.nc" Fnv=trim(Dir_Out)//"v.nc" Fnab=trim(Dir_Out)//"w.nc" Call Read_UVW (Fnu,Fnv,Fnab) endif Call Calc_3D_AB endif !-------------------------------------------------------- ! Write to NC files: Fnu=trim(Dir_Out)//"U_HW_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnv=trim(Dir_Out)//"V_HW_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt)//".nc" Fnab=trim(Dir_Out)//"AB_HW_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt)//".nc" Call Write_to_NC (Fnu,Fnv,Fnab) ! Calculate and write Homogenerous and Heterogenerous correlation maps: Call Homo_And_Hetero_Map Fname = trim(Dir_Out)//"Ho_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fname1 = trim(Dir_Out)//"He_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnu = trim(Dir_Out)//"pHo_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnv = trim(Dir_Out)//"pHe_SST_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Call Write_Homo_Hete_Maps_to_NC (Fname,Fname1,Fnu,Fnv) !-------------------------------------------------------- print*,"Wk", Wk print*,"Corr=", Corre print*,"Prob=",Probs print*,"Perc=",Percent print*,"Accum",Accum CONTAINS !-------------------------------------------------------------------- Subroutine Read_UVW_MatLab (Fu,Fv,Fw) Character (Len=*) :: Fu,Fv,Fw Real, Allocatable :: UU(:,:) Real, Allocatable :: xlo(:), xla(:) Integer :: nlon,nlat,i,j Real, parameter :: Valex=9.9692100E+36 Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S") print*, nlon,nlat,nx,ny !Do i=1,Ny; W(i) = SQRT(UU(i,i)); Enddo Do i=1,Ny; W(i) = UU(i,i); Enddo Where (W==Valex) W=0 Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"U") print*, nlon,nlat,nx,ny Do i=1,Nx U(i,1:Ny) = UU(i,1:Ny) !*W(i) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"V") print*, nlon,nlat,nx,ny Do i=1,Ny V(i,1:Ny) = UU(i,1:Ny) !*W(i) Enddo Deallocate (xlo,xla,UU) End Subroutine Read_UVW_MatLab !-------------------------------------------------------------------- Subroutine Read_UVW (Fu,Fv,Fw) Character (Len=*) :: Fu,Fv,Fw Real, Allocatable :: UU(:,:) Real, Allocatable :: xlo(:), xla(:) Integer :: nlon,nlat,i,j Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny Do i=1,Nx U(i,1:Ny) = UU(i,1:Ny) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny Do i=1,Ny V(i,1:Ny) = UU(1:Ny,i) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny W(1:Ny) = UU(1:Ny,1) Deallocate (xlo,xla,UU) End Subroutine Read_UVW !-------------------------------------------------------------------- Subroutine Calc_3D_AB ! Sign inversion of U and V: !U = -U; V = -V Wk(1:Nmode) = W(1:Nmode) ! Singular values Do i=1,Nmode Uk(:,i)=U(:,i); Vk(:,i)=V(:,i) Enddo ! Singular matrix: WW=0; Do j=1,Ny; WW(j,j)=W(j); Enddo ! Full WWk=0; Do j=1,Nmode; WWk(j,j)=W(j); Enddo ! Nmode ! Timeseries of Nmode Ak = MatMul(Transpose(Uk),Xc) Bk = MatMul(Transpose(Vk),Yc) ! Correlations: Do i=1,Nmode XX(1:Nt)=Ak(i,1:Nt); YY(1:Nt)=Bk(i,1:Nt) Call Correlation (XX,YY,Nt,Corr,Prob, Miss) Corre(i)=Corr; Probs(i)=prob Enddo ! Accum (%) of singular values: Sumx=Sum(W(1:Ny)); Accum(:)=0 Do i=1,Nmode Percent(i)=W(i)/Sumx*100 Accum(i)=Sum(Percent(1:i)) Enddo ! Convert Uk to U3D: U3D=Miss Do i=1,Nmode j=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)==1) then j=j+1 U3D(ilon,ilat,i) = Uk(j,i) endif Enddo Enddo Enddo ! Convert Vk to V3D: V3D=Miss Do i=1,Nmode j=0 Do ilat=1,Nlat_sst Do ilon=1,Nlon_sst if (Mask_sst(ilon,ilat)==1) then j=j+1 V3D(ilon,ilat,i) = Vk(j,i) endif Enddo Enddo Enddo End Subroutine Calc_3D_AB !-------------------------------------------------------------------- Subroutine Homo_And_Hetero_Map ! Calculate Homogenerous correlation map = r[X(t), a_k (t)] ! and Heterogenerous correlation map = r[X(t), b_k (t)] Integer :: ilon,ilat,it,i,j Real :: ctmp, ptmp Real, Allocatable :: Ho(:,:),He(:,:),pHo(:,:),pHe(:,:),X1d(:), A1d(:), B1d(:) Allocate (X1d(Nt),A1d(Nt),B1d(Nt),Ho(Nx,Nmode),He(Nx,Nmode)) Allocate (pHo(Nx,Nmode),pHe(Nx,Nmode)) Do i=1,Nmode A1d(1:Nt)=Ak(i,1:Nt); B1d(1:Nt)=Bk(i,1:Nt) Do j=1,Nx X1d(1:Nt)=X(j,1:Nt) Call Correlation(X1d,A1d,Nt,ctmp,ptmp,Miss) Ho(j,i)=ctmp; pHo(j,i)=ptmp Call Correlation(X1d,B1d,Nt,ctmp,ptmp,Miss) He(j,i)=ctmp; pHe(j,i)=ptmp Enddo Enddo ! Convert to 3D: Allocate (Homo(Nx,Nmode),p_Ho(Nx,Nmode)) Allocate (Hete(Nx,Nmode),p_He(Nx,Nmode)) Homo=Ho; Hete=He; p_Ho=pHo; p_He=pHe End Subroutine Homo_And_Hetero_Map !-------------------------------------------------------------------- Subroutine Write_Homo_Hete_Maps_to_NC (FHo,FHe,FpHo,FpHe) Character (Len=*) :: FHo,FHe,FpHo,FpHe Real, Allocatable :: xmod(:) Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Integer :: i,j,k Call Write_GrADS_Stn_Multivar_1 (Homo,Nst,Nmode,FHo) Call Write_GrADS_Stn_Multivar_1 (Hete,Nst,Nmode,FHe) Call Write_GrADS_Stn_Multivar_1 (p_Ho,Nst,Nmode,FpHo) Call Write_GrADS_Stn_Multivar_1 (p_He,Nst,Nmode,FpHe) Call Write_to_TXT (trim(FHo)//".txt",Homo,Nst,Nmode,"Homo") Call Write_to_TXT (trim(FHe)//".txt",Hete,Nst,Nmode,"Hete") Call Write_to_TXT (trim(FpHo)//".txt",p_Ho,Nst,Nmode,"p-val_Ho") Call Write_to_TXT (trim(FpHe)//".txt",p_He,Nst,Nmode,"p-val_He") End Subroutine Write_Homo_Hete_Maps_to_NC !-------------------------------------------------------------------- Subroutine Write_to_NC (Fnu,Fnv,Fnab) Character (Len=*) :: Fnu, Fnv, Fnab ! Filename of U, V and A,B Real, Allocatable :: xmod(:), X2D(:,:), xtime(:) Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Integer :: i,j,k Allocate (xmod(Nmode),X2D(Nmode,Ny),xtime(Nt)) Do i=1,Nmode; xmod(i)=i; Enddo Do i=1,Nt; xtime(i)=i; Enddo ! U matrix: t_unit="months since 1962-01-01 " Call Write_GrADS_Stn_Multivar_1 (Uk,Nst,Nmode,Fnu) Call Write_to_TXT (trim(Fnu)//".txt",Uk,Nst,Nmode,"Uk") ! V matrix: var_ln="MCA for drought index and SST - V-matrix" Call wr_ncdf_3d (Fnv,Lon_sst,Lat_sst,xmod,nlon_sst,nlat_sst,Nmode,V3D,"v", & t_unit, var_ln, "-", "lon", "lat", "mod", Miss) ! A & B matrix var_ln="Timeseries of MCA of drought index (Ak)" Call wr_ncdf_2d (Fnab,xmod,xtime,Nmode,Nt,Ak, "a", var_ln, & "-", "mod", "mon", Miss) var_ln="Timeseries of MCA of SST (Bk)" Call wr_ncdf_2d_apd (Fnab,Bk,"b",Nmode,Nt, var_ln, "-", & "mod", "mon", Miss) var_ln="Correlation coefficients between A & B" Call wr_ncdf_1d_apd (Fnab,Corre,"cor", Nmode, var_ln,"-","mod", Miss) var_ln="p-values of Correlation coefficients between A & B" Call wr_ncdf_1d_apd (Fnab,probs,"prob", Nmode, var_ln,"-","mod", Miss) var_ln="Percentage of contribution of singular values" Call wr_ncdf_1d_apd (Fnab,Percent,"per", Nmode, var_ln,"-","mod", Miss) var_ln="Accumulated Percent. of contribution of singular values" Call wr_ncdf_1d_apd (Fnab,Accum,"cum", Nmode, var_ln,"-","mod", Miss) End Subroutine Write_to_NC !--------------------------------------------------------------------- Subroutine Write_GrADS_Stn_Multivar_1 (X,Nst,Nmode,Grads_File) ! (Nst,NLen,ndr_idx) Integer :: imon, ivar, ist, ic, i, iidx,Nst,Nmode Character (Len=8) :: STID Character *(*) :: Grads_File Real :: X(Nst,Nmode) !Nst,Nmon,NLen,ndr_idx) ! (Nst,2,Nvar) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG !---------------------------------------------------------------------------------- Open (10,file=trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') TIM = 0.0 ! Co the bien doi tu -0.5 .. +0.5; thuong dat bang 0. NFLAG = 1 ! Khong doi Do i = 1,Nmode ! Nvar chi so, moi chi so duoc xem nhu 1 timeslice NLEV = 1 ! =0 o ban ghi cuoi cung cua 1 timeslice; =1 cho tat ca cac timeslice khac ic = 0 Do ist=1, Nst ic = ic+1 STID = Stname(ist)(1:8) ! Chi lay 8 ky tu dau tien RLAT = lat(ist); RLON = lon(ist) WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Var = X(ist,i) WRITE (10) Var ! Gia tri cua bien Enddo NLEV = 0 WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Enddo ! ivar Close(10) !---------------------------------------------------------------------! ! Write control file Station data ! !---------------------------------------------------------------------! open (2001,file=trim(Grads_File)//'.ctl') write (2001,'(2A)') 'dset ','^'//trim(Grads_File)//'.dat' write (2001,'(A)') 'dtype station' write (2001,'(2A)') 'stnmap ','^'//trim(Grads_File)//'.dat'//'.map' write (2001,'(A,F8.1)') 'undef ',Miss write (2001,'(A)') 'title Factor Analysis' write (2001,'(A,I3,1x,A)') 'tdef ',Nmode,' linear jan1980 1mo' write (2001,'(A,I3)') 'vars 1' ! ,Nidx write (2001,"(A)") 'uvar 0 99 Drough Chars' write (2001,'(A)') 'endvars' close (2001) End Subroutine Write_GrADS_Stn_Multivar_1 !-------------------------------------------------------------------- Subroutine Correlation(X,Y,N,Corr,prob,Miss) Integer :: N, N1, i, j Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0 Do i=1,N if (X(i)/=Miss.AND.Y(i)/=Miss) then n1=n1+1 Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i) Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i) Rxy=Rxy+X(i)*Y(i) endif Enddo if (N1>2) then Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1) Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb tmp=Rxy/Real(N1) Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy) t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) ) Prob = T_Dist (t, Real(N1-2)) else Rxy=Miss; prob=Miss endif Corr=Rxy End Subroutine Correlation !--------------------------------------------------------------------------- Subroutine Centered_Data (Normalize) Integer :: i,j Logical :: Normalize Real, Allocatable :: Std(:), X1d(:) ! Centering: Allocate (Std(Nx), X1d(Nt)) Do j=1,Nx X1d(1:Nt)=X(j,1:Nt) Xbar(j)=Ave (X1d,Nt,Miss) Std (j)=Stdev_1(X1d,Nt,Miss) Do i=1,Nt Xc(j,i)=X(j,i)-Xbar(j) if (Normalize) Xc(j,i)=Xc(j,i)/Std(j) !/Stdev_1(X(j,1:Nt),Nt,Miss) Enddo Enddo Deallocate (Std); Allocate (Std(Ny)) Do j=1,Ny X1d(1:Nt)=Y(j,1:Nt) Ybar(j)=Ave (X1d,Nt,Miss) Std (j)=Stdev_1(X1d,Nt,Miss) Do i=1,Nt Yc(j,i)=Y(j,i)-Ybar(j) if (Normalize) Yc(j,i)=Yc(j,i)/Std(j) ! Stdev_1(Y(j,1:Nt),Nt,Miss) Enddo Enddo End Subroutine Centered_Data !-------------------------------------------------------------------- Subroutine Read_sst_data (Fname,varn) Character (Len=*) :: Fname,varn !Integer, Parameter :: Yr_sta=1891 ! Nam bat dau chuoi SL SST Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:) Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:) Real :: Xtb(Nmon) Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat, imon, imon1 Integer :: M, Md2, it1, it2, ilon1, ilon2,ilat1,ilat2 Character (Len=255) :: Fn ! Lon_sst(:), Lat_sst(:), Mask_sst(:,:), SeaPoint, Nsst, Nlon_sst, Nlat_sst Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,XX,"sst") !trim(varn)) Where (IsNaN(XX)); XX=Miss; Endwhere Where (XX==Miss_sst); XX=Miss; Endwhere i=1; Do while (tlon(i)Slat1); j=j+1; Enddo; ilat1=j j=1; Do while (tlat(j)>Slat2); j=j+1; Enddo; ilat2=j it1=(Yr1-Yr_sta)*Nmon+1 !+ Nmon ! Bo nam dau tien it2=it1+Ntime-1 ! Lay them M/2 phan tu dau va M/2 phan tu cuoi cho Moving Average: M = Mov_Ave_Step Md2=M/2 it1=it1-Md2; it2=it2+Md2 Nt = Ntime ! Luu gia tri Ntime Ntime = it2-it1+1 ! Kich thuoc chuoi moi dung trong Moving Average Nlon_sst=ilon2-ilon1+1; Nlat_sst=ilat1-ilat2+1 ! Chu y: ilat1>ilat2 vi Latitude giam tu North to South Allocate (YY(Nlon_sst,Nlat_sst,Ntime),Lon_sst(Nlon_sst), Lat_sst(Nlat_sst)) Allocate (Mask_sst(Nlon_sst,Nlat_sst), X1d(Ntime)) Mask_sst(:,:)=Miss X1d(1:Ntime) = xtime(it1:it2) i=0 Do ilon=1,Nlon if (ilon>=ilon1.AND.ilon<=ilon2) then i=i+1; Lon_sst(i)=tlon(ilon) j=Nlat_sst+1 Do ilat=1,Nlat if (ilat>=ilat2.AND.ilat<=ilat1) then j=j-1; Lat_sst(j)=tlat(ilat) YY(i,j,1:Ntime) = XX(ilon,ilat,it1:it2) endif Enddo endif Enddo ! Processing for MCA: SeaPoint=0 Do i=1,Nlon_sst Do j=1,Nlat_sst X1d(1:Ntime) = YY(i,j,1:Ntime) if (.not.Ar1_Empty(X1d,Ntime)) then SeaPoint=SeaPoint+1; Mask_sst(i,j)=1 endif Enddo Enddo ! Moving Average: M = Mov_Ave_Step ! Buoc truot = 13 thang if (M>0) then Nt=Ntime-M+1 else Nt=Ntime ! Nt=(Yr2-Yr1+1)*Nmon endif Allocate (Y(SeaPoint,Nt),Y1d(Nt),tmp(Nt)) ipoint=0 Do ilat=1, Nlat_sst Do ilon=1, Nlon_sst if (Mask_sst(ilon,ilat)/=Miss) then ipoint=ipoint+1 X1d(1:Ntime)=YY(ilon,ilat,1:Ntime) if (IsAnomaly.AND.M>0) then Do i=1,Nt ! Trich doan thoi gian tu Jan/Yr1 den Dec/Yr2 tmp(i)=X1d(i+Md2) Enddo Call Anomaly (tmp,Nt,Xtb) ! Anomaly = Origin - Monthly mean Do i=1,Nt; X1d(i+Md2)=tmp(i); Enddo ! Gan lai cho X1d Do i=Md2,1,-1; X1d(i)=X1d(i)-Xtb(i+Md2); Enddo ! Doan Md2 phan tu dau cua chuoi Do i=1,Md2; X1d(Nt+i)=X1d(Nt+i)-Xtb(i); enddo ! Doan Md2 phan tu cuoi cua chuoi else Call Anomaly (X1d,Nt,Xtb) ! Nt=Ntime endif ! Moving average: if (M>0) then Call Moving_Ave (X1d,Ntime,M,Y1d) Y(ipoint,1:Nt)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime) else Y(ipoint,1:Nt)=X1d(1:Ntime) ! Nt=Ntime endif endif Enddo Enddo Ntime = Nt End Subroutine Read_sst_data !--------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N-m+1), tmp Do j=1,N-m+1 tmp=0 Do i=j,j+m-1 tmp=tmp+X(i) Enddo Y(j)=tmp/Real(m) Enddo End Subroutine Moving_Ave !-------------------------------------------------------------------- Subroutine Read_HW_data (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nst=102, Nyrs=Yr2-Yr1+1 Integer :: it, ist, iyr, imon, ireg, id, yr Real :: X3d(Nst,Nyrs,Nmon), X1d(Nyrs*Nmon), Xtb(Nmon), xlon,xlat Allocate (Lon(Nst),Lat(Nst),Regs(Nst),Stname(Nst)) Open (1,file=Fname,Status="Old") Read (1,*) ! Tieu de Do ist=1,Nst Do iyr=1,Nyrs Read (1,*) id, yr, (X3d(ist,iyr,imon),imon=1,Nmon), ireg, xlon, xlat if (iyr==1) then; Regs(ist)=ireg; Lon(ist)=xlon; Lat(ist)=xlat; endif Enddo !Regs(ist)=ireg; Lon(ist)=xlon; Lat(ist)=xlat Write(Stname(ist),"(A,I3.3)") "Stn_",ist Enddo Close (1) Landpoint=Nst Nt=Nyrs*Nmon; Ntime=Nt Allocate (X(Landpoint,Nt)) Do ist=1,Nst it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1; X(ist,it)=X3d(ist,iyr,imon) Enddo Enddo Enddo ! Anomaly: if (IsAnomaly) then Do ist=1,Nst X1d(:)=X(ist,:) Call Anomaly (X1d,Nt, Xtb) X(ist,:)=X1d(:) Enddo endif End Subroutine Read_HW_data !-------------------------------------------------------------------------------------- Subroutine Anomaly (X,N,Xtb) Integer :: N, it, j, imon, ic Real :: X(N), Xtb(Nmon), tmp Do imon=1,Nmon tmp=0; ic=0 Do it=1,N j=mod(it,Nmon); if (j==0) j=Nmon if (j==imon) then; ic=ic+1; tmp=tmp+X(it); endif Enddo Xtb(imon)=tmp/ic Enddo Do it=1,N j=mod(it,Nmon); if (j==0) j=Nmon X(it)=X(it)-Xtb(j) Enddo End Subroutine Anomaly !-------------------------------------------------------------------------------------- Subroutine Write_to_TXT (Fname,X,N,M,Varname) Character (Len=*) :: Fname,Varname Integer :: N,M,i,j Real :: X(N,M) print*,"Write_to_TXT X ",N,M,"--->", trim(Fname) Open (3,file=Fname,Status="Unknown") Write (3,"(A)") "Varname= "//trim(Varname) Do i=1,N Write (3,"(I4,50F10.4)") i,(X(i,j),j=1,M) ! print*,i Enddo Close (3) print*,"Writing ",trim(Varname)," ... Done !" End Subroutine Write_to_TXT !-------------------------------------------------------------------------------------- Logical Function Ar1_Empty(X,N) Integer :: N, i Real :: X(N) Logical :: Empty Empty=.false. Do i=1,N if (X(i)==Miss) then Empty=.true.; exit endif Enddo Ar1_Empty=Empty End Function Ar1_Empty !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Stdev_1.f90" !Include "Sub_MCA.f90" !-------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90"