! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Implicit None Real, Parameter :: Miss=-99.0, Lat_min=-11.0 ! Cat phan Autralia Integer, Parameter :: Nidx=16 ! 24 Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) !--------------------------------------------------------------------------------- Integer :: Yr1, Yr2, Nmon, Nyear ! Thoi ky Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Integer :: switch1, standard1, PC_Stard1, Nmode1 ! Clim Integer :: switch2, standard2, PC_Stard2, Nmode2 ! sst Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out Character (Len=250) :: Fnclim, pca_clim, eof_clim Character (Len=250) :: FnSST, pca_sst, eof_sst, pca_prec, eof_prec Real :: Slon1, Slon2, Slat1, Slat2, Miss_sst ! sst domain Integer :: Yr_sta Real*8 :: Factor Namelist /Model_Info/ Fname, MaskFname, pca_out, eof_out, & pca_prec, eof_prec, Yr1, Yr2, Nmon, & switch, standard, PC_Stard, Nmode Namelist /Clim_info/ Fnclim, pca_clim, eof_clim, & switch1, standard1, PC_Stard1, Nmode1 Namelist /SST_info/ FnSST, Yr_sta, Miss_sst, Slon1, Slon2, Slat1, & Slat2, pca_sst, eof_sst, & switch2, standard2, PC_Stard2, Nmode2 Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime, ilon, ilat Logical :: Status Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:) Real, allocatable :: ClimIDX(:,:) Real, allocatable :: Mask(:,:), Mask_all(:,:), Mask_dr(:,:) Real, allocatable :: Dr_IDX(:,:,:) !--------------------------------------------------------------------------------- ! SST: Real, allocatable :: lon_sst(:), lat_sst(:) Real, allocatable :: sst(:,:), Mask_sst(:,:) ! (Nt,Nx) Integer :: Nlon_sst, Nlat_sst, SeaPoint, Mov_Ave_Step !--------------------------------------------------------------------------------- ! Composite maps of SST: Real, allocatable :: sst_dry(:,:), sst_wet(:,:) ! (Nlon_sst, Nlat_sst): SST with Areal mean of SPEI<-0.5 and >-0.5 Real, allocatable :: sst_3d(:,:,:) ! (Nlon_sst, Nlat_sst,Ntime) Real, allocatable :: Areal_mean_Dr(:) ! (Ntime) !--------------------------------------------------------------------------------- ! Rainfall: Real, allocatable :: Prec(:,:), lon_prec(:), lat_prec(:), Mask_prec(:,:) Integer :: Ntime_prec, Nlon_prec, Nlat_prec, PrecPoint Integer :: Nst, Nst1, iidx Integer, Parameter :: NDr_idx=1 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"spe12"/) !,"spi12"/) Character (Len=255) :: Filename !--------------------------------------------------------------------------------- !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & ! "nino12", "nino3", "nino34", "nino4", "oni", "tni", & ! 6 !"whwp", "ammsst", "CAR_ersst", "NTA_ersst", "tsa", & ! 5 ! "whwp", "ammsst", "dmi", "NTA_ersst", "tsa", & ! 5 Bo CAR_ersst lay DMI ! "tna", "amonus", & ! "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 ! "rindo_slpa", "qbo30", "wp", "ao", "nao"/) ! 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"/) 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 !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_New.txt",Status="old") Read (1, Model_Info) Read (1, Clim_info) Read (1, SST_info) Close (1) Mov_Ave_Step = 13 ! Moving average !--------------------------------------------------------------------------------- ! Reading data: Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land") Do ilat=1,Nlat if (xlat(ilat) < Lat_min) Mask(:,ilat) = Miss Enddo print*,trim(Fname) Call Read_Drough_IDX (Fname) Nyear = Yr2-Yr1+1 ! So nam co SL Nt = Ntime; Nx = Nst !--------------------------------------------------------------------------------- if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) !--------------------------------------------------------------------------------- ! PCA for Dr_IDX: Loop over all drought indices: Do iidx=1,NDr_idx Xdat(1:Nt,1:Nx) = Dr_IDX(1:Nt,1:Nx,iidx) ! Tinh cac EOFs, PCs Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Deallocate (Mask_all) Allocate (Mask_all(Nlon,Nlat)) Mask_all(:,:)=Mask_dr(:,:) Call Write_EOF_PCA (pca_out,eof_out,iidx) ! For clustering only: !Filename=trim(IDX_nm(iidx))//"_"//trim(pca_out) !Call Write_PCA_for_Clustering (Filename) print*, "===========> Complete for ",IDX_nm(iidx) Enddo Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) !--------------------------------------------------------------------------------- ! PCA for rainfall: Nt = Ntime_prec; Nx = Nst1 Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) Factor = Sqrt(Real(Nt)) Xdat(1:Nt,1:Nx) = Prec(1:Nt,1:Nx) Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Deallocate (Mask_all) Allocate (Mask_all(Nlon,Nlat)) Mask_all(:,:)=Mask_prec(:,:) Call Write_EOF_PCA (pca_prec,eof_prec,0) ! 0 mean prec Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) print*, "===========> Complete for Rainfall..." !--------------------------------------------------------------------------------- ! PCA for Clim: !--------------------------------------------------------------------------------- Call Read_ClimIDX (Fnclim) Nt = Ntime; Nx = Nidx switch=switch1; standard=standard1; PC_Stard=PC_Stard1; Nmode=Nmode1 if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) Xdat(1:Nt,1:Nx) = ClimIDX(1:Nt,1:Nx) Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Call Write_EOF_PCA_Clim (pca_clim,eof_clim) Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) !print*,"Program is normal terminating ..." print*, "===========> Complete for ClimIDX" ! PCA for SST: !--------------------------------------------------------------------------------- Call Read_sst_data (FnSST) Nt = Ntime; Nx = SeaPoint switch=switch2; standard=standard2; PC_Stard=PC_Stard2; Nmode=Nmode2 if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) !print*,Nt,Nx,Mcount,Ncount; stop Xdat(1:Nt,1:Nx) = sst(1:Nt,1:Nx) Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) Deallocate (xlon, xlat, Mask_all) Nlon=Nlon_sst; Nlat=Nlat_sst Allocate (xlon(Nlon), xlat(Nlat), Mask_all(Nlon,Nlat)) xlon(:)=lon_sst(:); xlat(:)=lat_sst(:); Mask_all(:,:)=Mask_sst(:,:) Call Write_EOF_PCA (pca_sst,eof_sst,-1) ! -1 means sst Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) print*, "===========> Complete for SST" ! Composite SST maps: !--------------------------------------------------------------------------------- !Call SST_Composite_maps ("SST_composite_maps.nc") CONTAINS !================================================================================= Subroutine Write_PCA_for_Clustering (Fname) Character (Len=*) :: Fname Real, Allocatable :: PCA(:,:),xmode(:),xpoint(:) ! (Nx,Nmode) Integer :: i,j,k,iidx,it Allocate (PCA(Nx,Nmode),xmode(Nmode),xpoint(Nx)) Do i=1, Nmode; xmode(i)=i; Enddo Do i=1,Nx; xpoint(i)=i; Enddo Do i=1,Nmode; PCA(1:Nx,i)=PCs(1:Nx,i); Enddo !if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_2d (Fname,xpoint,xmode,Nx,Nmode,PCA,"pca", & "PC values of Drought index", "-", "point","mod", Miss) !else ! Call wr_ncdf_2d_apd (Fname,PCA,"???",Nx,Nmode, & ! trim(var_ln), "-", "mon", "mode", Miss) !endif End Subroutine Write_PCA_for_Clustering !------------------------------------------------------------------------- Subroutine Write_EOF_PCA_Clim (F_pca, F_eof) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln, t_unit Integer :: i,j,k,iidx,it Real :: tmp,sumx Real, Allocatable :: PCA(:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) Real, Allocatable :: Percent(:), Accum(:), xtime(:) Allocate (PCA(Nmode,Nx),xmode(Nmode), Xtmp(Nmode),EOF1(Nt,Nmode)) Allocate (Percent(Nmode),Accum(Nmode),X1(Nx),xtime(Nt)) !------------------------------------------------------------------------- Do i=1,Nmode; PCA(i,1:Nx)=PCs(1:Nx,i); Enddo Do i=1, Nmode; xmode(i)=i; Enddo Do i=1,Nx; X1(i)=i; Enddo ! Clim-idx Xtmp(1:Nmode) = Eigenvalues(1:Nmode) ! Eigen values ! Accum (%) of Eigen values values: Sumx=Sum(Eigenvalues(1:Mcount)); Accum(:)=0 Do i=1,Nmode Percent(i)=Xtmp(i)/Sumx*100 Accum(i)=Sum(Percent(1:i)) Enddo var_ln="PCA of ClimIDXs" Call wr_ncdf_2d (F_pca,xmode,X1,Nmode,Nx,PCA,"pca", & trim(var_ln), "-", "mod","idx", Miss) Call wr_ncdf_1d_apd (F_pca, Percent, & "per",Nmode, & "Percent Contribution of clim idx", & "-", "mod", Miss ) Call wr_ncdf_1d_apd (F_pca, Accum, & "cum",Nmode, & "Accumulated Contribution of clim idx", & "-", "mod", Miss ) Call wr_ncdf_1d_apd (F_pca, Xtmp, & "eig",Nmode, & "Eigen values of clim idx", & "-", "mod", Miss ) !------------------------------------------------------------------------- Do i=1, Nt; xtime(i)=i; EOF1(i,1:Nmode)=EOFs(i,1:Nmode); Enddo ! Months var_ln="Timeseries of eofs of ClimIDX" Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1,"eof", & trim(var_ln), "-", "mon","mod", Miss) End Subroutine Write_EOF_PCA_Clim !------------------------------------------------------------------------- Subroutine Write_EOF_PCA (F_pca, F_eof, iidx) Character (Len=*) :: F_eof, F_pca Character (Len=255) :: var_ln, t_unit, IDX_name, Vname Integer :: i,j,k,iidx,it Real :: tmp Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) if (iidx>=1.AND.iidx<=3) then IDX_name=IDX_nm(iidx); Vname=IDX_nm(iidx) else if (iidx==0) then IDX_name="Precipitation"; Vname="prec" else IDX_name="SST or other variable"; Vname="sst" endif Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode)) !, EOF1(Nmode,Nt)) Allocate (EOF1(Nt,Nmode)) !------------------------------------------------------------------------- ! Write to PCs file: PCA(:,:,:)=Miss; k=0 Do j=1, Nlat Do i=1, Nlon if (Mask_all(i,j)>0) then k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode) endif Enddo Enddo Do i=1, Nmode; xmode(i)=i; Enddo ! PCA's modes var_ln="Spatial distribution of PCAs of "//trim(IDX_name) write (t_unit,"(A,I4.4,A)") "months since ",Yr1,"-01-01 " if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_3d (F_pca,xlon,xlat,xmode,nlon,nlat,Nmode, PCA, & trim(Vname), trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) else Call wr_ncdf_3d_apd (F_pca,PCA,trim(Vname),nlon,nlat, & Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss) endif !------------------------------------------------------------------------- Xtmp(1:Nmode) = Eigenvalues(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(Vname)//"_eig", Nmode, & "Eigen values of "//trim(IDX_name), "-", "mode", Miss) tmp=0.0 Do i=1,Mcount tmp=tmp+Eigenvalues(i) if (I<=Nmode) then ECumulate(i)=0 Do j=1,I ECumulate(i)=ECumulate(i)+Eigenvalues(j) Enddo endif Enddo Do i=1,Nmode ECumulate(i)=ECumulate(i)/tmp; Xtmp(i)=Eigenvalues(i)/tmp Enddo Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(Vname)//"_var",Nmode, & "Percent Contribution of "//trim(IDX_name), & "-", "mode", Miss ) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_pca, Xtmp, & trim(Vname)//"_cum",Nmode, & "Cumuluted Eigenvalues of "//trim(IDX_name), & "-", "mode", Miss ) !------------------------------------------------------------------------- ! Write to EOF file: if (Allocated(xtime)) Deallocate (xtime); Allocate (xtime(1:Nt)) Do it=1, Nt; xtime(it)=it; Enddo ! Months Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(1:Nt,i); Enddo var_ln="Timeseries of eofs of "//trim(IDX_name) if (iidx==1.OR.iidx<1.OR.iidx>3) then Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, & trim(Vname), & trim(var_ln), "-", "mon","mode", Miss) else Call wr_ncdf_2d_apd (F_eof,EOF1, trim(Vname), Nt,Nmode, & trim(var_ln), "-", "mon", "mode", Miss) endif print*,"Write_EOF_PCA: Done for case: ",iidx End Subroutine Write_EOF_PCA !--------------------------------------------------------------------------------- Subroutine Read_Drough_IDX (Fname) Integer :: M,ist,ist1, ilon, ilat, it, Nt, imon, imon1, iyr Integer :: nlat1, nlat0 Character (Len=*) :: Fname Real :: tmp, xmean(12), ic(12) Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) Real, Allocatable :: Pr(:,:,:), X3d(:,:,:), X2d(:,:), X1d(:), Y1d(:) Real, Allocatable :: xlon(:), xlat(:), xtime(:) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,X3d,"Pr") ! Precipitation Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,SPEI12, & "spei12") Deallocate (xlon, xlat, xtime) Where (IsNaN(SPEI12)) SPEI12=Miss Where (IsNaN(X3d)) X3d=Miss Ntime = Nt - Nmon Allocate (Mask_all(nlon,nlat),Mask_dr(nlon,nlat),Mask_prec(nlon,nlat)) Mask_all=Mask Mask_dr=Miss; Mask_prec=Miss print*,"nlon, nlat, nt, Ntimei, Nmon: ",nlon, nlat, nt, Ntime, Nmon Allocate (X1(nlon, nlat, Ntime)) Allocate (Pr(nlon, nlat, Ntime)) !print*,Mask Do ilon=1,Nlon !print*,ilon Do ilat=1,Nlat Do it=1,Ntime X1(ilon,ilat,it) =SPEI12(ilon,ilat,Nmon+it) Pr(ilon,ilat,it) =X3d (ilon,ilat,Nmon+it) Enddo Enddo Enddo print*,"nlon, nlat, nt, Ntime: ",nlon, nlat, nt, Ntime Deallocate (SPEI12,X3d) !SPI12,X3d) Allocate (SPEI12(nlon,nlat,Ntime)) SPEI12=X1 Deallocate (X1) ist=0 ! To ensure all gridpoints do not contain missing values Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat)>0) then !if (.not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. & if (.not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then ist=ist+1 else Mask_all(ilon,ilat)=Miss; Mask(ilon,ilat)=Miss endif endif !print*,ist Enddo Enddo Nst=ist print*,"Nst= ",Nst Allocate (Dr_IDX(Ntime,Nst,NDr_idx),Prec(Ntime,Nst)) Dr_IDX(:,:,:) = Miss ist=0; ist1=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask_all(ilon,ilat)>0) then !if (.not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. & if (.not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then ist=ist+1 !Dr_IDX(1:Ntime,ist,2)=SPI12 (ilon,ilat,1:Ntime) Dr_IDX(1:Ntime,ist,1)=SPEI12(ilon,ilat,1:Ntime) Mask_dr(ilon,ilat)=1 endif if (.not.Empty_Ar1(Pr(ilon,ilat,1:Ntime), Ntime)) then ist1=ist1+1 Prec(1:Ntime,ist1)=Pr(ilon,ilat,1:Ntime) Mask_prec(ilon,ilat)=1 endif endif Enddo Enddo Nst1=ist1 Deallocate (SPEI12,Pr) !SPI12,Pr) ! Prec anomaly OR moving average: Allocate (X1d(Ntime)) Do ist=1,Nst1 X1d(1:Ntime)=Prec(1:Ntime,ist) Call Anomaly (X1d,Ntime) Prec(1:Ntime,ist)=X1d(1:Ntime) Enddo ! Moving Average for Prec: !M=Mov_Ave_Step ! Buoc truot = 13 thang !Nt=Ntime-M+1 !--> Da bo sung !Allocate (X2d(Nt,Nst), Y1d(Nt), X1d(Ntime)) !Do ist=1,Nst ! X1d(1:Ntime)=Prec(1:Ntime,ist) ! Call Moving_Ave (X1d,Ntime,M,Y1d) ! X2d(1:Nt,ist)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime) !Enddo !Deallocate (Prec) !Allocate (Prec(Nt,Nst)) !Prec(1:Nt,1:Nst) = X2d(1:Nt,1:Nst) Ntime_prec = Ntime Deallocate (X1d) !,Y1d,X2d) print*,"Reading data ... Done, Ntime_prec= ",Ntime_prec End Subroutine Read_Drough_IDX !------------------------------------------------------------------------------ Subroutine Anomaly (X,Ntime) Integer :: Ntime Real :: X(Ntime), Y(Ntime), TB(Nmon) Integer :: ic, it, imon, imon1 Real :: tmp if (mod(Ntime,Nmon)/=0) print*,"Something wrong !" Do imon=1,Nmon tmp=0; ic=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then if (X(it)/=Miss) then ic=ic+1; tmp=tmp+X(it) endif endif Enddo if (ic>0) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif Enddo Do it=1,Ntime imon=mod(it,Nmon); if (imon==0) imon=Nmon Y(it)=X(it)-TB(imon) Enddo X(:) = Y(:) End Subroutine Anomaly !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k Yr11=Yr1 !+1 ! Bo hoac lay ca nam dau tien Ntime=(Yr2-Yr11+1)*Nmon Allocate (ClimIDX(Ntime,Nidx)) !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; ClimIDX(j,iidx)=tmp(i,k) Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) !print "(12F7.3)",Y(iidx,:) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------- Subroutine SST_Composite_maps (Fout) Character (Len=*) :: Fout Character (Len=255) :: Fname Real, Allocatable :: X2d(:,:), xreg(:), time(:) Real, Allocatable :: Pval(:,:) Real, parameter :: Dr_thres=-0.25 Integer :: Nregs, Nt Integer :: ist, i, j, ilon, ilat, ipoint, it, nlon, nlat Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2 Real :: T_Dist, tmp1, tmp2, ntmp1, ntmp2, Sx, Sy ! Real, allocatable :: sst_dry(:,:), sst_wet(:,:) ! (Nlon_sst, Nlat_sst): SST with Areal mean of SPEI<-0.5 and >-0.5 ! Real, allocatable :: Areal_mean_Dr(:) ! (Ntime) Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_Sta/" & //"Areal_mean_and_Trend_APH_spei_12.nc" Call rd_ncdf_2d (Fname,xreg,time,Nregs,Nt,X2d,"dridx") Allocate (Areal_mean_Dr(Nt)) Areal_mean_Dr(1:Nt) = X2d(Nregs,1:Nt) ! Chi lay trung binh toan SEA Deallocate (X2d,xreg,time) Allocate (sst_dry(Nlon_sst,Nlon_sst),sst_wet(Nlon_sst,Nlon_sst)) Allocate (Pval(Nlon_sst,Nlon_sst)) sst_dry=Miss; sst_wet=Miss; Pval=Miss Do ilon=1,Nlon_sst Do ilat=1,Nlat_sst if (Mask_sst(ilon,ilat)==1) then tmp1=0; tmp2=0; ntmp1=0; ntmp2=0; Sx=0; Sy=0 Do it=1,Nt if (Areal_mean_Dr(it)<=Dr_thres) then ntmp1=ntmp1+1 tmp1=tmp1+sst_3d(ilon,ilat,it) Sx=Sx+sst_3d(ilon,ilat,it)*sst_3d(ilon,ilat,it) else if (Areal_mean_Dr(it)>=ABS(Dr_thres)) then ntmp2=ntmp2+1 tmp2=tmp2+sst_3d(ilon,ilat,it) Sy=Sy+sst_3d(ilon,ilat,it)*sst_3d(ilon,ilat,it) endif Enddo tmp1=tmp1/ntmp1; tmp2=tmp2/ntmp2 Sx=Sx/(ntmp1-1); Sy=Sy/(ntmp2-1) Sx=SQRT(Sx-tmp1*tmp1); Sy=SQRT(Sy-tmp2*tmp2) sst_dry(ilon,ilat)=tmp1 sst_wet(ilon,ilat)=tmp2 Pval(ilon,ilat)=P_value(tmp1, tmp2, Sx, Sy, ntmp1, ntmp2) endif Enddo Enddo ! Lon_sst(Nlon_sst), Lat_sst(Nlat_sst) Call wr_ncdf_2d (Fout,Lon_sst,Lat_sst, Nlon_sst,Nlat_sst,sst_dry, & "sst_dry","Mean SST in case Dr_idx<-0.5", & "C", "lon","lat", Miss) Call wr_ncdf_2d_apd (Fout,sst_wet, "sst_wet", Nlon_sst,Nlat_sst, & "Mean SST in case Dr_idx>0.5", "C", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,Pval,"pval",Nlon_sst,Nlat_sst, & "P-value of testing if sst_dry /= sst_wet","-","lon", "lat", Miss) End Subroutine SST_Composite_maps !--------------------------------------------------------------------------------- Real Function P_value(Xtb, Ytb, Sx, Sy, Nx, Ny) ! Testing whether Xmean and Ymean are different Real :: t, Xtb, Ytb, Sx, Sy, Nx, Ny, df, Sx2, Sy2 Real :: T_Dist, A, tmp1, tmp2, tmp3, tmp4, pval df=nx+ny-2 if (Xtb/=Miss .AND. Ytb/=Miss .AND. Sx>0 .AND. Sy>0 .AND. Nx>1 .AND. Ny>1) then Sx2=Sx*Sx; Sy2=Sy*Sy ! S^2 tmp1=1.0/real(nx)+1.0/real(ny); tmp1=1.0/SQRT(tmp1) ! 1/Sqrt(1/n1 + 1/n2) tmp2=SQRT(df) tmp3=SQRT((nx-1)*Sx2 + (ny-1)*Sy2) tmp4=tmp3/tmp2 A=tmp1/tmp4 t=ABS(A*(Xtb-Ytb)) pval=T_Dist (t, df) else pval=Miss endif P_value = pval End Function P_value !-------------------------------------------------------------------- Subroutine Read_sst_data (Fname) Character (Len=*) :: Fname !Integer, Parameter :: Yr_sta=1891 ! Nam bat dau chuoi SL SST Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:) Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:), xtime(:) Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2 !Character (Len=255) :: Fn print*, "Reading SST ..." Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,XX,"sst") !trim(varn)) Where (IsNaN(XX)); XX=Miss; Endwhere Where (XX==Miss_sst); XX=Miss; Endwhere i=1; Do while (tlon(i)Slat1); j=j+1; Enddo; ilat1=j j=1; Do while (tlat(j)>Slat2); j=j+1; Enddo; ilat2=j-1 it1=(Yr1-Yr_sta)*Nmon+1 ! + Nmon ! Bo nam dau tien it2=it1+Ntime-1 ! Lay them M/2 phan tu dau va M/2 phan tu cuoi cho Moving Average: !M = Mov_Ave_Step !Md2=M/2 !it1=it1-Md2; it2=it2+Md2 Nt = Ntime ! Luu gia tri Ntime Ntime = it2-it1+1 ! Kich thuoc chuoi moi dung trong Moving Average Nlon_sst=ilon2-ilon1+1; Nlat_sst=ilat1-ilat2+1 ! Chu y: ilat1>ilat2 vi Latitude giam tu North to South Allocate (YY(Nlon_sst,Nlat_sst,Ntime),Lon_sst(Nlon_sst), Lat_sst(Nlat_sst)) Allocate (sst_3d(Nlon_sst,Nlat_sst,Nt)) ! Luu lai de su dung trong composite maps, Nt = do dai chuoi chua them Md2 Allocate (Mask_sst(Nlon_sst,Nlat_sst), X1d(Ntime)) Mask_sst(:,:)=Miss !X1d(1:Ntime) = xtime(it1:it2) i=0 Do ilon=1,Nlon if (ilon>=ilon1.AND.ilon<=ilon2) then i=i+1; Lon_sst(i)=tlon(ilon) j=Nlat_sst+1 Do ilat=1,Nlat if (ilat>=ilat2.AND.ilat<=ilat1) then j=j-1; Lat_sst(j)=tlat(ilat) YY(i,j,1:Ntime) = XX(ilon,ilat,it1:it2) endif Enddo endif Enddo sst_3d(:,:,1:Nt) = YY(:,:,1:Nt) ! YY(:,:,Md2+1:Ntime-Md2) ! Luu lai !print*,"check: Nt,Ntime=",Nt,Ntime,Md2+1,"-",Ntime-Md2; pause ! Processing for MCA: SeaPoint=0 Do i=1,Nlon_sst Do j=1,Nlat_sst X1d(1:Ntime) = YY(i,j,1:Ntime) if (.not.Empty_Ar1(X1d,Ntime)) then SeaPoint=SeaPoint+1; Mask_sst(i,j)=1 endif Enddo Enddo ! Anomaly: print*,"SST: ",Nt,SeaPoint Allocate (sst(Nt,SeaPoint), Y1d(Nt)) ipoint=0 Do ilat=1, Nlat_sst Do ilon=1, Nlon_sst if (Mask_sst(ilon,ilat)/=Miss) then X1d(1:Ntime)=YY(ilon,ilat,1:Ntime) Call Anomaly (X1d,Ntime) ipoint=ipoint+1 sst(1:Nt,ipoint)=X1d(1:Nt) endif Enddo Enddo print*,"SST: ",ipoint ! Moving Average: !M=Mov_Ave_Step ! Buoc truot = 13 thang !Nt=Ntime-M+1 !--> Da bo sung !Allocate (sst(Nt,SeaPoint), Y1d(Nt)) !ipoint=0 !Do ilat=1, Nlat_sst ! Do ilon=1, Nlon_sst ! if (Mask_sst(ilon,ilat)/=Miss) then ! ipoint=ipoint+1 ! X1d(1:Ntime)=YY(ilon,ilat,1:Ntime) ! Call Moving_Ave (X1d,Ntime,M,Y1d) ! sst(1:Nt,ipoint)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime) ! endif ! Enddo !Enddo !Ntime = Nt Deallocate (YY, XX, X1d,Y1d,tlat,tlon,xtime) End Subroutine Read_sst_data !--------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N-m+1), tmp Do j=1,N-m+1 tmp=0 Do i=j,j+m-1 tmp=tmp+X(i) Enddo Y(j)=tmp/Real(m) Enddo End Subroutine Moving_Ave !--------------------------------------------------------------------------------- Logical Function Empty_Ar1(X,N) Integer :: N,i Real :: X(N) Logical :: Empty Empty=.False. Do i=1,N if (X(i)==Miss.OR.IsNaN(X(i)).OR.ABS(X(i))>1000) then Empty=.True.; Exit endif Enddo Empty_Ar1=Empty End Function Empty_Ar1 !--------------------------------------------------------------------------------- !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" END !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/Student_Prob.f90"