! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1960, Yr2=2019, Nidx=16 ! 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"/) Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & "nino3", "nino34", "nino4", "oni", "censo", & ! "tni", & ! 6 "rindo_slpa","soi", "reqsoi", "repac_slpa", & !"whwp", "ammsst", "pdo", "NTA_ersst", "dmi", & ! 5 Thay CAR_ersst bang pdo, thay tsa bang dmi "whwp", "ammsst", "NTA_ersst", "tna", "amonus", & "pdo", "dmi"/) ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 Character (Len=200) :: Clim_ln(Nidx)=(/ & "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", & "BEST (censo): Bivariate ENSO Timeseries", & "rindo_slpa: Indonesia SLP (Standardized Anomalies)", & "soi: Southern Oscillation Index", & "reqsoi: Equatorial SOI", & "repac_slpa: Equatorial Eastern Pacific SLP "// & "(Standardized Anomalies)", & "whwp: Western Hemisphere warm pool (SST>28.5C", & "ammsst: Atlantic Meridional Mode (AMM) SST Index", & "NTA_ersst: North Tropical Atlantic SST Index", & "tna: Tropical Northern Atlantic Index", & "amonus: AMO (Atlantic Multidecadal Oscillation) Index", & "pdo: Pacific Decadal Oscillation", & "DMI: HadISST1.1."/) ! "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" ! For MCA: !Real, Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:) Real, Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:) Real, Allocatable :: XX(:), YY(:), Dr_IDX(:,:,:) ! (Nlon,Nlat,Ntime) Real, Allocatable :: ClimIDX(:,:) ! (Nidx,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 :: Sensl(:),PvSen(:) ! Correlation matrix between ClimIDXs: Real, Allocatable :: RClim(:,:), PClim(:,:) ! (Nidx,Nidx) ! 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(:,:,:) ! Annual Real, Allocatable :: T2mon(:,:,:),Precmon(:,:,:) ! Monthly Real, Allocatable :: SenT2m(:,:), SenPrec(:,:),ProT2m(:,:),ProPrec(:,:) Integer :: Nyear Real, Allocatable :: Prec_ano(:,:),Cor_Dr_Pr_Reg(:),P_Dr_Pr_Reg(:) Real, Allocatable :: spe12(:,:) 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/2021/NC/" !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_"//trim(Dat_src)//".nc" Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/"// & "SPEI_SPI_GPCC_CRU_59_19_New.nc" Fname1="/home/tanpv/Progs/2021/NC/kmeans_4_modified.nc" ! Calculated by R Call Read_data (Fname,Fname1,varn) Fname="/work/projects/G09/SEACLID/ClimIDX/Clim_IDX_2020.nc" Call Read_ClimIDX (Fname) !,varn) Fname="/home/tanpv/Progs/2021/NC/eof_clim.nc" Call Read_eof_clim (Fname) !if (trim(Dat_src)=="APH") then Fname="/home/tanpv/Progs/2021/NC/eof_spei.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") Fname="/home/tanpv/Progs/2021/NC/eof_sst.nc" Call Read_eof_sst (Fname) 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 (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: Call Corr_Drou_Clim ! Calculating Trend of drought and T2m, Prec: Call Trend_drought_and_T_R ! Calculating areal mean and trend: Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/"// & "SPEI_SPI_GPCC_CRU_Reg_59_19_New.nc" Call Areal_Mean_of_Drought_and_Trend (Fname) Fname1=trim(Dir_Out)//"Dr_Prec_and_Corr_Revise.nc" Call Corr_Dr_Prec (Fname, Fname1) !-------------------------------------------------------- ! Write to NC files: Fname=trim(Dir_Out)//"Corr_Dr_Clim_Revise.nc" Fname1=trim(Dir_Out)//"Corr_Dr_EOF_Clim_Revise.nc" Call Write_Corr_Matrix (Fname,Fname1) Fname=trim(Dir_Out)//"Trend_spei_Revise.nc" Call Write_Trend (Fname,varn) Fname=trim(Dir_Out)//"Corr_eof_Dr_eof_clim_Revise.txt" Call Corr_eof_Clim_and_eof_drou (Fname,varn) Fname=trim(Dir_Out)//"Areal_mean_and_Trend_Revise.nc" Call Write_Mean_and_Trend_Reg (Fname,varn) Fname=trim(Dir_Out)//"Corr_ClimIDX-ClimIDX_Revise.txt" Call Corr_ClimIDX_ClimIDX (Fname) !-------------------------------------------------------- CONTAINS !-------------------------------------------------------------------- Subroutine Corr_Dr_Prec (Fname, Fname1) ! For Minor revision Character (Len=*) :: Fname,Fname1 Integer :: imon,imon1,i,i1,j,it,ireg,N1,NumYear,nt,nregs,nscal,M_truot Real :: tmp,tmp1,ntmp, cor, prob, tb Real, Allocatable :: xreg(:),xscal(:),xmon(:),X3d(:,:,:),Y3d(:,:,:) Real, Allocatable :: X1d(:),Y1d(:) Real, Allocatable :: R_ano(:,:),Corr(:),Pval(:) ! Prec_ano(:,:),Cor_Dr_Pr_Reg(:),P_Dr_Pr_Reg(:),spe12(:,:) M_truot=12 ! Trung binh truot 12 thang Call rd_ncdf_3d (Fname,xreg,xscal,xmon,nregs,nscal,nt,X3d,"pr") Deallocate (xreg,xscal,xmon) Call rd_ncdf_3d (Fname,xreg,xscal,xmon,nregs,nscal,nt,Y3d,"spe12") Allocate (Prec_ano(nregs,nt),spe12(nregs,nt)) Allocate (X1d(Nt),Y1d(Nt)) Allocate (Cor_Dr_Pr_Reg(Nregs),P_Dr_Pr_Reg(Nregs)) Allocate (R_ano(nregs,nt),Corr(Nregs),Pval(Nregs)) Prec_ano=Miss; Cor_Dr_Pr_Reg=Miss; P_Dr_Pr_Reg=Miss R_ano=Miss; Corr=Miss; Pval=Miss Do ireg=1,Nregs spe12(ireg,1:Nt)=Y3d(ireg,1,1:Nt) X1d(1:Nt)=X3d(ireg,1,1:Nt) ! Monthly Ano: Do imon=1,Nmon tmp=0; ntmp=0 Do it=1,Nt imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon.AND.X1d(it)/=Miss) then ntmp=ntmp+1; tmp=tmp+X1d(it) endif Enddo if (ntmp/=0) then; tmp=tmp/ntmp; else; tmp=Miss; endif ! TB cua thang imon Do it=1,Nt ! Ano cua thang imon imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon.AND.X1d(it)/=Miss.AND.tmp/=Miss) then R_ano(ireg,it)=X1d(it)-tmp endif Enddo Enddo Call Moving_Ave_12_mon (X1d,Nt,M_truot,Y1d) tmp=Ave(Y1d,Nt,Miss) Prec_ano(ireg,1:Nt)=Y1d(1:Nt)-tmp X1d(1:Nt)=Prec_ano(ireg,1:Nt) Y1d(1:Nt)=spe12(ireg,1:Nt) Call Correlation (X1d,Y1d,Nt,cor,prob,Miss) Cor_Dr_Pr_Reg(ireg)=cor P_Dr_Pr_Reg(ireg) =prob X1d(1:Nt)=R_ano(ireg,1:Nt) Y1d(1:Nt)=spe12(ireg,1:Nt) Call Correlation (X1d,Y1d,Nt,cor,prob,Miss) Corr(ireg)=cor; Pval(ireg)=prob Enddo ! Write to file: Call wr_ncdf_2d (Fname1,xreg,xmon,Nregs,Nt,spe12,"spe12", & "SPEI_12 drought index","-","reg","mon",Miss) Call wr_ncdf_2d_apd (Fname1,Prec_ano,"pr_ano",Nregs,Nt, & ! Ann Ano "12-mon moving aver rainfall anomaly","mm/mon","reg","mon",Miss) Call wr_ncdf_2d_apd (Fname1,R_ano,"r_ano",Nregs,Nt, & ! Monthly Ano "Monthly rainfall anomaly","mm/mon","reg","mon",Miss) Call wr_ncdf_1d_apd (Fname1,Cor_Dr_Pr_Reg,"cor",Nregs, & "Corr between Dr and Prec","-","reg",Miss) Call wr_ncdf_1d_apd (Fname1,P_Dr_Pr_Reg,"pval",Nregs, & "p-value of Corr between Dr and Prec","-","reg",Miss) Call wr_ncdf_1d_apd (Fname1,Corr,"cor1",Nregs, & "Corr between Dr and Prec","-","reg",Miss) Call wr_ncdf_1d_apd (Fname1,Pval,"pval1",Nregs, & "p-value of Corr between Dr and Prec","-","reg",Miss) End Subroutine Corr_Dr_Prec !--------------------------------------------------------------------------------- Subroutine Moving_Ave_12_mon (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N), tmp ! Trung binh m phan tu gan vao gia tri cuoi cua chuoi m phan tu do: 1:12->12; 2:13->13,...,N-m+1:N->N Y=Miss Do j=m,N ! 1,N-m+1 tmp=0 Do i=j-m+1,j ! j,j+m-1 tmp=tmp+X(i) Enddo Y(j)=tmp/Real(m) Enddo End Subroutine Moving_Ave_12_mon !-------------------------------------------------------------------- Subroutine Corr_ClimIDX_ClimIDX (Fname) Character (Len=*) :: Fname Integer :: i,j,k Real :: ctmp,ptmp Real, Allocatable :: X1d(:), Y1d(:) Allocate (X1d(Nt), Y1d(Nt),RClim(Nidx,Nidx),PClim(Nidx,Nidx)) Do i=1,Nidx X1d(1:Nt)=ClimIDX(i,1:Nt) Do j=1,Nidx Y1d(1:Nt)=ClimIDX(j,1:Nt) Call Correlation (X1d,Y1d,Nt,ctmp,ptmp, Miss) RClim(i,j)=ctmp; PClim(i,j)=ptmp Enddo Enddo Deallocate (X1d,Y1d) Open (3,file=Fname,Status="Unknown") write (3,*) "Correlation matrix between Climate indices" write (3,"(20A15)") "Cor ", (Clim_nm(i),i=1,Nidx) Do i=1,Nidx write (3,"(A15,20F15.3)") Clim_nm(i), (RClim(i,j),j=1,Nidx) Enddo write (3,*) write (3,"(20A15)") "p-value ", (Clim_nm(i),i=1,Nidx) Do i=1,Nidx write (3,"(A15,20F15.3)") Clim_nm(i), (PClim(i,j),j=1,Nidx) Enddo Close (3) End Subroutine Corr_ClimIDX_ClimIDX !-------------------------------------------------------------------- Subroutine Areal_Mean_of_Drought_and_Trend (Fname) Character (Len=*) :: Fname Integer :: i,i1,j,it,ireg,N1,NumYear,nt,nregs,nscal Real :: tmp,tmp1, SenSlope, prob Real, Allocatable :: xreg(:),xscal(:),xmon(:),X3d(:,:,:) Call rd_ncdf_3d (Fname,xreg,xscal,xmon,nregs,nscal,nt,X3d,"spe12") Do ireg=1,Nreg+1 Drou_Mean(ireg,1:Nt)=X3d(ireg,1,1:Nt) Enddo Deallocate (xreg,xscal,xmon,X3d) 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, parameter :: Mov_ave=13 Integer :: i,j,NumYear,N1,Nt1,ilon,ilat, md2 Real :: SenSlope, prob, X1d(Nyear), Ptb, XX(Nt), YY1(Nt), YY(Nt-Mov_ave+1) md2=Mov_ave/2; Nt1=Nt-Mov_ave+1 ! MK_Sen(Y,N,NumYear,Miss,SenSlope,prob,N1) NumYear=20 ! Minimum sample print*,"Calculating Trend drought..." print*,"Nt=",Nt; pause 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 !if (SenSlope/=Miss) 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 ! XX(1:Nt)=T2mon(ilon,ilat,1:Nt) !if (.not.Ar1_Empty(XX,Nt)) then ! Call Moving_Ave (XX,Nt,Mov_ave,YY1) ! YY(1:Nt1)=YY1(md2+1:Nt-md2) ! Call MK_Sen(YY,Nt1,NumYear,Miss,SenSlope,prob,N1) ! SenT2m(ilon,ilat)=SenSlope; ProT2m(ilon,ilat)=prob !endif 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 ! XX(1:Nt)=Precmon(ilon,ilat,1:Nt) !if (.not.Ar1_Empty(XX,Nt)) then ! Call Moving_Ave (XX,Nt,Mov_ave,YY1) ! YY(1:Nt1)=YY1(md2+1:Nt-md2) ! Ptb=Ave(YY,Nt1,Miss) ! Call MK_Sen(YY,Nt1,NumYear,Miss,SenSlope,prob,N1) ! SenPrec(ilon,ilat)=SenSlope/Ptb*100; ProPrec(ilon,ilat)=prob !endif endif Enddo Enddo End Subroutine Trend_drought_and_T_R !--------------------------------------------------------------------------------- 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 !--------------------------------------------------------------------------------- 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 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) 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, Nt1, it, nlo, nla, imon, imon1 ! Read in subregions: 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: Nt1=Nt-Nmon ! Bo nam dau tien Allocate (T2mon(Nlo,Nla,Nt1),Precmon(Nlo,Nla,Nt1)) T2mon(:,:,1:Nt1) =XX(:,:,Nmon+1:Nt) Precmon(:,:,1:Nt1)=YY(:,:,Nmon+1:Nt) ! De tinh trend Nyear = Nt/Nmon ! So nam Allocate (T2m(Nlo,Nla,Nyear),Prec(Nlo,Nla,Nyear)) Allocate (X1d(Nmon),Y1d(Nmon)) print*,"nlo,nla,nt,Nyear,Nt1=",nlo,nla,nt,Nyear,Nt1 Where (IsNaN(XX)); XX=Miss; Endwhere Where (IsNaN(YY)); YY=Miss; Endwhere ! Annual T2m and Prec: 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)= T2mon(ilon,ilat,imon:imon1) !XX(ilon,ilat,imon:imon1) Y1d(1:Nmon)=Precmon(ilon,ilat,imon:imon1) !YY(ilon,ilat,imon:imon1) T2m(ilon,ilat,it)=Ave(X1d,Nmon,Miss) Prec(ilon,ilat,it)=Ave(Y1d,Nmon,Miss) ! TB thang, *Nmon neu Tong nam. 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,"spei12") 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*,"Landpoint, ipoint: ",Landpoint, ipoint; stop !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),ClimIDX(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) ClimIDX(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*,trim(Clim_nm(iidx)) ! print "(12F7.3)",Y(iidx,:) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !------------------------------------------------------------------------------ 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"