Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Parameter :: Dr_Lev = 3 Real, parameter :: Dr_The(Dr_Lev)=(/-0.5, -1.0, -1.5/) Character (Len=4), Parameter :: Dr_Lnm(Dr_Lev)=(/"mild","mode","seve"/) Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6 Integer, Parameter :: Ncase = 5 ! "obs","mod","cli","ann","mlr" Integer, Parameter :: Ncase_1 = Ncase-1 ! "mod","cli","ann","mlr" Integer, Parameter :: Ncase_2 = 7 ! "ann_reg","mlr_reg","ann_clw","mlr_clw","ens","ens1","a_m" ("a_m" = f("ann","mlr") Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"mod","cli","ann","mlr"/) Character (Len=3), Parameter :: Case_2(Ncase_2)=(/"an1","mr1","an2","mr2","ens","en1","a_m"/) Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","mod","cli","ann","mlr"/) !Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","cli","ann"/) Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:), T2m(:), PET(:) ! For calculating SPI, SPEI ! Distribution function's parameters: Real, Allocatable :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:) Real :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2) Logical :: use_baseline Character (Len=5) :: IDX_nm ! Name of Drought Index Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Pr_obs(:,:,:), T2_obs(:,:,:) ! Nlon,Nlat,Ntime Real, Allocatable :: Pr_All(:,:,:,:,:,:,:), T2_All(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel Real, Allocatable :: xmem(:), xlead(:), xtime(:) Real, Allocatable :: SPEI_All(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel Real, Allocatable :: SPEI_M1 (:,:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_1,Nmodel) - Mem 1 only Real, Allocatable :: SPEI_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Ensemble mean over all Nmem*Ncase*Nmodel Real, Allocatable :: SPEI_ens1(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Ensemble Median over all Nmem*Ncase*Nmodel Real, Allocatable :: SPEI_obs(:,:,:) ! Nlon,Nlat,Ntime Real, Allocatable :: SPEI_Out(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel,Nlen Real, Allocatable :: SPEI_sel(:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_2 Real, Allocatable :: SPEI_sel_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Ensemble mean OF SPEI_sel !------------------------------------------------------------------------------ ! MLR: Real, Allocatable :: SPEI_mlr(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - MLR using SPEI_ens and SPEI_ens Real, Allocatable :: Xorg(:), Xhq(:) Real, Allocatable :: Beta(:,:,:,:) ! (Nlon,Nlat,Nlead,Mhq) Integer :: Mhq, Mvar, Nhq !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:,:,:) :: PCs, PODs, FARs, FBIs, TSs, HSSs ! Nlon,Nlat,Nlead,Nmon,Dr_Lev Real, Allocatable, Dimension(:,:,:,:) :: PCall, PODall, FARall, FBIall, TSall, HSSall ! Nlon,Nlat,Nlead,Dr_Lev ! For all timeseries Real, Allocatable, Dimension(:,:,:,:) :: Geo_F ! Nreg,Dr_Lev,Nlead,Ntime Real, Allocatable, Dimension(:,:,:) :: Geo_O ! Nreg,Dr_Lev,Ntime Real, Allocatable, Dimension(:,:,:,:,:) :: Corre, Probe ! Nlon,Nlat,Nmodel,Ncase_1,Nlead !------------------------------------------------------------------------------ Real, Allocatable :: SPEI_All_Reg(:,:,:,:,:,:) ! Nreg,Nmem,Nlead,Ntime,Ncase,Nmodel Real, Allocatable :: SPEI_ens_Reg(:,:,:) ! Nreg,Nlead,Ntime - Median of Ensemble mean Real, Allocatable :: SPEI_ens1_Reg(:,:,:) ! Nreg,Nlead,Ntime - Median of Ensemble Median Real, Allocatable :: SPEI_mlr_Reg(:,:,:) ! Nreg,Nlead,Ntime - MLR Real, Allocatable :: SPEI_obs_Reg(:,:) ! Nreg,Ntime - Median !------------------------------------------------------------------------------ ! Dr probability: Real, Allocatable :: Prob_mild(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Prob occurrence of mild drought Real, Allocatable :: Prob_mode(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Prob occurrence of moderate drought Real, Allocatable :: Prob_seve(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - Prob occurrence of seve drought !------------------------------------------------------------------------------ ! Ensemble: Real, Allocatable :: Pr_ens (:,:,:,:,:), T2_ens(:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_1 Real, Allocatable :: Pet_ens (:,:,:,:,:) !Real, Allocatable :: SPEI_ens(:,:,:,:,:) ! Nlon,Nlat,Ntime,Ncase_1,Nlen Real, Allocatable :: Pr_Reg (:,:,:,:), T2_Reg(:,:,:,:) ! (Nreg,Nlead,Ncase_1,Ntime) - Obs va Ens Real, Allocatable :: Pet_Reg(:,:,:,:) Real, Allocatable :: SPI_Reg(:,:,:,:,:) ! (Nreg,Nlead,Ncase_1,Ntime,Nlen) !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File Character (Len=3) :: TVar, RVar Character (Len=10),Allocatable :: Model_nm(:), IDXname(:) Logical :: status, Empty Integer :: Yr1, Yr2, Nmodel, NYear, maxyrs, Nidx Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat Integer :: Nt1 !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, xspi !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Nmodel, NLen, Nidx Namelist /Model_name/ Model_nm, IDXname Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_All_Dr_IDXs_1983-2020.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Len_SPI(NLen)); Allocate (Model_nm(Nmodel)); Allocate (IDXname(Nidx)) Read (1, Model_name) Read (1, Dr_Timescale) Close (1) ! Inp_Dir = "/home/tanpv/Data/562_VNU/0.25_VN/2024/Combined/" ! Out_Dir = "/home/tanpv/Data/562_VNU/0.25_VN/Dr_IDXs_Ope/" !------------------------------------------------------------------------ Ntime = (Yr2-Yr1+1)*Nmon; maxyrs = Yr2-Yr1+1; NYear = (Yr2 - Yr1 + 1) Call Read_Regs print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime Call Reinit ! Call Read_all_data ! Read all data of "mod","cli","ann","mlr" ! Call Test !stop ! Out_Dir = "/home/tanpv/Data/562_VNU/0.25_VN/Dr_IDXs_Ope/" Allocate (SPEI_Out(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel,Nlen)); SPEI_Out=Miss Do idx=1,Nidx ! "spi", "spei" Do ilen=1,Nlen ! 1, 3, 6 Call Read_all_Dr_IDXs (trim(IDXname(idx)),ilen) ! Call Calc_Correlation ! Call Write_Corre (trim(IDXname(idx)),ilen) Call Calc_IDX_Ens ! Ens over cli,ann and mlr. Outputs: SPEI_ens,SPEI_ens1 Call Calc_HQTT ! HQ over "ann" and "mlr". Outputs: SPEI_mlr. Based on SPEI_All Call Selected_IDXs ! Combine "ann" and "mlr" of eg and clw, and SPEI_ens,SPEI_ens1,SPEI_mlr (7) Call Write_Add_DrIDXs (trim(IDXname(idx)),ilen) Call Calc_Dr_Reg ! Based on outputs of Calc_IDXs Call Calc_Dr_events ! Based on outputs of Calc_IDXs Call Verification ! Based on outputs of Calc_Dr_events Call Write_Verification (trim(IDXname(idx)),ilen) Call Calc_Geo_Extent ! Based on outputs of Calc_Dr_events Call Write_GeoExt (trim(IDXname(idx)),ilen) Call Write_Dr_Reg_ens (trim(IDXname(idx)),ilen) Enddo Enddo !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_HQTT Integer :: i,j,k, ilead,imon,ireg,ilon,ilat ! SPEI_mlr(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime if (Allocated (SPEI_mlr)) Deallocate (SPEI_mlr) Allocate (SPEI_mlr(Nlon,Nlat,Nlead,Ntime)) Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Call HQTT (ilon,ilat,ilead) SPEI_mlr(ilon,ilat,ilead,1:Ntime)=Xhq(1:Ntime) !print*,ilead,ilon,ilat,Xhq(10) endif Enddo Enddo Enddo End Subroutine Calc_HQTT !------------------------------------------------------------------------------ Subroutine HQTT (ilon,ilat,ilead) ! SPEI_All(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase_1,Nmodel ! SPEI_mlr(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime - MLR using SPEI_ens and SPEI_ens Integer :: i,j,k,it, ilead,imon,ireg,ilon,ilat,imod,icas,imem Real :: S,Q,U,RR,F, tmp Real, Allocatable :: X2d(:,:), Check(:), A(:) Real :: X_inp(Ntime,Nmodel*Ncase_1-2),Y_inp(Ntime), xtmp(Nmodel*Ncase_1) Logical :: OK ! Allocate (SPEI_ens(Nlon,Nlat,Nlead,Ntime)); SPEI_ens=Miss ! Ens Mean ! Allocate (SPEI_ens1(Nlon,Nlat,Nlead,Ntime)); SPEI_ens1=Miss ! Ens Median Mvar=Nmodel*(Ncase_1-2) ! =2*2=4 (RegCM,clWRF, "ann","mlr") Do it=1,Ntime Y_inp(it) = SPEI_obs (ilon,ilat,it) j=0 Do imod=1,Nmodel Do icas=3,Ncase_1 ! Use only "ann" and "mlr" j=j+1 X_inp(it,j) = SPEI_M1(ilon,ilat,ilead,it,icas,imod) Enddo Enddo Enddo ! Xhq(:), Beta(:) Mhq = Mvar+1 ! Tong cac bien Doc lap va 1 bien phu thuoc Allocate (X2d(Ntime,Mhq),Check(Ntime)) icount=0; Check=0 Do i=1,Ntime xtmp(1:Mvar)=X_inp(i,1:Mvar); OK=Valid_Ar1(xtmp,Mvar,Miss) if (Y_inp(i)/=Miss.AND.OK) then icount=icount+1; Check(i)=1 X2d(icount,1)=Y_inp(i); X2d(icount,2:Mhq)=X_inp(i,1:Mvar) endif Enddo ! Whether do HQTT or not? if (icount>10) then Nhq = Mhq*icount ! N x M if (Allocated (Xorg)) Deallocate (Xorg) ! 1-D array X(N x M) Allocate (Xorg(Nhq)) i=0 Do j=1,Mhq Do it=1,icount ! Sample size i=i+1; Xorg(i)=X2d(it,j) ! Xorg(icount x Mhq) Enddo Enddo Nhq = icount ! Sample size Deallocate (X2d) Allocate (A(Mhq)) ! Regression Coefficiences Call REGRE(Xorg,Nhq,Mhq,A,S,Q,U,RR,F) if (.not.Allocated(Beta)) Allocate (Beta(Nlon,Nlat,Nlead,Mhq)) ! To be written out to file Beta=Miss Beta(ilon,ilat,ilead,1:Mhq) = A(1:Mhq) ! Luu he so HQ ! Tinh gia tri HQ: if (Allocated (Xhq)) Deallocate (Xhq); Allocate (Xhq(Ntime)); Xhq=Miss Do i=1,Ntime if (Check(i)==1) then tmp=0. Do j=2,Mhq tmp=tmp+A(j)*X_inp(i,j-1) Enddo Xhq(i)=A(1) + tmp endif Enddo Deallocate (A,Check) else if (Allocated (Xhq)) Deallocate (Xhq); Allocate (Xhq(Ntime)); Xhq=Miss Deallocate (Check) endif End Subroutine HQTT !------------------------------------------------------------------------------ Subroutine Reinit if (Allocated(T2_All)) Deallocate (T2_All) if (Allocated(Pr_All)) Deallocate (Pr_All) if (Allocated (Rain)) Deallocate (Rain) if (Allocated (T2m)) Deallocate (T2m) if (Allocated (PET)) Deallocate (PET) if (Allocated (xspi)) Deallocate (xspi) if (Allocated (Pr_Reg)) Deallocate (Pr_Reg) if (Allocated (T2_Reg)) Deallocate (T2_Reg) if (Allocated (Pet_Reg)) Deallocate (Pet_Reg) if (Allocated (beta_spi)) Deallocate (beta_spi) if (Allocated (gamm)) Deallocate (gamm) if (Allocated (pzero)) Deallocate (pzero) if (Allocated (SPEI_ens)) Deallocate (SPEI_ens) !if (Allocated (Lat)) Deallocate (Lat) if (Allocated (X1d)) Deallocate (X1d) Allocate (Rain(Ntime),T2m(Ntime),PET(Ntime),xspi(Ntime)) !DatBL(Ntime,Nst),DatTMP(Ntime,Nst)) Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth)) ! Allocate (SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen)) ! Nlon,Nlat,Ntime,Ncase_1,Nlen ! Allocate (Pr_Reg(Nreg,Nlead,Ncase_1,Ntime),T2_Reg(Nreg,Nlead,Ncase_1,Ntime)) ! Allocate (Pet_Reg(Nreg,Nlead,Ncase_1,Ntime)) ! Allocate (SPI_Reg(Nreg,Nlead,Ncase_1,Ntime,Nlen)) Allocate (X1d(Ntime)) End Subroutine Reinit !------------------------------------------------------------------------------ Integer Function Length_series (X,Ninp,Miss) Integer :: Ninp,i,j Real :: X(Ninp), Miss j=0 Do i=Ninp,1,-1 if (X(i)==Miss) then j=j+1 else exit endif Enddo !print*,Ninp,j Length_series=Ninp-j End Function Length_series !------------------------------------------------------------------------------ Subroutine Read_Regs Integer :: ilon,ilat, ireg, ic,j Real :: tmp Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg") Nst=0; Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Nst=Nst+1 endif Enddo Enddo !print*,"Nst=",Nst; pause End Subroutine Read_Regs !-------------------------------------------------------------------------------- Integer Function Get_Reg_Points (ireg) Integer :: ireg,ilon,ilat,Npoint Npoint=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==ireg) Npoint=Npoint+1 Enddo Enddo Get_Reg_Points = Npoint End Function Get_Reg_Points !------------------------------------------------------------------------------ Subroutine Write_Corre (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=3), parameter :: Modelname(2)=(/"reg",'clw'/) Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg Real :: xlead(Nlead),xcase(Ncase_1) Real, Allocatable :: X4d(:,:,:,:),Y4d(:,:,:,:) write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Corr_Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" Nz=Nmem*Ncase_1*Nmodel ! 7*3*2 = 42 ("cli","ann","mlr") if (Allocated(X4d)) Deallocate (X4d) if (Allocated(Y4d)) Deallocate (Y4d) Allocate (X4d(Nlon,Nlat,Ncase_1,Nlead),Y4d(Nlon,Nlat,Ncase_1,Nlead)) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="case"; t_nm="lead" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do i=1,Ncase_1; xcase(i)=i; Enddo Where (IsNaN(Corre)) Corre=Miss Where (IsNaN(Probe)) Probe=Miss Do imod=1,Nmodel Do icas=1,Ncase_1 Do ilead=1,Nlead X4d(:,:,icas,ilead) = Corre(:,:,imod,icas,ilead) Y4d(:,:,icas,ilead) = Probe(:,:,imod,icas,ilead) Enddo Enddo if (imod==1) then Var_ln="Correlation between "//trim(Modelname(imod))//" and obs for the cases of mod, cli, ann, and mlr" Var_nm=trim(Modelname(imod)) Call wr_ncdf_4d (Fname,Lon,Lat,xcase,xlead, Nlon,Nlat,Ncase_1,Nlead,X4d,trim(Var_nm), & t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="P-value of Correlation between "//trim(Modelname(imod))//" and obs for the cases of mod, cli, ann, and mlr" Var_nm="p_"//trim(Modelname(imod)) !//"_"//trim(Case_1(icas)) Call wr_ncdf_4d_apd (Fname,Y4d,trim(Var_nm),Nlon,Nlat,Ncase_1,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) else Var_ln="Correlation between "//trim(Modelname(imod))//" and obs for the cases of mod, cli, ann, and mlr" Var_nm=trim(Modelname(imod)) !//"_"//trim(Case_1(icas)) Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Ncase_1,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="P-value of Correlation between "//trim(Modelname(imod))//" and obs for the cases of mod, cli, ann, and mlr" Var_nm="p_"//trim(Modelname(imod)) Call wr_ncdf_4d_apd (Fname,Y4d,trim(Var_nm),Nlon,Nlat,Ncase_1,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo End Subroutine Write_Corre !------------------------------------------------------------------------------ Subroutine Calc_Correlation Integer :: it,i,ilen,ilon,ilat,ilead,imod,imem,icas, it1, it2, icpr, ict, icpet Real :: X1d(Ntime), Y1d(Ntime) Real :: cor,prob if (Allocated (Corre)) Deallocate (Corre) if (Allocated (Probe)) Deallocate (Probe) Allocate (Corre(Nlon,Nlat,Nmodel,Ncase_1,Nlead)) Allocate (Probe(Nlon,Nlat,Nmodel,Ncase_1,Nlead)) Corre=Miss; Probe=Miss imem=1 ! Use member 1 only Do imod=1,Nmodel Do icas=1,Ncase_1 Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Ntime X1d(it) = SPEI_obs(ilon,ilat,it) Y1d(it) = SPEI_All(ilon,ilat,imem,ilead,it,icas,imod) Enddo Call Correlation (X1d,Y1d,Ntime,Cor,prob,Miss) Corre(ilon,ilat,imod,icas,ilead) = Cor Probe(ilon,ilat,imod,icas,ilead) = Prob !print*,Corre(ilon,ilat,imod,icas,ilead) endif Enddo Enddo !print*,"Model:",imod," ",Case_1(icas) Enddo !pause Enddo Enddo End Subroutine Calc_Correlation !------------------------------------------------------------------------------ Subroutine Read_all_Dr_IDXs (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=20) :: Var_nm Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real, Allocatable :: X4d(:,:,:,:),xlon(:),xlat(:),xlead(:),xtime(:) if (Allocated (SPEI_M1)) Deallocate (SPEI_M1) Allocate (SPEI_M1(Nlon,Nlat,Nlead,Ntime,Ncase_1,Nmodel)) ! SPEI_obs(Nlon,Nlat,Ntime) write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" i=0 Do imod=1,Nmodel !Do imem=1,Nmem Do icas=1,Ncase_1 i=i+1 !write (Var_nm,"(A,I2.2)") "idx_",i Var_nm=trim(Model_nm(imod))//"_"trim(Case_1(icas)) Call rd_ncdf_4d (Fname,xlon,xlat,xlead,xtime,Nlo,Nla,Nle,Nt,X4d,trim(Var_nm)) ! X4d(lon,lat,lead,time) Do ilead=1,Nlead Do it=1,Ntime SPEI_M1(:,:,imem,ilead,it,icas,imod) = X4d(:,:,ilead,it) Enddo Enddo Deallocate (xlon,xlat,xlead,xtime,X4d) Enddo !Enddo Enddo if (Allocated(SPEI_obs)) Deallocate (SPEI_obs) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,SPEI_obs,"obs") Deallocate (xlon,xlat,xtime) Where (IsNaN(SPEI_M1)) SPEI_M1=Miss Where (IsNaN(SPEI_obs)) SPEI_obs=Miss Where (SPEI_M1>10.0) SPEI_M1=Miss Where (SPEI_obs>10.0) SPEI_obs=Miss End Subroutine Read_all_Dr_IDXs !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Sub_SPI_One_Stn_2023.f90" Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_2023.f90" Include "/home/tanpv/Mylibs/PET_thornthwaite.f90" Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_PET_2023.f90" !Include "/home/tanpv/Mylibs/Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI !Include "/home/tanpv/Mylibs/Sub_SPEI.f90" ! Chuong trinh tinh chi so han SPEI Include "/home/tanpv/Mylibs/Sub_Percentile.f90" Include "/home/tanpv/Mylibs/Sub_Verification.f90" Include "/home/tanpv/Mylibs/Sub_HQTT.f90" Include "/home/tanpv/Mylibs/Function_Valid_Ar1.f90" Include "/home/tanpv/Mylibs/Correlation.f90" !-------------------------------------------------------------------- END ! Program Include "/home/tanpv/Mylibs/Student_Prob.f90"