! 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=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(:,:) 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_cli(:,:) ! (Nlon_sst, Nlat_sst): climatological SST 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, iidx Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"pdsi","spi12","spe12"/) Character (Len=255) :: Filename, Dir_Out !--------------------------------------------------------------------------------- !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"/) !--------------------------------------------------------------------------------- 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 Dir_Out="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/" !--------------------------------------------------------------------------------- ! 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 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) Call Write_EOF_PCA (pca_out,eof_out,iidx) Filename="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/pca_APH_"// & trim(IDX_nm(iidx))//".txt" Call Write_PCA_to_TXT (Filename,iidx) ! For clustering only: Filename="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/pca_APH_clu.nc" Call Write_PCA_for_Clustering (Filename,iidx) print*, "===========> Complete for ",IDX_nm(iidx) Enddo Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) !--------------------------------------------------------------------------------- ! PCA for rainfall: Nt = Ntime_prec; Nx = Nst 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) Call Write_EOF_PCA (pca_prec,eof_prec,0) ! 0 mean prec Deallocate (xdat,ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate) !--------------------------------------------------------------------------------- ! 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)) 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 print*, "===========> Complete for SST" ! Composite SST maps: !--------------------------------------------------------------------------------- Call SST_Composite_maps ( & "/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/SST_composite_maps.nc") CONTAINS !================================================================================= Subroutine Write_PCA_to_TXT (Fname,iidx) Character (Len=*) :: Fname Integer :: i,j,k,imod,iidx Open (3,file=Fname,Status="Unknown") !------------------------------------------------------------------------- ! Write to PCs file: IDX_nm k=0 Do j=1, Nlat Do i=1, Nlon if (Mask_all(i,j)>0) then k=k+1; write (3,"(20F10.3)") xlon(i),xlat(j),(PCs(k,imod),imod=1,Nmode) endif Enddo Enddo Close(3) End Subroutine Write_PCA_to_TXT !--------------------------------------------------------------------------------- Subroutine Write_PCA_for_Clustering (F_pca,iidx) !Character (Len=*) :: Fname !Real, Allocatable :: PCA(:,:),xmode(:),xpoint(:) ! (Nx,Nmode) !Integer :: i,j,k,iidx,it Character (Len=*) :: F_pca Character (Len=255) :: var_ln, t_unit, IDX_name, Vname Integer :: i,j,k,iidx,it Real :: tmp, lonmax, latmax Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:) Real, Allocatable :: lon_std(:), lat_std(:) ! Standardize lon and lat: lonmax=xlon(Nlon); latmax=xlat(Nlat) Allocate (lon_std(Nlon),lat_std(Nlat)) Do i=1,Nlon; lon_std(i)=xlon(i)/lonmax; enddo Do j=1,Nlat; lat_std(j)=xlat(j)/latmax; enddo 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+1), xmode(Nmode+1), Xtmp(Nmode+1)) !, 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) PCA(i,j,Nmode+1) = Sqrt(lon_std(i)*lon_std(i)+lat_std(j)*lat_std(j)) endif Enddo Enddo Do i=1, Nmode+1; 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+1, 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+1, trim(var_ln), "-", "lon", "lat", "mode", Miss) endif !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 !!!" 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 (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, X3d, "pr") ! Precipitation Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "ori_pdsi") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "spei_12") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "spi_12") Deallocate (xlon, xlat) !Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land") !Call rd_ncdf_2d (trim(MaskFname), xlon, xlat, nlon,nlat,Mask,"land") Where (IsNaN(PDSI)); PDSI=Miss; Endwhere Where (IsNaN(SPEI12)); SPEI12=Miss; Endwhere Where (IsNaN(SPI12)); SPI12=Miss; Endwhere Where (IsNaN(X3d)); X3d=Miss; Endwhere Ntime = Nt - Nmon Allocate (Mask_all(nlon,nlat)) Mask_all=Mask Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime)) Allocate (X3(nlon, nlat, Ntime),Pr(nlon, nlat, Ntime)) !X1(:,1:nlat,1:Ntime)=PDSI(:,Nlat0+1:Nlat1,Nmon+1:Nt) !X2(:,1:nlat,1:Ntime)=SPEI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) !X3(:,1:nlat,1:Ntime)=SPI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) X1(:,1:nlat,1:Ntime) =PDSI (:,1:Nlat,Nmon+1:Nt) X3(:,1:nlat,1:Ntime) =SPEI12(:,1:Nlat,Nmon+1:Nt) X2(:,1:nlat,1:Ntime) =SPI12 (:,1:Nlat,Nmon+1:Nt) Pr(:,1:nlat,1:Ntime) =X3d (:,1:Nlat,Nmon+1:Nt) Deallocate (PDSI,SPEI12,SPI12,X3d) Allocate (PDSI(nlon,nlat,Ntime),SPEI12(nlon,nlat,Ntime)) Allocate (SPI12(nlon,nlat,Ntime)) !,X3d(nlon,nlat,Ntime)) PDSI=X1; SPEI12=X3; SPI12=X2; ! X3d=Pr Deallocate (X1, X2, X3) ist=0;ist1=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(PDSI (ilon,ilat,1:Ntime), Ntime).AND. & ! .not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. & ! .not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then !if (.not.Empty_Ar1(PDSI (ilon,ilat,1:Ntime), Ntime).AND. & ! .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 ist1=ist1+1 endif endif Enddo Enddo Nst=ist print*,"Nst= ",Nst," Nst1=",ist1," Sum=",Nlat*Nlon; pause Allocate (Dr_IDX(Ntime,Nst,NDr_idx),Prec(Ntime,Nst)) Dr_IDX(:,:,:) = Miss ist=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask_all(ilon,ilat)>0) then !if (.not.Empty_Ar1(PDSI (ilon,ilat,1:Ntime), Ntime).AND. & ! .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,1)=PDSI (ilon,ilat,1:Ntime) Dr_IDX(1:Ntime,ist,2)=SPI12 (ilon,ilat,1:Ntime) Dr_IDX(1:Ntime,ist,3)=SPEI12(ilon,ilat,1:Ntime) Prec (1:Ntime,ist) =Pr (ilon,ilat,1:Ntime) endif endif Enddo Enddo Deallocate (PDSI,SPEI12,SPI12,Pr) ! 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 = Nt Deallocate (X1d,Y1d,X2d) print*,"Reading data ... Done, Ntime_prec= ",Ntime_prec End Subroutine Read_Drough_IDX !-------------------------------------------------------------------------------------- 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 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.1 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, tmp0, tmp1, tmp2, ntmp0, 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/Post_APH/NC1/" & //"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 (sst_cli(Nlon_sst,Nlon_sst)) Allocate (Pval(Nlon_sst,Nlon_sst)) sst_dry=Miss; sst_wet=Miss; Pval=Miss sst_cli=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 tmp0=0; ntmp0=0 ! All cases Do it=1,Nt ntmp0=ntmp0+1; tmp0=tmp0+sst_3d(ilon,ilat,it) 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 tmp0=tmp0/ntmp0 ! Climatology 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 sst_cli(ilon,ilat)=tmp0 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) Call wr_ncdf_2d_apd (Fout,sst_cli, "sst_cli", Nlon_sst,Nlat_sst, & "Climatological Mean SST", "C", "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 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(:,:,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 ! 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"