! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF USE nrtype USE nr Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1961, Yr2=2005, Nidx=24 Integer, Parameter :: Nmode=6, Nreg=6 Real, Allocatable :: Lon(:), Lat(:), xtime(:) Real, Allocatable :: Mask(:,:), xlon(:), xlat(:) Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 ! "pna", "censo", "ammsst", "CAR_ersst", "NTA_ersst", & ! 5 ! "soi", "reqsoi", "whwp", "amonus", "np", "repac_slpa", &! 6 ! "rindo_slpa", "tsa", "tna", & ! 3 ! "qbo30", "wp", "ao", "nao"/) ! 4 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 "whwp", "ammsst", "CAR_ersst", "NTA_ersst", "tsa", & ! 5 "tna", "amonus", & "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 "rindo_slpa", "qbo30", "wp", "ao", "nao"/) ! "pna", "censo", "ammsst", "CAR_ersst", "NTA_ersst", & ! 5 ! "soi", "reqsoi", "whwp", "amonus", "np", "repac_slpa", &! 6 ! "rindo_slpa", "tsa", "tna", & ! 3 ! "qbo30", "wp", "ao", "nao"/) ! 4 Character (Len=200) :: Clim_ln(Nidx)=(/ & "nino12: Extreme Eastern Tropical Pacific SST (0-10S, 90W-80W)", & "nino3: Eastern Tropical Pacific SST (5N-5S,150W-90W)", & "nino34: East Central Tropical Pacific SST (5N-5S)(170-120W)", & "nino4: Central Tropical Pacific SST *(5N-5S) (160E-150W)", & "oni: Oceanic Nino Index", & "tni: Indices of El Nino evolution (gradient in SST anomalies between the central and eastern equatorial Pacific)", & "whwp: Western Hemisphere warm pool (SST>28.5C", & "ammsst: Atlantic Meridional Mode (AMM) SST Index", & "CAR_ersst: Caribbean SST Index: SST anomalies averaged over the the Caribbean", & "NTA_ersst: North Tropical Atlantic SST Index", & "tsa: Tropical Southern Atlantic Index", & "tna: Tropical Northern Atlantic Index", & "amonus: AMO (Atlantic Multidecadal Oscillation) Index", & "pna: Pacific North American Index (500 hPa)", & "BEST (censo): Bivariate ENSO Timeseries", & "soi: Southern Oscillation Index", & "reqsoi: Equatorial SOI", & "np: North Pacific pattern is the area-weighted sea "// & "level pressure over the region 30N-65N, 160E-140W", & "repac_slpa: Equatorial Eastern Pacific SLP "// & "(Standardized Anomalies)", & "rindo_slpa: Indonesia SLP (Standardized Anomalies)", & "qbo30: Quasi-Biennial Oscillation", & "wp: Western Pacific Index (a primary mode of "// & "low-frequency variability)", & "ao: The loading pattern of AO (AAO) at 1000 (NH) or 700hPa (SH)", & "nao: North Atlantic Oscillation"/) !---------------------------------------------------------------- !Character (Len=15) :: varn="ori_pdsi" !"spi_12" Character (Len=15) :: varn="spei_12" ! For MCA: Real (SP), Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:) Real (SP), Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:) Real (SP), Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:) Real (SP), 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(:,:,:) ! Uk in 3D array Integer :: Nt, Nx, Ny Real :: Sumx,Corr,prob !-------------------------------------------------------------------- ! Correlation matrix between Drought index and ClimIDX: Real, Allocatable :: Rxy(:,:), Pval(:,:) ! (Nst,Nidx) Real, Allocatable :: eof_cli(:,:), Rcl(:,:), Pvcl(:,:) ! EOF Clim Real, Allocatable :: eof_drou(:,:) ! eof of drought Real, Allocatable :: Reof(:,:),Peof(:,:) ! Corr between eofs of Clim and Drou Index Real, Allocatable :: Sensl(:),PvSen(:) ! Average and Trend of drought on sub-reg Real, Allocatable :: Reg(:), Reg2D(:,:) ! Mask for subregs Real, Allocatable :: Drou_Mean(:,:) ! (Nreg,Nx) Real, Allocatable :: SenReg(:), Preg(:) ! (Nreg): Trend !-------------------------------------------------------------------- ! Trend of T2m and Prec: Real, Allocatable :: T2m(:,:,:),Prec(:,:,:) Real, Allocatable :: SenT2m(:,:), SenPrec(:,:),ProT2m(:,:),ProPrec(:,:) Integer :: Nyear Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab !Character (Len=15) :: Dat_src="Ens" ! "APH" ! "Ens" Character (Len=15) :: Dat_src="APH" ! "Ens" Real, Allocatable :: xx1d(:), yy1d(:) !-------------------------------------------------------------------- !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_"//trim(Dat_src)//".nc" Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_"//trim(Dat_src)//".nc" Fname1= & "/home/tanpv/Progs/Dr_SEA/Cluster/Clus_BL_"//trim(Dat_src)// & "_spei_12_06.nc" Call Read_data (Fname,Fname1,varn) Fname="/home/tanpv/Data/ClimIDX/Clim_IDX_2020.nc" Call Read_ClimIDX (Fname) !,varn) Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/eof_clim.nc" Call Read_eof_clim (Fname) if (trim(Dat_src)=="APH") then Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/eof_APH.nc" else if (trim(Dat_src)=="Ens") then Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/eof_BL_Ens.nc" endif Call Read_eof_Drou (Fname,"spe12") Nt=Ntime; Nx=Landpoint; Ny=Nidx Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(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)) Allocate (XX(Nt),YY(Nt)) Allocate (Rxy(Nx,Ny),Pval(Nx,Ny)) Allocate (Rcl(Nx,Nmode),Pvcl(Nx,Nmode)) Allocate (Reof(Nmode,Nmode),Peof(Nmode,Nmode)) Allocate (Drou_Mean(Nreg+1,Nt),SenReg(Nreg+1),Preg(Nreg+1)) ! Trends of drought index: Allocate (Sensl(Nx),PvSen(Nx)) Allocate (SenT2m(Nlon,Nlat), SenPrec(Nlon,Nlat)) Allocate (ProT2m(Nlon,Nlat),ProPrec(Nlon,Nlat)) !-------------------------------------------------------------------- Call Centered_Data ! Covariance matrix: C = MatMul(Xc,Transpose(Yc))/Real(Nt-1) !print"(30F6.2)",C, "C"; pause !-------------------------------------------------------------------- ! 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 ("C_Matrix.nc",xx1d,yy1d,Nx,Ny,C,"cor", "-", "-", "x", "y", Miss) !Deallocate (xx1d,yy1d) !stop !-------------------------------------------------------------------- ! SVD: !Call svdcmp (C, Nx, Ny, U, W, V) !Call SVDsrt (U, Nx, Ny, V, W) ! copy original matrix into u U(:,:) = C(:,:) ! perform decomposition ! call svdcmp(u(1:m,1:n),w(1:n),v(1:n,1:n)) Call svdcmp (U(1:Nx,1:Ny),W(1:Ny),V(1:Ny,1:Ny)) Call SVDsrt(W,U,Ny) Call SVDsrt(W,V,Ny) !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(Xc),Uk) !Bk = MatMul(Transpose(Yc),Vk) Ak = MatMul(Transpose(Uk),Xc) Bk = MatMul(Transpose(Vk),Yc) ! Correlations: Do i=1,Nmode !XX(1:Nt)=Ak(1:Nt,i); YY(1:Nt)=Bk(1:Nt,i) 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 !Do j=1,i Accum(i)=Sum(Percent(1:i)) !Enddo 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 ! Calculating corr matx between drou index and climIDX, and drou index and eof Clim: !! Call Corr_Drou_Clim ! Calculating Trend of drought and T2m, Prec: !! Call Trend_drought_and_T_R ! Calculating areal mean and trend: !! Call Areal_Mean_of_Drought_and_Trend !-------------------------------------------------------- ! Write to NC files: Fnu="U_"//trim(Dat_src)//"_"//trim(varn)//".nc" Fnv="V_"//trim(Dat_src)//"_"//trim(varn)//".nc" Fnab="AB_"//trim(Dat_src)//"_"//trim(varn)//".nc" Call Write_to_NC (Fnu,Fnv,Fnab) Fname="Corr_Dr_Clim_"//trim(Dat_src)//"_"//trim(varn)//".nc" Fname1="Corr_Dr_EOF_Clim_"//trim(Dat_src)//"_"//trim(varn)//".nc" !! Call Write_Corr_Matrix (Fname,Fname1) Fname="Trend_"//trim(Dat_src)//"_"//trim(varn)//".nc" !! Call Write_Trend (Fname,varn) Fname="Corr_eof_Dr_eof_clim_"//trim(Dat_src)//"_"// & trim(varn)//".txt" !! Call Corr_eof_Clim_and_eof_drou (Fname,varn) Fname="Areal_mean_and_Trend_"//trim(Dat_src)// & "_"//trim(varn)//".nc" !! Call Write_Mean_and_Trend_Reg (Fname,varn) !-------------------------------------------------------- print*,"Wk", Wk print*,"Corr=", Corre print*,"Prob=",Probs print*,"Perc=",Percent print*,"Accum",Accum CONTAINS !-------------------------------------------------------------------- Subroutine Areal_Mean_of_Drought_and_Trend Integer :: i,i1,j,it,ireg,N1,NumYear Real :: tmp,tmp1, SenSlope, prob ! Drou_Mean(Nreg,Nt),SenReg(Nreg),Preg(Nreg) Do it=1,Nt i1=0;tmp1=0 Do ireg=1,Nreg i=0; tmp=0 Do j=1,Nx ! Search over all gridpoints for ireg-th subreg if (Int(Reg(j))>=1.AND.Int(Reg(j))<=Nreg) then i1=i1+1; tmp1=tmp1+X(j,it) endif if (Int(Reg(j))==ireg) then i=i+1; tmp=tmp+X(j,it) endif Enddo if (i>0) Drou_Mean(ireg,it)=tmp/Real(i) Enddo if (i1>0) Drou_Mean(Nreg+1,it)=tmp1/Real(i1) Enddo NumYear=30 ! Minimum sample Do ireg=1,Nreg+1 XX(1:Nt)=Drou_Mean(ireg,1:Nt) !! Call MK_Sen(XX,Nt,NumYear,Miss,SenSlope,prob,N1) SenReg(ireg)=SenSlope; Preg(ireg)=prob Enddo End Subroutine Areal_Mean_of_Drought_and_Trend !-------------------------------------------------------------------- Subroutine Write_Mean_and_Trend_Reg (Fname,varn) Character (Len=*) :: Fname, varn Real, Allocatable :: xtime(:),xreg(:) Integer :: i,j,ilon,ilat Character (Len=150) :: t_unit, var_ln Allocate (xtime(Nt),xreg(Nreg+1)) Do i=1,Nt; xtime(i)=i; Enddo Do i=1,Nreg+1; xreg(i)=i; Enddo var_ln="Areal mean of drought index over sub-regs" Call wr_ncdf_2d (Fname,xreg,xtime,Nreg+1,Nt,Drou_Mean,"dridx", & var_ln,"-", "reg", "mon", Miss) var_ln="Trend (Sen slopes) of the drought index ("//trim(varn)//")" Call wr_ncdf_1d_apd (Fname,SenReg,"sen", Nreg+1, var_ln, & "-","reg", Miss) var_ln="p-values of Sen slopes" Call wr_ncdf_1d_apd (Fname,Preg,"pval", Nreg+1, var_ln,"-","reg", Miss) Deallocate (xtime,xreg) End Subroutine Write_Mean_and_Trend_Reg !-------------------------------------------------------------------- Subroutine Trend_drought_and_T_R Integer :: i,j,NumYear,N1,ilon,ilat Real :: SenSlope, prob, X1d(Nyear), Ptb NumYear=30 ! Minimum sample print*,"Calculating Trend drought..." Do i=1,Nx XX(1:Nt)=X(i,1:Nt) !! Call MK_Sen(XX,Nt,NumYear,Miss,SenSlope,prob,N1) Sensl(i)=SenSlope; PvSen(i)=prob !print*,XX !print*,SenSlope,prob;pause Enddo print*,"Calculating Trend T & R..." SenT2m=Miss; SenPrec=Miss; ProT2m=Miss; ProPrec=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)>0) then X1d(1:Nyear)=T2m(ilon,ilat,1:Nyear) !! Call MK_Sen(X1d,Nyear,20,Miss,SenSlope,prob,N1) SenT2m(ilon,ilat)=SenSlope; ProT2m(ilon,ilat)=prob X1d(1:Nyear)=Prec(ilon,ilat,1:Nyear) Ptb=Ave(X1d,Nyear,Miss) !! Call MK_Sen(X1d,Nyear,20,Miss,SenSlope,prob,N1) SenPrec(ilon,ilat)=SenSlope/Ptb*100; ProPrec(ilon,ilat)=prob endif Enddo Enddo End Subroutine Trend_drought_and_T_R !-------------------------------------------------------------------- Subroutine Corr_eof_Clim_and_eof_drou (Fname,varn) Character (Len=*) :: Fname,varn Integer :: i,j Real :: ctmp, ptmp ! Correlation between eof of drought index and eof clim index Do i=1,Nmode XX(1:Nt)=eof_drou(1:Nt,i) Do j=1,Nmode YY(1:Nt)=eof_cli(1:Nt,j) !! Call Correlation (XX,YY,Nt,ctmp,ptmp, Miss) Reof(i,j)=ctmp; Peof(i,j)=ptmp !print*,"--->",ctmp, ptmp Enddo Enddo Open (3,file=Fname,status="Unknown") write (3,*) "Correlation between eof of "//trim(varn)// & " and eof of clim indices" write (3,"(2A8,20I8)") "Dr \"," Cli",(i,i=1,Nmode) Do i=1, Nmode write(3,"(A8,I8,20F8.3)") "R",i,(Reof(i,j),j=1,Nmode) !write(3,"(A8,I8,20F8.3)") "p",i,(Peof(i,j),j=1,Nmode) Enddo write (3,*) Do i=1, Nmode write(3,"(A8,I8,20F8.3)") "p",i,(Peof(i,j),j=1,Nmode) Enddo Close(3) End Subroutine Corr_eof_Clim_and_eof_drou !-------------------------------------------------------------------- Subroutine Corr_Drou_Clim Integer :: i,j Real :: ctmp, ptmp Do i=1,Nx ! Correlation between drought index and clim index XX(1:Nt)=X(i,1:Nt) Do j=1,Ny YY(1:Nt)=Y(j,1:Nt) !! Call Correlation (XX,YY,Nt,ctmp,ptmp, Miss) Rxy(i,j)=ctmp; Pval(i,j)=ptmp Enddo ! Correlation between drought index and eof of clim indices Do j=1,Nmode YY(1:Nt)=eof_cli(1:Nt,j) !! Call Correlation (XX,YY,Nt,ctmp,ptmp, Miss) Rcl(i,j)=ctmp; Pvcl(i,j)=ptmp Enddo Enddo End Subroutine Corr_Drou_Clim !-------------------------------------------------------------------- Subroutine Write_Trend (Fname,varn) Character (Len=*) :: Fname, varn Real, Allocatable :: X2D(:,:), Y2D(:,:) Integer :: i,j,ilon,ilat Character (Len=150) :: t_unit, var_ln Allocate (X2D(Nlon,Nlat),Y2D(Nlon,Nlat)) X2D=Miss; Y2D=Miss j=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)==1) then j=j+1; X2D(ilon,ilat)=Sensl(j); Y2D(ilon,ilat)=PvSen(j) endif Enddo Enddo var_ln="Trend (Sen slopes) of the drought index ("//trim(varn)//")" Call wr_ncdf_2d (Fname,xlon,xlat,Nlon,Nlat, X2D, "sen", var_ln, & "1/mon", "lon", "lat", Miss) var_ln="p-values of Sen slopes" Call wr_ncdf_2d_apd (Fname,Y2D,"pval",nlon,nlat,var_ln, & "-", "lon", "lat", Miss) var_ln="Trend (Sen slopes) of "//trim(Dat_src)//" 2m-temp" Call wr_ncdf_2d_apd (Fname, SenT2m,"sent",nlon,nlat,var_ln, & "C/year", "lon", "lat", Miss) var_ln="Trend (Sen slopes) of "//trim(Dat_src)//" 2m-temp" Call wr_ncdf_2d_apd (Fname, SenPrec,"senp",nlon,nlat,var_ln, & "%/year", "lon", "lat", Miss) var_ln="p-values of Sen slopes of Temp" Call wr_ncdf_2d_apd (Fname,ProT2m,"p_t",nlon,nlat,var_ln, & "-", "lon", "lat", Miss) var_ln="p-values of Sen slopes of Prec" Call wr_ncdf_2d_apd (Fname,ProT2m,"p_p",nlon,nlat,var_ln, & "-", "lon", "lat", Miss) Deallocate (X2D,Y2D) ! SenT2m=Miss; SenPrec=Miss; ProT2m=Miss; ProPrec=Miss End Subroutine Write_Trend !-------------------------------------------------------------------- Subroutine Write_Corr_Matrix (Fname,Fname1) Character (Len=*) :: Fname,Fname1 Real, Allocatable :: X3D(:,:,:), Y3D(:,:,:), climidx(:) Integer :: i,j,ilon,ilat Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Allocate (X3D(Nlon,Nlat,Ny), Y3D(Nlon,Nlat,Ny), climidx(Ny)) Do i=1,Ny; climidx(i)=i; Enddo ! Ny=Nidx X3D=Miss; Y3D=Miss Do i=1,Ny j=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)==1) then j=j+1 X3D(ilon,ilat,i) = Rxy(j,i); Y3D(ilon,ilat,i) = Pval(j,i) endif Enddo Enddo Enddo t_unit="months since 1977-01-01 " var_ln="Correlation between Drought IDX and ClimIDX, and p-valu" Call wr_ncdf_3d (Fname,xlon,xlat,climidx,nlon,nlat,Ny,X3D,"r", & t_unit, var_ln, "-", "lon", "lat", "idx", Miss) Call wr_ncdf_3d_apd (Fname,Y3D,"pval",nlon,nlat,Ny, var_ln, "-","lon", "lat", "idx", Miss) !---------------------------------------------------------- Deallocate (X3D,Y3D,climidx) Allocate (X3D(Nlon,Nlat,Nmode), Y3D(Nlon,Nlat,Nmode), climidx(Nmode)) !---------------------------------------------------------- Do i=1,Nmode; climidx(i)=i; Enddo X3D=Miss; Y3D=Miss Do i=1,Nmode j=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)==1) then j=j+1 X3D(ilon,ilat,i) = Rcl(j,i); Y3D(ilon,ilat,i) = Pvcl(j,i) endif Enddo Enddo Enddo var_ln="Correlation between Drought IDX and EOF ClimIDX" Call wr_ncdf_3d (Fname1,xlon,xlat,climidx,nlon,nlat,Nmode,X3D,& "r",t_unit, var_ln, "-", "lon", "lat", "mod", Miss) Call wr_ncdf_3d_apd (Fname1,Y3D,"pval",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) Deallocate (X3D,Y3D,climidx) End Subroutine Write_Corr_Matrix !-------------------------------------------------------------------- Subroutine Write_to_NC (Fnu,Fnv,Fnab) Character (Len=*) :: Fnu, Fnv, Fnab ! Filename of U, V and A,B Real, Allocatable :: xmod(:), X2D(:,:), 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),climidx(Nidx),xtime(Nt)) Do i=1,Nmode; xmod(i)=i; Enddo Do i=1,Nidx; climidx(i)=i; Enddo Do i=1,Nt; xtime(i)=i; Enddo ! U matrix: t_unit="months since 1977-01-01 " var_ln="MCA for drought index - U-matrix" Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,U3D,"u", & 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 climate indices (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) End Subroutine Write_to_NC !-------------------------------------------------------------------- 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 Integer :: i,j ! Centering: Do j=1,Nx Xbar(j)=Ave(X(j,1:Nt),Nt,Miss) Do i=1,Nt Xc(j,i)=X(j,i)-Xbar(j) Enddo Enddo Do j=1,Ny Ybar(j)=Ave(Y(j,1:Nt),Nt,Miss) Do i=1,Nt Yc(j,i)=Y(j,i)-Ybar(j) Enddo Enddo End Subroutine Centered_Data !-------------------------------------------------------------------- Subroutine Read_data (Fname,Fname1,varn) Character (Len=*) :: Fname,Fname1,varn Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:) Real, Allocatable :: X1d(:), Y1d(:) Integer :: ist, ilon, ilat, ipoint, nt, it, nlo, nla, imon, imon1 ! Read in subregions: Call rd_ncdf_2d (Fname1,xlon,xlat,Nlo,Nla, Reg2D, "reg") Deallocate (xlon,xlat) ! Read in T2m and Prec: print*,"Reading T2m..." Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,XX,"t2m") Deallocate (xlon,xlat,xtime) print*,"Reading Prec...",Nt Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,YY,"pr") Deallocate (xlon,xlat,xtime) ! Process T & P: Nyear = Nt/Nmon ! So nam Allocate (T2m(Nlo,Nla,Nyear),Prec(Nlo,Nla,Nyear)) Allocate (X1d(Nmon),Y1d(Nmon)) print*,"nlo,nla,nt,Nyear=",nlo,nla,nt,Nyear Where (IsNaN(XX)); XX=Miss; Endwhere Where (IsNaN(YY)); YY=Miss; Endwhere Do ilat=1,Nla Do ilon=1,Nlo Do it=1,Nyear !print*,ilat,ilon,it imon=(it-1)*Nmon+1; imon1=imon+Nmon-1 X1d(1:Nmon)=XX(ilon,ilat,imon:imon1) Y1d(1:Nmon)=YY(ilon,ilat,imon:imon1) T2m(ilon,ilat,it)=Ave(X1d,Nmon,Miss) Prec(ilon,ilat,it)=Ave(Y1d,Nmon,Miss)*Nmon Enddo Enddo Enddo Deallocate (XX,YY,X1d,Y1d) ! Read in drought index: print*,"Reading Drought..." Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,XX,trim(varn)) Where (IsNaN(XX)); XX=Miss; Endwhere Ntime=Nt-Nmon; it=Nmon+1 Allocate (PDSI(Nlon,Nlat,Ntime),Dr_IDX(Nlon,Nlat,Ntime)) PDSI(:,:,1:Ntime) = XX(:,:,it:Nt) ! Bo Nmon thang dau tien Nst=Nlon*Nlat Allocate (Mask(Nlon,Nlat)) Allocate (lon(Nst), lat(Nst), tmp(Ntime)) ist=0; Landpoint=0; Mask(:,:)=Miss Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1; Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat) tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime) if (.not.Ar1_Empty(tmp,Ntime).AND.Reg2D(ilon,ilat)/=Miss) then Landpoint=Landpoint+1; Mask(ilon,ilat)=1 endif Enddo Enddo Allocate (X(Landpoint,Ntime),Reg(Landpoint)) ipoint=0 Do ilat=1, Nlat Do ilon=1, Nlon !tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime) if (Mask(ilon,ilat)/=Miss) then !if (.not.Ar1_Empty(tmp,Ntime)) then ipoint=ipoint+1 Reg(ipoint)=Reg2D(ilon,ilat) ! 1D Reg X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime) Dr_IDX(ilon,ilat,1:Ntime)=PDSI(ilon,ilat,1:Ntime) ! For further use endif Enddo Enddo !print*,X; pause Deallocate (PDSI,XX,xtime) End Subroutine Read_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)) !if (trim(varn)=="spi_12".OR.trim(varn)=="spei_12") then Yr11=Yr1+1 !else ! Yr11=Yr1 !endif !print*,Nidx,Ntime,Yr1,Yr11,Yr2; pause Do iidx=1,Nidx 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 !j=j+1; Y(iidx=tmp(i,:) !print "(F5.0,12F8.3)",Year(i),ClimIDX(iidx,j,:) endif Enddo Enddo Deallocate (Year,Month,tmp) !print "(12F7.3)",Y(iidx,:) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !------------------------------------------------------------------------------ Subroutine Read_eof_clim (Fname) Character (Len=*) :: Fname Real, Allocatable :: xmod(:),xtime(:) Integer :: Nmn,Nmod Call rd_ncdf_2d (Fname,xtime,xmod,Nmn,Nmod,eof_cli, "eof") Deallocate(xmod,xtime) End Subroutine Read_eof_clim !------------------------------------------------------------------------------ Subroutine Read_eof_Drou (Fname,varn) Character (Len=*) :: Fname,varn Real, Allocatable :: xmod(:),xtime(:) Integer :: Nmn,Nmod Call rd_ncdf_2d (Fname,xtime,xmod,Nmn,Nmod,eof_drou, trim(varn)) Deallocate(xmod,xtime) End Subroutine Read_eof_Drou !------------------------------------------------------------------------------ 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/Progs/Dr_SEA/MCA/SUB90/NR.f90" !Include "/home/tanpv/Progs/Dr_SEA/MCA/SUB90/NrType.f90" !Include "/home/tanpv/Progs/Dr_SEA/MCA/SUB90/NrUtil.f90" Include "/home/tanpv/Progs/Dr_SEA/MCA/SUB90/SVDcmp.f90" Include "/home/tanpv/Progs/Dr_SEA/MCA/SUB90/SVDsrt.f90" !Include "Sub_MCA.f90" !-------------------------------------------------------------------- END !Include "/home/tanpv/Mylibs/Student_Prob.f90"