! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1981, Yr2=2020, Nidx=16 Integer, Parameter :: Nmode=6 !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40 Real, Parameter :: Slon1=70, Slon2=290, Slat1=-20, Slat2=40 Real, Allocatable :: Lon(:), Lat(:), xtime(:) 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 !---------------------------------------------------------------- Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! Use indices same as in J. Clim (2022) "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 "whwp", "pdo", "dmi", "tna", "amonus", & "pna", "censo", "soi", "reqsoi", "repac_slpa"/)!, & ! 6 ! "nino3", "nino34", "nino4", "oni", "censo", & ! "tni", & ! 6 ! "rindo_slpa","soi", "reqsoi", "repac_slpa", & ! "whwp", "ammsst", "NTA_ersst", "tna", "amonus", & ! "pdo", "dmi"/) ! "nino3", "nino34", "nino4", "oni", "censo", & ! "rindo_slpa", "soi", "reqsoi", "repac_slpa", "whwp", & ! "ammsst", "pna","tna", "amonus", "pdo", "dmi" /) ! Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 ! "whwp", "ammsst", "pdo", "NTA_ersst", "dmi", & ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi ! "tna", "amonus", & ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 ! "rindo_slpa", "qbo30", "wp", "ao", "nao"/) !"rindo_slpa", "tsa", "wp", "ao", "nao"/) !Character (Len=15) :: varn="ori_pdsi" !"spi_12" !Character (Len=15) :: varn="spei_12" ! 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 :: Nt, Nx, Ny, isys, ilen Integer :: Mov_Ave_Step ! Buoc truot trong Moving average Real :: Sumx,Corr,prob Real, Allocatable :: xx1d(:), yy1d(:) Logical :: Calc_mtx, Normalize, SDV_MatLab !-------------------------------------------------------------------- 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="spi" 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 !-------------------------------------------------------------------- Nyear = Yr2-Yr1+1 Ntime = Nyear*Nmon Mov_Ave_Step = 0 ! Smooth with moving average ilen = 3 ! 1=3, 2=6,3=12 Fname="/home/tanpv/Progs/Dr_4Books/Res/" & //"spi_spei_stn_1980-2020.nc" Call Read_Dr_IDX (Fname,varn,ilen) !Call Read_drou_data (Fname,varn,ilen) print*,"Reading dought index ... done !" Fname="/home/tanpv/Data/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" Call Read_ClimIDX (Fname) !,varn) Nt=Ntime; Nx=Nst; Ny=Nidx 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 = .true. SDV_MatLab = .true. Calc_mtx=.true. !-------------------------------------------------------------------- Call Centered_Data (Normalize) ! Su dung de tinh Ak, Bk Dir_Out="/home/tanpv/Progs/Dr_4Books/Res/" 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 (SDV_MatLab) then 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 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 R if (SDV_MatLab) then 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 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 ! print test whether U and V are orthogonal: print"(6F10.2)", MatMul(Transpose(Vk),Vk) endif !-------------------------------------------------------- ! Write to NC files: !Fnu=trim(Dir_Out)//"U_Dr_Cli_"//trim(Dat_src)//"_"//trim(varn)//".nc" Fnu=trim(Dir_Out)//"U_Dr_Cli_"//trim(varn) Fnv=trim(Dir_Out)//"V_Dr_Cli_"//trim(varn)//".nc" Fnab=trim(Dir_Out)//"AB_Dr_Cli_"//trim(varn)//".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_Cli_"//trim(Dat_src)//"_"//trim(varn) Fname1 = trim(Dir_Out)//"He_Cli_"//trim(Dat_src)//"_"//trim(varn) Fnu = trim(Dir_Out)//"pHo_Cli_"//trim(Dat_src)//"_"//trim(varn) Fnv = trim(Dir_Out)//"pHe_Cli_"//trim(Dat_src)//"_"//trim(varn) Call Write_Homo_Hete_Maps_to_NC (Fname,Fname1,Fnu,Fnv) !Fnu =trim(Dir_Out)//"Ho_He_Cli_"//trim(varn)//".nc" !Call Write_Homo_Hete_Maps_to_NC (Fnu) !-------------------------------------------------------- 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 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 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) !Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S") !print*, nlon,nlat,nx,ny !Do i=1,Ny; W(i) = UU(i,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 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 ! Do i=1,Nmode ! j=0 ! Do ilat=1,Nlat ! Do ilon=1,Nlon ! if (Mask(ilon,ilat)==1) then ! j=j+1 ! Homo(ilon,ilat,i)=Ho(j,i); p_Ho(ilon,ilat,i)=pHo(j,i) ! Hete(ilon,ilat,i)=He(j,i); p_He(ilon,ilat,i)=pHe(j,i) ! endif ! Enddo ! Enddo ! Enddo 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) !Allocate (xmod(Nmode)) !Do i=1,Nmode; xmod(i)=i; Enddo !t_unit="months since 1962-01-01 " !var_ln="Homogenerous correlation map r[Dr(t), a_k (t)]" !Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,Homo,"ho", & ! t_unit, var_ln, "-", "lon", "lat", "mod", Miss) !var_ln="p-values of Homo" !Call wr_ncdf_3d_apd (Fnu,p_Ho,"pho",nlon,nlat,Nmode, var_ln, & ! "-", "lon", "lat", "mod", Miss) !var_ln="Heterogenerous correlation map r[Dr(t), b_k (t)]" !Call wr_ncdf_3d_apd (Fnu,Hete,"he",nlon,nlat,Nmode, var_ln, & ! "-", "lon", "lat", "mod", Miss) !var_ln="p-values of Hete" !Call wr_ncdf_3d_apd (Fnu,p_He,"phe",nlon,nlat,Nmode, var_ln, & ! "-", "lon", "lat", "mod", Miss) 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(:), climidx(:) 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),climidx(Nidx)) Do i=1,Nmode; xmod(i)=i; Enddo Do i=1,Nt; xtime(i)=i; Enddo Do i=1,Nidx; climidx(i)=i; Enddo ! U matrix: !t_unit="months since 1962-01-01 " !var_ln="MCA for drought index and SST - U-matrix" !Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,U3D,"u", & ! t_unit, var_ln, "-", "lon", "lat", "mod", Miss) Call Write_GrADS_Stn_Multivar_1 (Uk,Nst,Nmode,Fnu) ! 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) ! V matrix: Do i=1,Nmode X2D(i,:)=Vk(:,i) Enddo var_ln="MCA for climate indices" Call wr_ncdf_2d (Fnv,xmod,climidx,Nmode,Nidx, X2D, "v", var_ln, & "-", "mod", "idx", Miss) !var_ln="Percentage of contribution of singular values" !Call wr_ncdf_1d_apd (Fnv,Percent,"per", Nmode, var_ln,"-","mod", Miss) !var_ln="Accumulated Percent. of contribution of singular values" !Call wr_ncdf_1d_apd (Fnv,Accum,"cum", Nmode, var_ln,"-","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 !Do imon=1,Nidx ! Var = X(ist,imon,i,iidx) Var = X(ist,i) WRITE (10) Var ! Gia tri cua bien !Enddo 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' !Do imon=1,Nidx !write (2001,'(A,I2.2,A)') trim(IDX_nm(iidx)),imon,' 0 99 Drough Chars' !Enddo 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_ClimIDX (Fname) !,varn) Character (Len=*) :: Fname !,varn Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k Allocate (Y(Nidx,Ntime)) Yr11=Yr1 ! Lay tu Yr1-Yr2 Do iidx=1,Nidx print*,Ntime,iidx,trim(Clim_nm(iidx)) Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) j=0 Do iyr=Yr11,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do k=1,Nmon j=j+1; Y(iidx,j)=tmp(i,k) Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) Enddo End Subroutine Read_ClimIDX !--------------------------------------------------------------------------------- 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_Dr_IDX (Fname,varn,ilen) Character (Len=*) :: Fname,varn Real, allocatable :: X3d(:,:,:),stnid(:),time(:),xlen(:) Integer :: Nt,Nle,it1 it1=Nmon+1 Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spi") !Allocate (Dr_IDX(Nst,Ntime,Nlen,Ndr_idx)) ! (Nst,Ntime,Nlen,NDr_idx) Where (IsNaN(X3d)) X3d=Miss Allocate (X(Nst,Ntime)) !,Y1d(Nt),X1d(Ntime)) X(1:Nst,1:Ntime) = X3d(1:Nst,it1:Nt,ilen) ! Chi lay timescale=ilen, tru Nmon dau tien print*,Nst,Ntime,it1,ilen; pause !Dr_IDX(1:Nst,1:Ntime,1:Nlen,1) = X3d(1:Nst,it1:Nt,1:Nlen) !Deallocate (stnid,time,xlen,X3d) !Call rd_ncdf_3d (Fname,stnid,time,xlen,Nst,Nt,Nle,X3d,"spei") !Where (IsNaN(X3d)) X3d=Miss !Dr_IDX(1:Nst,1:Ntime,1:Nlen,2) = X3d(1:Nst,it1:Nt,1:Nlen) Deallocate (stnid,time,xlen,X3d) !Where (Dr_IDX(:,:,:,2)==Miss) Dr_IDX(:,:,:,2)=Dr_IDX(:,:,:,1) Call rd_ncdf_1d (Fname,stnid,Nst,Lon,"lon") Deallocate (stnid) Call rd_ncdf_1d (Fname,stnid,Nst,Lat,"lat") Deallocate (stnid) ! Call rd_ncdf_1d (Fname,stnid,Nst,Reg,"reg") ! Deallocate (stnid) Call rd_ncdf_1d_Char (Fname,stnid,Nst,stname,"stname") Deallocate (stnid) print*,"Calling Read_Dr_IDX Done" End Subroutine Read_Dr_IDX !-------------------------------------------------------------------- Subroutine Read_drou_data (Fname,varn,ilen) Character (Len=*) :: Fname,varn Real, Allocatable :: SPI(:,:), tmp(:), XX(:,:,:),YY(:,:,:) Real, Allocatable :: X1d(:), Y1d(:), StnID(:), xmon(:), xlen(:) Integer :: M, ilen, ist, ipoint, nt, it, Nlen, imon, imon1 ! Read in drought index: print*,"Reading Drought..." Call rd_ncdf_3d (Fname,StnID,xmon,xlen,Nst,Nt,Nlen,XX,trim(varn)) Deallocate (StnID) Call rd_ncdf_1d (Fname, StnID, Nst, Lon, "lon") Deallocate (StnID) Call rd_ncdf_1d (Fname, StnID, Nst, Lat, "lat") Deallocate (StnID) Call rd_ncdf_1d_Char (Fname,StnID,Nst,Stname,"stname") Where (IsNaN(XX)); XX=Miss; Endwhere Ntime=Nt-ilen; it=ilen+1 Allocate (SPI(Nst,Ntime)) SPI(:,1:Ntime) = XX(:,it:Nt,ilen) ! Bo Nmon thang dau tien Allocate (tmp(Ntime)) Landpoint=Nst ! Moving Average: M = Mov_Ave_Step ! Buoc truot = 13 thang if (M>0) then Nt=Ntime-M+1 else Nt=Ntime endif Allocate (X(Landpoint,Nt),Y1d(Nt),X1d(Ntime)) !,Reg(Landpoint)) print*,Ntime,Nt,Nst Do ist=1, Nst X1d(1:Ntime)=SPI(ist,1:Ntime) if (M>0) then Call Moving_Ave (X1d,Ntime,M,Y1d) X(ist,1:Nt)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime) else X(ist,1:Nt)=X1d(1:Ntime) ! Nt=Ntime endif Enddo Ntime=Nt !print*,X; pause Deallocate (SPI,XX) End Subroutine Read_drou_data !-------------------------------------------------------------------------------------- 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"