! Chuong trinh tinh EOF cua 14 ClimIDX cho cac Lag tu 3-12 thang Use NetCDF ! Su dung SL cua cac thang tu 1979 den dau 2010 de Du bao cho cac thang tu 1/1980 den 12/2010 (Training) ! Real, Parameter :: Miss=-99.0 !Integer, Parameter :: Yr11=1979, Yr1=1980, Yr2=2010, Nyear=Yr2-Yr1+1, Nmon=12, Ntime=Nyear*Nmon, Ntime_BL=Ntime+Nmon !Integer, Parameter :: Yf11=2009, Yf1=2010, Yf2=2020, Nyfc= Yf2-Yf1+1, Ntime_fc=(Nyfc-1)*Nmon, Nlead=6 Integer, Parameter :: Yr1=1980, Yr2=2020, Nmon=12, Nyear=Yr2-Yr1, Ntime=Nyear*Nmon, Ntime0=Ntime+Nmon Integer, Parameter :: Yr1_BL=1981, Yr2_BL=2010, Nyr_BL=Yr2_BL-Yr1_BL+1 ! Thoi ky baseline de tinh trung binh khi hau !Integer, Parameter :: Nt0=(Yf1-Yr1)*Nmon+1 Integer, Parameter :: Nlag=6, Nlag0=2 ! Tre 6 thang, bat dau tu thang thu Nlag0 ke tu thang Du bao (VD: Leadtime=0 tai thang 4, thi Nlag0 tu thang 2) Integer, Parameter :: Nlen=9 ! R, T2m, R_ano, T_ano, R_std, T_std, SPEI_01, 03, 06 Character (Len=7), Parameter :: Dr_nm(Nlen)=(/"R","T2m","R_ano","T_ano","R_std","T_std","spei_01","spei_03","spei_06"/) Character (Len=50), Parameter :: Dr_ln(Nlen)=(/"Monthly rainfall (1/1981-12/2020)","Monthly 2m-temperature (1/1981-12/2020)", & "Monthly rainfall anomaly (1/1981-12/2020)","Monthly 2m-temperature anomaly (1/1981-12/2020)", & "Monthly rainfall standardization (1/1981-12/2020)", "Monthly 2m-temperature standardization (1/1981-12/2020)", & "The 1-month drought index - spei_01","The 3-month drought index - spei_03","The 6-month drought index - spei_06"/) Character (Len=7), Parameter :: Dr_unit(Nlen)=(/"mm/mon","C deg.","mm/mon","C deg.","-","-","-","-","-"/) ! Nidx=24: !Integer, Parameter :: Nidx=24 !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "oni","tni",& ! "mei", "censo", "ammsst", "NTA_ersst", & ! "soi", "reqsoi", "whwp", "amonus", "repac_slpa", & ! "rindo_slpa", "tsa", "tna","zwnd20", "wpac85", "epac85", & ! "cpac85", "qbo30", "wp"/) ! Cac chi so co do tre thoi gian 1 thang: ! Nino12, nino3, nino34, nino4, pna, best (censo), ammsst, nta_sst, soi, whwp, amonus, tna, qbo, wp, ao, nao, ReqSOI, ! Repac_slpa, Rindo_slpa, mei_v2, zwnd200, wpac850 !Integer, Parameter :: Nidx=22 !Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino12", "nino3", "nino34", "nino4", "pna", & !5 ! "censo", "ammsst", "NTA_ersst", & !3 ! "soi", "whwp", "amonus", "tna", "qbo30", "wp", "ao", & !8 ! "nao", "reqsoi", "repac_slpa", & !3 ! "rindo_slpa", "mei", "zwnd20", "wpac85"/) ! , "epac85", & !4 Integer, Parameter :: Nidx=14 Character (Len=15), Parameter :: Clim_nm(Nidx)=(/"nino3", "nino34", "nino4", "pna", & !5 "censo", & "soi", "qbo30", "wp", & "reqsoi", "repac_slpa", & !3 "rindo_slpa", "mei", "epac85", & !4 "cpac85"/) !, "qbo50"/) !2 Real, Allocatable :: Clim_Dat(:,:), Y_Clim(:,:) Real, Allocatable :: Clim_Fcs(:,:), Y_Clim_Fcs(:,:) Character (Len=200) :: FnameTC, FnameIDX, Fname, Out_Dir !-------------------------------------------------------------------------------------- ! Tinh cac EOF cua cac ClimIDX va tuong quan giua cac EOF voi TC_Char Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Real*8 :: Factor Real, Allocatable :: Xtb_Clim(:), Sx_Clim(:), PCs_Clim(:,:), Eig_Clim(:) Real, Allocatable :: lamda(:,:), PC1(:,:) Real, Allocatable :: EOF_Clim(:,:), ClimIDX(:,:) ! (Nidx,Ntime+Nmon) Loai bo CPAC850 va WPC850 vi SL ngan Real, Allocatable :: EOF_Fcs(:,:), EOF1(:,:), ClimFcs(:,:) Integer :: imon, iyr, iidx, imod, Nt_Fcs Real, Allocatable :: Corr_all_lag(:,:,:),Prob_all_lag(:,:,:) ! Nmon,Nmode,NTC_Cha !-------------------------------------------------------------------------------------- ! So bien duoc chon boi Stepwise Regression Integer, Allocatable :: Mvar (:,:) ! (NTC_Cha,Nmon) Integer, Allocatable :: IDX_Var(:,:,:) ! (NTC_Cha,Nx,Nmon) Integer, Parameter :: Mvar_X=9 Real, Allocatable :: X_Clim(:,:,:,:), X_Fcs(:,:,:,:) Integer :: NMax, MMax, N, M Real :: PCT, ANS(11) Integer :: NSTEP(5), i,j, NN, KSTEP Integer, Allocatable :: IDX(:),L(:) Real, Allocatable :: X(:),XBAR(:),STD(:),D(:),B(:),RX(:),R(:),T(:) Integer, Allocatable :: NSTEPS(:,:),LS(:,:) Real, Allocatable :: BS(:,:), ANSS(:,:) Real, Allocatable :: XMEAN(:),DX(:), RR(:,:), YORG(:), YHQ(:) !-------------------------------------------------------------------------------------- Real, Allocatable :: DrIDX(:,:,:), Dr_BL(:,:,:,:), Dr_Fcs(:,:,:,:) Real, Allocatable :: Mean(:,:,:), Std0(:,:,:) ! (Nst,Nmon,2) Real, Allocatable :: Lon(:), Lat(:) Character (Len=20),Allocatable :: Stname(:) Integer :: Nst, ist, Ntime_dr !-------------------------------------------------------------------------------------- Real, Allocatable :: Beta(:,:,:,:) ! Cac he so HQ Real, Allocatable :: Yout(:,:,:,:) ! Y-estimated Real, Allocatable :: Yfcs(:,:,:,:) ! Y-estimated for Independent forecast Real, Allocatable, Dimension(:,:,:) :: ME0, MAE0, RMSE0, RME0, Corr0, Pvalue0 ! Training Real, Allocatable, Dimension(:,:,:) :: ME1, MAE1, RMSE1, RME1, Corr1, Pvalue1 ! Forecast !-------------------------------------------------------------------------------------- Out_Dir="/home/tanpv/Progs/2022/Dr_VNU/RES/" !FnameTC="/home/tanpv/Progs/2022/Dr_VNU/RES/SPEI_Stn_VN_1980-2020.nc" FnameTC="/home/tanpv/OBS_Stn/Monthly/NC/R_T2m_1980-2020_mon_3D_VN.nc" Fname="/home/tanpv/Progs/2022/Dr_VNU/RES/SPEI_Stn_VN_1980-2020.nc" FnameIDX="/work/projects/G09/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" print*, "Calling Read_ClimIDX" Call Read_ClimIDX (FnameIDX) Call Read_Predictand (FnameTC) Call Read_Dr_IDX (Fname) Call Multi_lag Fname=trim(Out_Dir)//"Predictors_Clim.nc" Call Write_Predictors_Clim (Fname) Fname=trim(Out_Dir)//"Predictands.nc" FnameTC=trim(Out_Dir)//"R_T_Climatology.nc" Call Write_Predictands (Fname,FnameTC) stop !print*, "Calling Multi_lag" !Call Multi_lag ! Return: Clim_Dat switch=0; Standard=1; PC_Stard=0; Nmode=6 Nt=Ntime; Nx=Nidx*Nlag !-------------------------------------------------------------------------------------- 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)) Allocate (EOF_Clim(Nt,Nmode)) Allocate (Xtb_Clim(Nx), Sx_Clim(Nx), PCs_Clim(Nx,Nmode), Eig_Clim(Nmode)) !-------------------------------------------------------------------------------------- ! PCA: Do iyr=1,Nt !=Ntime Do iidx=1,Nx ! =Nidx*Nlag xdat(iyr,iidx)=Clim_Dat(iidx,iyr) Enddo Enddo Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) print"(A,6(I2,F8.2),F10.2)","Eigen - mode ",(imod, Eigenvalues(imod)/sum(Eigenvalues)*100,imod=1,Nmode),sum(Eigenvalues(1:Nmode))/sum(Eigenvalues)*100. !-------------------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------------- Subroutine Write_Predictands (Fname,Fname1) Character (Len=*) :: Fname,Fname1 Real, Allocatable :: Stn(:),xtime(:), X2d(:,:) Integer :: ist, it, iidx Allocate (Stn(Nst), xtime(Ntime_dr), X2d(Nst,Ntime_dr)) Do ist=1,Nst; Stn(ist)=ist; Enddo Do it=1,Ntime_dr; xtime(it)=it; Enddo Do iidx=1,Nlen Do ist=1,Nst Do it=1,Ntime_dr X2d(ist,it)=DrIDX(ist,it,iidx) Enddo Enddo if (iidx==1) then Call wr_ncdf_2d (Fname,Stn,xtime,Nst,Ntime_dr,X2d,trim(Dr_nm(iidx)),trim(Dr_ln(iidx)),trim(Dr_unit(iidx)),"Stn_id","month",Miss) else Call wr_ncdf_2d_apd (Fname,X2d,trim(Dr_nm(iidx)),Nst,Ntime_dr,trim(Dr_ln(iidx)),trim(Dr_unit(iidx)),"Stn_id","month",Miss) endif Enddo Call wr_ncdf_1d_apd (Fname,Lon,"lon",Nst,"Longtitude of Stations","Deg. East","Stn_id",Miss) Call wr_ncdf_1d_apd (Fname,Lat,"lat",Nst,"Latitude of Stations","Deg. North","Stn_id",Miss) Call wr_ncdf_1d_apd_char (Fname,Stname,"stname",Nst,"Name of stations","Character","Stn_id",Miss) ! Write Mean and Std Climatology (1981-2010): Mean(:,:,:), Std0(:,:,:) ! (Nst,Nmon,2) Deallocate (X2d,xtime) Allocate (xtime(Nmon),X2d(Nst,Nmon)) Do it=1,Nmon; xtime(it)=it; Enddo X2d(1:Nst,1:Nmon)=Mean(1:Nst,1:Nmon,1) Call wr_ncdf_2d (Fname1,Stn,xtime,Nst,Nmon,X2d,"r_mean","R mean 1981-2010","mm/mon","Stn_id","month",Miss) X2d(1:Nst,1:Nmon)=Mean(1:Nst,1:Nmon,2) Call wr_ncdf_2d_apd (Fname1,X2d,"t_mean",Nst,Nmon,"T mean 1981-2010","C deg.","Stn_id","month",Miss) X2d(1:Nst,1:Nmon)=Std0(1:Nst,1:Nmon,1) Call wr_ncdf_2d_apd (Fname1,X2d,"r_std",Nst,Nmon,"R std 1981-2010","mm/mon","Stn_id","month",Miss) X2d(1:Nst,1:Nmon)=Std0(1:Nst,1:Nmon,2) Call wr_ncdf_2d_apd (Fname1,X2d,"t_std",Nst,Nmon,"T std 1981-2010","C deg.","Stn_id","month",Miss) End Subroutine Write_Predictands !-------------------------------------------------------------------------------------- Subroutine Write_Predictors_Clim (Fname) Character (Len=*) :: Fname Real, Allocatable :: Clim_ID(:),xtime(:) ! Clim_Dat(Nidx*Nlag,Ntime) Integer :: ist, it, iidx, Nvar Nvar=Nidx*Nlag Allocate (Clim_ID(Nvar),xtime(Ntime)) Do iidx=1,Nvar; Clim_ID(iidx)=iidx; Enddo Do it=1,Ntime; xtime(it)=it; Enddo Call wr_ncdf_2d (Fname,Clim_ID,xtime,Nvar,Ntime,Clim_Dat,"clim","Climate indices","-","idx","month",Miss) End Subroutine Write_Predictors_Clim !-------------------------------------------------------------------------------------- Subroutine Read_Dr_IDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Stn(:),xmon(:), tmp(:,:) Integer :: Nstn, Nmn, icha, iyr, i,j, it1, it2 ! DrIDX, Dr_BL(:,:,:), Dr_Fcs(:,:,:) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_01") Nst = Nstn ! Num of Stn it1=Nmon+1; it2=Nmn ! Bo nam dau tien Ntime_dr = it2-it1+1 !Allocate (DrIDX(Nst,Ntime_dr,Nlen)) print*,"Nst,Ntime_dr,Nlen=",Nst,Ntime_dr,Nlen;pause DrIDX(:,1:Ntime_dr,7) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_03") DrIDX(:,1:Ntime_dr,8) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_06") DrIDX(:,1:Ntime_dr,9) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) !Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"spei_12") !DrIDX(:,1:Ntime_dr,4) = tmp(:,it1:it2); Deallocate (Stn,xmon,tmp) !Call rd_ncdf_1d (Fname,Stn,Nstn,Lon,"lon"); Deallocate (Stn) !Call rd_ncdf_1d (Fname,Stn,Nstn,Lat,"lat"); Deallocate (Stn) !Call rd_ncdf_1d_Char (Fname,Stn,Nstn,Stname,"stname"); Deallocate (Stn) End Subroutine Read_Dr_IDX !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) Character (Len=*) :: Fname Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Nyr, Nmn, iidx, iyr, i,j, k, imon !Yr11=Yr1-1 ! Lui lai 1 nam Allocate (ClimIDX(Nidx,Ntime0)) !, ClimFcs(Nidx,Ntime_fc+Nmon)) Do iidx=1,Nidx Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) print*,"Reading ",iidx," ",trim(Clim_nm(iidx)) j=0 Do iyr=Yr1,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do imon=1,Nmon j=j+1; ClimIDX(iidx,j)=tmp(i,imon) !print "(F5.0,12F8.3)",Year(i) !,ClimIDX(iidx,j) print*,"j=",j Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) !print*,Clim_nm(iidx);pause Enddo End Subroutine Read_ClimIDX !-------------------------------------------------------------------------------------- Subroutine Multi_lag Integer :: it, imon, iidx, ilag, imon1, it1, it2, k, Nlg0 Integer :: mn,yrs Real, Allocatable :: X(:), Y(:) Nlg0=Nmon-Nlag ! Lt=0 tai 01/1980, ClimIDX tinh tu Nmon-Lag0+1 = 11/1979 Allocate (Clim_Dat(Nidx*Nlag,Ntime)) !,Clim_Fcs(Nidx*Nlag,Ntime_fc)) k=0 Do iidx=1, Nidx ! imon=1, Nmon Do ilag=1,Nlag k=k+1; it1=0 !print*,"Nlg0+ilag-1,Nlg0+ilag-2+Ntime=",Nlg0+ilag-1,Nlg0+ilag-2+Ntime,(Nlg0+ilag-2+Ntime)-(Nlg0+ilag-1); pause Do it=Nlg0+ilag-1,Nlg0+ilag-2+Ntime it1=it1+1; Clim_Dat(k,it1)=ClimIDX(iidx,it) !print*,k,it1,it !mn=mod(it,Nmon); if (mn==0) mn=Nmon !yrs=it/Nmon; if (yrs*Nmon/=it) yrs=yrs+1; yrs=yrs+Yr11-1 !print"(A,6I6,F10.4)","it,yrs,mn,iidx,ilag,it1=",it,yrs,mn,iidx,ilag,it1,Clim_Dat(k,it1) Enddo !pause Enddo Enddo print*,"Nidx*Nlag,k=",Nidx*Nlag,k; pause End Subroutine Multi_lag !-------------------------------------------------------------------------------------- Subroutine Read_Predictand (Fname) Character (Len=*) :: Fname Real, Allocatable :: Stn(:),xmon(:), xyear(:), tmp1(:,:,:), tmp2(:,:,:), X1d(:),Y1d(:) Integer :: ist,imon, Nstn, Nmn, Nys, iyr, i,j, it1, it2 ! DrIDX, Dr_BL(:,:,:), Dr_Fcs(:,:,:) !Call rd_ncdf_2d (Fname,Stn,xmon,Nstn,Nmn,tmp,"R") Call rd_ncdf_3d (Fname,Stn,xyear,xmon,Nstn,Nys,Nmn,tmp1,"R"); Deallocate (Stn,xmon,xyear) Call rd_ncdf_3d (Fname,Stn,xyear,xmon,Nstn,Nys,Nmn,tmp2,"T2m"); ! Calculating Mean and Anomaly: !Allocate (X1d(Nys),Mean(Nstn,Nmon,Nlen)) ! (Nst,Nmon,Nlen)) !Do ist=1,Nstn ! Do imon=1,Nmon ! X1d(1:Nys)=tmp1(ist,1:Nys,imon); Mean(ist,imon,1)=Ave(X1d,Nys,Miss) ! X1d(1:Nys)=tmp2(ist,1:Nys,imon); Mean(ist,imon,2)=Ave(X1d,Nys,Miss) ! Enddo ! Do iyr=1,Nys ! Do imon=1,Nmon ! tmp1(ist,iyr,imon)=tmp1(ist,iyr,imon)-Mean(ist,imon,1) ! tmp2(ist,iyr,imon)=tmp2(ist,iyr,imon)-Mean(ist,imon,2) ! Enddo ! Enddo !Enddo Allocate (X1d(Nyr_BL),Y1d(Nyr_BL),Mean(Nstn,Nmon,2),Std0(Nstn,Nmon,2)) Do ist=1,Nstn Do imon=1,Nmon i=0 Do iyr=1,Nys if (Int(xyear(iyr))>=Yr1_BL.AND.Int(xyear(iyr))<=Yr2_BL) then i=i+1; X1d(i)=tmp1(ist,iyr,imon); Y1d(i)=tmp2(ist,iyr,imon) Mean(ist,imon,1)=Ave(X1d,Nyr_BL,Miss) ! R mean Mean(ist,imon,2)=Ave(Y1d,Nyr_BL,Miss) ! T2m mean Std0(ist,imon,1)=Stdev(X1d,Nyr_BL,Miss) Std0(ist,imon,2)=Stdev(Y1d,Nyr_BL,Miss) endif Enddo Enddo Enddo Nst = Nstn ! Num of Stn !it1=1; it2=Nmn !it1=Nmon+1; it2=Nmn ! Bo nam dau tien Ntime_dr = Nmn*(Nys-1) ! Bo nam dau tien (1980) Allocate (DrIDX(Nst,Ntime_dr,Nlen)) Do ist=1,Nst i=0 Do iyr=1,Nys-1 ! Tinh tu 1981 Do j=1,Nmon i=i+1; DrIDX(ist,i,1)=tmp1(ist,iyr,j); DrIDX(ist,i,2)=tmp2(ist,iyr,j) ! Origin values DrIDX(ist,i,3)=tmp1(ist,iyr,j)-Mean(ist,j,1) ! R Anomaly DrIDX(ist,i,4)=tmp2(ist,iyr,j)-Mean(ist,j,2) ! T Anomaly DrIDX(ist,i,5)=(tmp1(ist,iyr,j)-Mean(ist,j,1))/Std0(ist,j,1) ! R Standardisation DrIDX(ist,i,6)=(tmp2(ist,iyr,j)-Mean(ist,j,2))/Std0(ist,j,2) ! T Standardisation Enddo Enddo Enddo print*,"Nst,Ntime_dr,Nlen=",Nst,Ntime_dr,Nlen;pause Deallocate (Stn,xmon,xyear,tmp1,tmp2) Call rd_ncdf_1d (Fname,Stn,Nstn,Lon,"Lon"); Deallocate (Stn) Call rd_ncdf_1d (Fname,Stn,Nstn,Lat,"Lat"); Deallocate (Stn) Call rd_ncdf_1d_Char (Fname,Stn,Nstn,Stname,"Stname"); Deallocate (Stn) End Subroutine Read_Predictand !-------------------------------------------------------------------------------------- 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" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Stdev.f90" !Include "/home/tanpv/Mylibs/Sub_HQTT.f90" !Include "/home/tanpv/Mylibs/Sub_HQTB.f90" !-------------------------------------------------------------------------------------- END !-------------------------------------------------------------------------------------- !Include "/home/tanpv/Mylibs/Sub_HQTB.f90" Include "/home/tanpv/Mylibs/Sub_HQTT.f90" Include "/home/tanpv/Mylibs/Student_Prob.f90" !--------------------------------------------------------------------------------------