! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1961, Yr2=2015, Nidx=24 Integer, Parameter :: Nmode=6, Nreg=4 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", "pdo", "NTA_ersst", "dmi", & ! 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, 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(:,:,:) ! 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 :: eof_sst(:,:), Rsst(:,:),Pvsst(:,:) ! eof of sst and Corr between eofs of sst and Drou Index Real, Allocatable :: moist_flux(:,:), Rmoi(:,:), Pmoi(:,:) ! Moisture flux(Ntime,2): 1-month, 2-13month aver 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, Dir_Out !Character (Len=15) :: Dat_src="Ens" ! "APH" ! "Ens" Character (Len=15) :: Dat_src="APH" ! "Ens" !-------------------------------------------------------------------- Dir_Out="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/" !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_"//trim(Dat_src)//".nc" Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_APH_SEA_61-15.nc" !BL_APH_SEA.nc" Fname1= & ! "/home/tanpv/Progs/Dr_SEA/Cluster/Clus_BL_"//trim(Dat_src)// & ! "_spe12_06.nc" "/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/kmeans_4.nc" ! Calculated by R Call Read_data (Fname,Fname1,varn) Fname="/home/tanpv/Data/ClimIDX/Clim_IDX_2020.nc" Call Read_ClimIDX (Fname) !,varn) print*, "Calling Read_ClimIDX done !" Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/eof_clim.nc" Call Read_eof_clim (Fname) print*, "Calling Read_eof_clim done !" !if (trim(Dat_src)=="APH") then Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/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") print*, "Calling Read_eof_Drou done !" Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/eof_sst.nc" Call Read_eof_sst (Fname) print*, "Calling Read_eof_sst done !" ! Moisture Flux: Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/Moisture_flux.nc" Call Read_moist_flux (Fname) print*, "Calling Read_moist_flux done !" Nt=Ntime; Nx=Landpoint; Ny=Nidx Allocate (XX(Nt),YY(Nt)) Allocate (Rxy(Nx,Ny),Pval(Nx,Ny)) ! Corr between ClimIDX and Dr_idx Allocate (Rcl(Nx,Nmode),Pvcl(Nx,Nmode)) Allocate (Reof(Nmode,Nmode),Peof(Nmode,Nmode)) Allocate (Rsst(Nx,Nmode),Pvsst(Nx,Nmode)) ! Corr between sst and Dr_idx Allocate (Rmoi(Nx,Nmode),Pmoi(Nx,Nmode)) ! Corr between moi and Dr_idx 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)) !-------------------------------------------------------------------- ! Calculating corr matx between drou index and climIDX, and drou index and eof Clim, and Moi Flx: 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: Fname=trim(Dir_Out)//"Corr_Dr_Clim_"//trim(Dat_src)//"_"//trim(varn)//".nc" Fname1=trim(Dir_Out)//"Corr_Dr_EOF_Clim_"//trim(Dat_src)//"_"//trim(varn)//".nc" Call Write_Corr_Matrix (Fname,Fname1) Fname=trim(Dir_Out)//"Trend_"//trim(Dat_src)// & "_"//trim(varn)//".nc" Call Write_Trend (Fname,varn) Fname=trim(Dir_Out)//"Corr_eof_Dr_eof_clim_"//trim(Dat_src)//"_"// & trim(varn)//".txt" Call Corr_eof_Clim_and_eof_drou (Fname,varn) Fname=trim(Dir_Out)//"Areal_mean_and_Trend_"//trim(Dat_src)// & "_"//trim(varn)//".nc" Call Write_Mean_and_Trend_Reg (Fname,varn) !-------------------------------------------------------- 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) !if (X(j,it)<0) then ! Chi tinh nhung gridpoit co Dr<0 ! i1=i1+1; tmp1=tmp1+X(j,it) !endif endif if (Int(Reg(j))==ireg) then i=i+1; tmp=tmp+X(j,it) !if (X(j,it)<0) then ! Chi tinh nhung gridpoit co Dr<0 ! i=i+1; tmp=tmp+X(j,it) !endif 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 ! Correlation between drought index and eof of SST Do j=1,Nmode YY(1:Nt)=eof_sst(1:Nt,j) Call Correlation (XX,YY,Nt,ctmp,ptmp, Miss) Rsst(i,j)=ctmp; Pvsst(i,j)=ptmp Enddo ! Correlation between drought index and Moisture Fluxes: ! moist_flux(:,:), Rmoi(:,:), Pmoi(:,:) Do j=1,2 ! Monthly and 13-mon moving average YY(1:Nt)=moist_flux(1:Nt,j) Call Correlation (XX,YY,Nt,ctmp,ptmp, Miss) Rmoi(i,j)=ctmp; Pmoi(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) var_ln="p-value of Corr. between Drought IDX and EOF ClimIDX" Call wr_ncdf_3d_apd (Fname1,Y3D,"pval",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) !---------------------------------------------------------- !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) = Rsst(j,i); Y3D(ilon,ilat,i) = Pvsst(j,i) endif Enddo Enddo Enddo var_ln="Correlation between Drought IDX and EOF SST" Call wr_ncdf_3d_apd (Fname1,X3D,"rsst",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) var_ln="p-value of Correlation between Drought IDX and EOF SST" Call wr_ncdf_3d_apd (Fname1,Y3D,"psst",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) !---------------------------------------------------------- ! Cor Dr and moi: X3D=Miss; Y3D=Miss Do i=1,2 ! Nmode j=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)==1) then j=j+1 X3D(ilon,ilat,i) = Rmoi(j,i); Y3D(ilon,ilat,i) = Pmoi(j,i) endif Enddo Enddo Enddo var_ln="Correlation between Drought IDX and Moisture Fluxes:" & //" 1: Monthly, 2: 13-month moving average" Call wr_ncdf_3d_apd (Fname1,X3D,"rmoi",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) var_ln="p-value of Correlation between Drought IDX and moi flx" Call wr_ncdf_3d_apd (Fname1,Y3D,"pmoi",nlon,nlat,Nmode, var_ln,& "-","lon", "lat", "mod", Miss) Deallocate (X3D,Y3D,climidx) End Subroutine Write_Corr_Matrix !-------------------------------------------------------------------- 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 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") Call rd_ncdf_2d (Fname1,xlon,xlat,Nlo,Nla, Reg2D, "cluster") ! Calculated by R Where (Reg2D==-9999.) Reg2D=Miss ! Replate R's Missing values to Miss 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*,iidx," ",trim(Clim_nm(iidx)) Enddo End Subroutine Read_ClimIDX !------------------------------------------------------------------------------ Subroutine Read_moist_flux (Fname) Character (Len=*) :: Fname Real, Allocatable :: xmod(:),xtime(:) Integer :: Nmn,Nmod Call rd_ncdf_2d (Fname,xtime,xmod,Nmn,Nmod,moist_flux, "qflux") Deallocate(xmod,xtime) End Subroutine Read_moist_flux !------------------------------------------------------------------------------ Subroutine Read_eof_sst (Fname) Character (Len=*) :: Fname Real, Allocatable :: xmod(:),xtime(:) Integer :: Nmn,Nmod Call rd_ncdf_2d (Fname,xtime,xmod,Nmn,Nmod,eof_sst, "sst") Deallocate(xmod,xtime) End Subroutine Read_eof_sst !------------------------------------------------------------------------------ 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 "Sub_MCA.f90" !-------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90"