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_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") imem=1 ! Use only mem=1 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_All(ilon,ilat,imem,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 Selected_IDXs Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: xtmp(Ncase_2) ! SPEI_sel((:,:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_2 ! SPEI_All(ilon,ilat,imem,ilead,1:Nt1,icas,imod) ! SPEI_ens (ilon,ilat,ilead,it) if (Allocated (SPEI_sel)) Deallocate (SPEI_sel) Allocate (SPEI_sel(Nlon,Nlat,Nlead,Ntime,Ncase_2)) SPEI_sel=Miss imem=1 ! Use only mem 1 Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Ntime SPEI_sel(ilon,ilat,ilead,it,1)=SPEI_All(ilon,ilat,imem,ilead,it,3,1) ! "an1" SPEI_sel(ilon,ilat,ilead,it,2)=SPEI_All(ilon,ilat,imem,ilead,it,4,1) ! "mr1" SPEI_sel(ilon,ilat,ilead,it,3)=SPEI_All(ilon,ilat,imem,ilead,it,3,2) ! "an2" SPEI_sel(ilon,ilat,ilead,it,4)=SPEI_All(ilon,ilat,imem,ilead,it,4,2) ! "mr2" SPEI_sel(ilon,ilat,ilead,it,5)=SPEI_ens (ilon,ilat,ilead,it) SPEI_sel(ilon,ilat,ilead,it,6)=SPEI_ens1(ilon,ilat,ilead,it) SPEI_sel(ilon,ilat,ilead,it,7)=SPEI_mlr (ilon,ilat,ilead,it) Enddo endif Enddo Enddo Enddo ! SPEI_sel_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime if (Allocated (SPEI_sel_ens)) Deallocate (SPEI_sel_ens) Allocate (SPEI_sel_ens(Nlon,Nlat,Nlead,Ntime)) SPEI_sel_ens=Miss Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Ntime ispl=0; Do icas=1,Ncase_2 ispl=ispl+1; xtmp(ispl)=SPEI_sel(ilon,ilat,ilead,it,icas) Enddo SPEI_sel_ens(ilon,ilat,ilead,it)=percentile(xtmp, ispl,0.5,Miss) Enddo endif Enddo Enddo Enddo End Subroutine Selected_IDXs !------------------------------------------------------------------------------ Subroutine Write_Verification (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str 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) Real, Allocatable :: X4d(:,:,:,:),xlev(:) Logical :: append write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Veri_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(xlev)) Deallocate (xlev) Allocate (X4d(Nlon,Nlat,Dr_Lev,Nlead),xlev(Dr_Lev)) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="dr_lev"; t_nm="lead" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do i=1,Dr_Lev; xlev(i)=i; Enddo write (Str,"(A,I2,A)") " of drought forecast by "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " month(s). z=1..3: Drought Level - mild, moderate, extreme ; t=1..6 months: Leadtime" append=.false. ! Writing to new file Do imon=1,Nmon X4d(:,:,:,:) = PCs(:,:,:,:,imon) ! Proportion Correction Var_ln="Proportion Correction (PC)"//trim(Str) write (Var_nm,"(A,I2.2)") "pc_",imon if (.not.append) then Call wr_ncdf_4d (Fname,Lon,Lat,xlev,xlead, Nlon,Nlat,Dr_Lev,Nlead,X4d,trim(Var_nm), & t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) append=.true. ! Append to current file else Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) endif X4d(:,:,:,:) = PODs(:,:,:,:,imon) ! Probability of Detection Var_ln="Probability of Detection (POD)"//trim(Str) write (Var_nm,"(A,I2.2)") "pod_",imon Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) X4d(:,:,:,:) = TSs(:,:,:,:,imon) ! Threat score (critical success index) Var_ln="Threat Score (TS)"//trim(Str) write (Var_nm,"(A,I2.2)") "ts_",imon Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Heidke Skill Score (HSS)"//trim(Str) write (Var_nm,"(A,I2.2)") "hss_",imon Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Enddo ! Allocate (PCall(Nlon,Nlat,Dr_Lev,Nlead), PODall(Nlon,Nlat,Dr_Lev,Nlead)) ! Allocate (TSall(Nlon,Nlat,Dr_Lev,Nlead), HSSall(Nlon,Nlat,Dr_Lev,Nlead)) ! PCall=Miss; PODall=Miss; TSall=Miss; HSSall=Miss Var_ln="Proportion Correction (PC) calculated over all timeseries "//trim(Str) Var_nm="pc" Call wr_ncdf_4d_apd (Fname,PCall,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Probability of Detection (POD) calculated over all timeseries "//trim(Str) Var_nm="pod" Call wr_ncdf_4d_apd (Fname,PODall,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Threat Score (TS) calculated over all timeseries"//trim(Str) Var_nm="ts" Call wr_ncdf_4d_apd (Fname,TSall,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Heidke Skill Score (HSS) calculated over all timeseries"//trim(Str) Var_nm="hss" Call wr_ncdf_4d_apd (Fname,HSSall,trim(Var_nm),Nlon,Nlat,Dr_Lev,Nlead,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Deallocate (X4d,xlev) End Subroutine Write_Verification !------------------------------------------------------------------------------ Subroutine Calc_Dr_events ! Set Dr_IDXs to 0 for all IDX values > -0.5 (no drought) Where (SPEI_obs >-0.5) SPEI_obs = 0 Where (SPEI_ens >-0.25) SPEI_ens = 0 Where (SPEI_ens1>-0.5) SPEI_ens1 = 0 Where (SPEI_sel_ens >-0.5) SPEI_sel_ens = 0 ! Set Dr_IDXs to 3 for all IDX values <= -1.5 (extreme drought) Where (SPEI_obs <=-1.5.AND.SPEI_obs/=Miss) SPEI_obs = 3 Where (SPEI_ens <=-0.75.AND.SPEI_ens/=Miss) SPEI_ens = 3 Where (SPEI_ens1<=-1.5.AND.SPEI_ens1/=Miss) SPEI_ens1 = 3 Where (SPEI_sel_ens <=-1.5.AND.SPEI_sel_ens /=Miss) SPEI_sel_ens = 3 ! Set Dr_IDXs to 2 for all IDX values <= -1.0 (moderate drought) Where (SPEI_obs <=-1.0.AND.SPEI_obs/=Miss) SPEI_obs = 2 Where (SPEI_ens <=-0.50.AND.SPEI_ens/=Miss) SPEI_ens = 2 Where (SPEI_ens1<=-1.0.AND.SPEI_ens1/=Miss) SPEI_ens1 = 2 Where (SPEI_sel_ens <=-1.0.AND.SPEI_sel_ens /=Miss) SPEI_sel_ens = 2 ! Set Dr_IDXs to 1 for all IDX values <= -0.5 (mild drought) Where (SPEI_obs <=-0.5.AND.SPEI_obs/=Miss) SPEI_obs = 1 Where (SPEI_ens <=-0.25.AND.SPEI_ens/=Miss) SPEI_ens = 1 Where (SPEI_ens1<=-0.5.AND.SPEI_ens1/=Miss) SPEI_ens1 = 1 Where (SPEI_sel_ens <=-0.5.AND.SPEI_sel_ens /=Miss) SPEI_sel_ens = 1 End Subroutine Calc_Dr_events !------------------------------------------------------------------------------ Subroutine Verification Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, ict1 Real :: nct, Cate, a,b,c,d Real :: O(Ntime), F(Ntime) ! ------------------------------------------------------------- if (Allocated (PCs)) Deallocate (PCs) if (Allocated (PODs)) Deallocate (PODs) if (Allocated (TSs)) Deallocate (TSs) if (Allocated (HSSs)) Deallocate (HSSs) if (Allocated (PCall)) Deallocate (PCall) if (Allocated (PODall)) Deallocate (PODall) if (Allocated (TSall)) Deallocate (TSall) if (Allocated (HSSall)) Deallocate (HSSall) ! ------------------------------------------------------------- Allocate (PCs(Nlon,Nlat,Dr_Lev,Nlead,Nmon), PODs(Nlon,Nlat,Dr_Lev,Nlead,Nmon)) Allocate (TSs(Nlon,Nlat,Dr_Lev,Nlead,Nmon), HSSs(Nlon,Nlat,Dr_Lev,Nlead,Nmon)) PCs=Miss; PODs=Miss; TSs=Miss; HSSs=Miss Allocate (PCall(Nlon,Nlat,Dr_Lev,Nlead), PODall(Nlon,Nlat,Dr_Lev,Nlead)) Allocate (TSall(Nlon,Nlat,Dr_Lev,Nlead), HSSall(Nlon,Nlat,Dr_Lev,Nlead)) PCall=Miss; PODall=Miss; TSall=Miss; HSSall=Miss ! ------------------------------------------------------------- Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ilead=1,Nlead Do iidx=1,Dr_Lev Cate = iidx ! Drought level Do imon=1,Nmon ict=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then ict=ict+1 O(ict)=SPEI_obs (ilon,ilat,it); F(ict)=SPEI_sel_ens(ilon,ilat,ilead,it) endif Enddo Call Calc_abcd (O,F,ict,Cate,a,b,c,d,Miss) PCs (ilon,ilat,iidx,ilead,imon)=PC (a,b,c,d) PODs(ilon,ilat,iidx,ilead,imon)=POD (a,b,c,d) TSs (ilon,ilat,iidx,ilead,imon)=TS (a,b,c,d) HSSs(ilon,ilat,iidx,ilead,imon)=HSS (a,b,c,d) Enddo ! All months Do it=1,Ntime O(it)=SPEI_obs (ilon,ilat,it); F(it)=SPEI_sel_ens(ilon,ilat,ilead,it) Enddo Call Calc_abcd (O,F,Ntime,Cate,a,b,c,d,Miss) PCall (ilon,ilat,iidx,ilead)=PC (a,b,c,d) PODall(ilon,ilat,iidx,ilead)=POD (a,b,c,d) TSall (ilon,ilat,iidx,ilead)=TS (a,b,c,d) HSSall(ilon,ilat,iidx,ilead)=HSS (a,b,c,d) Enddo Enddo endif Enddo Enddo End Subroutine Verification !------------------------------------------------------------------------------ Subroutine Write_GeoExt (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str 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),xreg(Nreg),xlev(Dr_Lev),xtime(Ntime) Logical :: append write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Geo_Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="reg"; lat_nm="dr_lev"; z_nm="lead"; t_nm="time" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do i=1,Dr_Lev; xlev(i)=i; Enddo Do i=1,Nreg; xreg(i)=i; Enddo Do i=1,Ntime; xtime(i)=i; Enddo Var_ln="Model drought geographic extent. reg=1..7 subregs; dr_lev=1..3: mild,moder,extreme; lead=1..6 months; time=1/1983-12/2020" Var_nm="geof" Call wr_ncdf_4d (Fname,xreg,xlev,xlead,xtime, Nreg,Dr_Lev,Nlead,Ntime,Geo_F,trim(Var_nm), & t_unit,Var_ln,"%",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Observed drought geographic extent. reg=1..7 subregs; dr_lev=1..3: mild,moder,extreme; lead=1..6 months; time=1/1983-12/2020" Var_nm="geob" Call wr_ncdf_3d_apd (Fname,Geo_O,trim(Var_nm),Nreg,Dr_Lev,Ntime,Var_ln,"%",lon_nm,lat_nm,t_nm,Miss) End Subroutine Write_GeoExt !------------------------------------------------------------------------------ Subroutine Calc_Geo_Extent Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, ict1 Real :: nct, Cate, a,b,c,d Real :: igrdf, igrdo,idrf, idro ! ------------------------------------------------------------- if (Allocated (Geo_F)) Deallocate (Geo_F) if (Allocated (Geo_O)) Deallocate (Geo_O) Allocate (Geo_F(Nreg,Dr_Lev,Nlead,Ntime),Geo_O(Nreg,Dr_Lev,Ntime)) Geo_F=Miss; Geo_O=Miss ! ------------------------------------------------------------- Do ilead=1,Nlead Do ilev=1,Dr_Lev Do it=1,Ntime Do ireg=1,Nreg igrdf=0; igrdo=0; idrf=0; idro=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then if (SPEI_sel_ens(ilon,ilat,ilead,it)/=Miss) then igrdf=igrdf+1 if (SPEI_sel_ens(ilon,ilat,ilead,it)==Real(ilev)) idrf=idrf+1 endif if (ilead==1) then if (SPEI_obs(ilon,ilat,it)/=Miss) then igrdo=igrdo+1 if (SPEI_obs(ilon,ilat,it)==Real(ilev)) idro=idro+1 endif endif endif Enddo Enddo if (igrdf>0) Geo_F(ireg,ilev,ilead,it)=idrf/igrdf*100.0 if (ilead==1) then if (igrdo>0) Geo_O(ireg,ilev,it)=idro/igrdo*100.0 endif Enddo Enddo Enddo Enddo End Subroutine Calc_Geo_Extent !------------------------------------------------------------------------------ Subroutine Calc_abcd (O,F,N,Category,a,b,c,d,Miss) Real :: a,b,c,d,Category,Miss ! Category=1,2,3 Integer :: N,i Real :: O(N),F(N) a=0; b=0; c=0; d=0 Do i=1,N if (O(i)/=Miss.AND.F(i)/=Miss) then if (O(i)==Category.AND.F(i)==Category) a=a+1 if (O(i)/=Category.AND.F(i)==Category) b=b+1 if (O(i)==Category.AND.F(i)/=Category) c=c+1 if (O(i)/=Category.AND.F(i)/=Category) c=c+1 endif Enddo End Subroutine Calc_abcd !------------------------------------------------------------------------------ Subroutine Calc_Dr_Reg ! Real, Allocatable :: SPEI_All_Reg(:,:,:,:,:,:) ! Nreg,Nmem,Nlead,Ntime,Ncase,Nmodel ! Real, Allocatable :: SPEI_ens_Reg(:,:,:) ! Nreg,Nlead,Ntime - Ensemble mean ! Real, Allocatable :: SPEI_ens1_Reg(:,:,:) ! Nreg,Nlead,Ntime - Ensemble Median ! SPEI_obs_Reg(:,:) ! Nreg,Ntime ! ! SPEI_sel((:,:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_2 Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, latx Real :: xtmp(Nlon*Nlat*Ncase_2),xtmp1(Nlon*Nlat),xtmp2(Nlon*Nlat) ! if (Allocated (SPEI_All_Reg)) Deallocate (SPEI_All_Reg) if (Allocated (SPEI_ens_Reg)) Deallocate (SPEI_ens_Reg) if (Allocated (SPEI_ens1_Reg)) Deallocate (SPEI_ens1_Reg) ! if (Allocated (SPEI_mlr_Reg)) Deallocate (SPEI_mlr_Reg) if (Allocated (SPEI_obs_Reg)) Deallocate (SPEI_obs_Reg) ! Allocate (SPEI_All_Reg(Nreg,Nmem,Nlead,Ntime,Ncase_1,Nmodel)) Allocate (SPEI_ens_Reg(Nreg,Nlead,Ntime)) Allocate (SPEI_ens1_Reg(Nreg,Nlead,Ntime)) ! Allocate (SPEI_mlr_Reg(Nreg,Nlead,Ntime)) Allocate (SPEI_obs_Reg(Nreg,Ntime)) SPEI_All_Reg=Miss; SPEI_ens_Reg=Miss; SPEI_ens1_Reg=Miss; SPEI_mlr_Reg=Miss; SPEI_obs_Reg=Miss ! SPEI_ens_Reg, and SPEI_ens1_Reg(Nreg,Nlead,Ntime): SPEI_ens(Nlon,Nlat,Nlead,Ntime) Do ilead=1,Nlead Do it=1,Ntime Do ireg=1,Nreg ispl=0; xtmp=Miss; Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==Real(ireg)) then Do icas=1,Ncase_2 ispl=ispl+1 xtmp (ispl)=SPEI_sel(ilon,ilat,ilead,it,icas) Enddo endif Enddo Enddo SPEI_ens_Reg (ireg,ilead,it) = Ave (xtmp, ispl,Miss) SPEI_ens1_Reg(ireg,ilead,it) = percentile(xtmp, ispl,0.4,Miss) Enddo Enddo Enddo ! SPEI_obs(Nlon,Nlat,Ntime): Do it=1,Ntime Do ireg=1,Nreg ispl=0; xtmp=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==Real(ireg)) then ispl=ispl+1 xtmp (ispl)=SPEI_obs(ilon,ilat,it) endif Enddo Enddo !SPEI_obs_Reg(ireg,it) = percentile(xtmp, ispl,0.5,Miss) SPEI_obs_Reg(ireg,it) = Ave (xtmp, ispl,Miss) Enddo Enddo End Subroutine Calc_Dr_Reg !------------------------------------------------------------------------------ Subroutine Write_Dr_Reg_ens (IDX_nm,ilen) Character (Len=*) :: IDX_nm 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 :: xreg(Nreg),xlead(Nlead) Real, Allocatable :: X3d(:,:,:),xtime(:) ! SPEI_ens1_Reg(Nreg,Nlead,Ntime) write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Reg_ens_Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" if (Allocated(xtime)) Deallocate (xtime) if (Allocated(X3d)) Deallocate (X3d) Allocate (xtime(Ntime)); Allocate (X3d(Nreg,Nlead,Ntime)) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="reg"; lat_nm="lead"; t_nm="time" Do i=1,Ntime; xtime(i)=i; Enddo Do i=1,Nreg; xreg(i) =i; Enddo Do i=1,Nlead; xlead(i) =i; Enddo Var_ln="Regional Median values of ensemble mean of drought index" ! Based on ensemble mean Call wr_ncdf_3d (Fname,xreg,xlead,xtime,Nreg,Nlead,Ntime,SPEI_ens_Reg,"ens",t_unit, & var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_nm="Regional Median values of ensemble median of drought index" ! Based on ensemble median Call wr_ncdf_3d_apd (Fname,SPEI_ens1_Reg,"ens1",Nreg,Nlead,Ntime,Var_ln,"-",lon_nm,lat_nm,t_nm,Miss) !Var_nm="Regional Median values of MLR of drought index" ! Based on ensemble MLR !Call wr_ncdf_3d_apd (Fname,SPEI_mlr_Reg,"mlr",Nreg,Nlead,Ntime,Var_ln,"-",lon_nm,lat_nm,t_nm,Miss) Do ireg=1,Nreg Do ilead=1,Nlead Do it=1,Ntime X3d(ireg,ilead,it) = SPEI_obs_Reg(ireg,it) Enddo Enddo Enddo Var_ln="Regional Median values of observed drought index" Call wr_ncdf_3d_apd (Fname,X3d,"obs",Nreg,Nlead,Ntime,Var_ln,"-",lon_nm,lat_nm,t_nm,Miss) End Subroutine Write_Dr_Reg_ens !------------------------------------------------------------------------------ Subroutine Calc_IDXs (IDX_nm,ilen) ! IDX_nm="spi", "spei" Character (Len=*) :: IDX_nm Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, latx ! Real :: xens(Nmodel*Nmem*Ncase_1) ! For calculating SPI, SPEI: T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel) ! Pr_All(:,:,imem,ilead,it,icas,imod) if (Allocated (SPEI_All)) Deallocate (SPEI_All) Allocate (SPEI_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase_1,Nmodel)) SPEI_All=Miss if (Allocated (SPEI_obs)) Deallocate (SPEI_obs) Allocate (SPEI_obs(Nlon,Nlat,Ntime)); SPEI_obs=Miss Nt=Ntime if (Allocated (Rain)) Deallocate (Rain) if (Allocated (T2m)) Deallocate (T2m) Allocate (Rain(Nt),T2m(Nt)) ! Obs drought IDX Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then use_baseline=.false.; xspi=Miss; latx=Lat(ilat) Rain(1:Nt)=Pr_obs(ilon,ilat,1:Nt) T2m (1:Nt)=T2_obs(ilon,ilat,1:Nt) Nt1=Min(Length_series (Rain,Nt,Miss), Length_series (T2m,Nt,Miss)) if (Nt1>0) then if (IDX_nm=="spi") then Call Cal_SPI_One_Stn (Rain,Nt1,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss) SPEI_obs(ilon,ilat,1:Nt1)=xspi(1:Nt1) else if (IDX_nm=="spei") then Call Calc_SPEI_One_Stn (T2m,Rain,latx,Nt1,Nmon,Len_SPI(ilen),logLogisticParams,use_baseline,xspi,Miss) SPEI_obs(ilon,ilat,1:Nt1)=xspi(1:Nt1) endif endif endif Enddo Enddo ! Model Drought IDX Do imod=1,Nmodel Do ilead=1,Nlead Nt=Ntime if (Allocated (Rain)) Deallocate (Rain) if (Allocated (T2m)) Deallocate (T2m) Allocate (Rain(Nt),T2m(Nt)) Do imem=1,Nmem Do icas=1,Ncase_1 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then use_baseline=.false.; xspi=Miss; latx=Lat(ilat) Rain(1:Nt)=Pr_All(ilon,ilat,imem,ilead,1:Nt,icas,imod) T2m (1:Nt)=T2_All(ilon,ilat,imem,ilead,1:Nt,icas,imod) Nt1=Min(Length_series (Rain,Nt,Miss), Length_series (T2m,Nt,Miss)) if (Nt1>0) then if (IDX_nm=="spi") then Call Cal_SPI_One_Stn (Rain,Nt1,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss) SPEI_All(ilon,ilat,imem,ilead,1:Nt1,icas,imod)=xspi(1:Nt1) else if (IDX_nm=="spei") then Call Calc_SPEI_One_Stn (T2m,Rain,latx,Nt1,Nmon,Len_SPI(ilen),logLogisticParams,use_baseline,xspi,Miss) SPEI_All(ilon,ilat,imem,ilead,1:Nt1,icas,imod)=xspi(1:Nt1) endif !SPEI_Out(ilon,ilat,imem,ilead,icas,imod,ilen)=SPEI_All(ilon,ilat,imem,ilead,Nt,icas,imod) endif endif Enddo Enddo Enddo Enddo Enddo Enddo ! Outputs: SPEI_obs, SPEI_All,SPEI_ens,SPEI_ens1,Prob_mild,Prob_mode,Prob_seve ! Ens: SPEI_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime ! SPEI_ens1 End Subroutine Calc_IDXs !------------------------------------------------------------------------------ Subroutine Calc_IDX_Ens Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, latx Real :: xens(Nmodel*Ncase_1) ! xens(Nmodel*Nmem*Ncase_1) if (Allocated(SPEI_ens)) Deallocate (SPEI_ens) Allocate (SPEI_ens(Nlon,Nlat,Nlead,Ntime)); SPEI_ens=Miss ! Ens Mean if (Allocated(SPEI_ens1)) Deallocate (SPEI_ens1) Allocate (SPEI_ens1(Nlon,Nlat,Nlead,Ntime)); SPEI_ens1=Miss ! Ens Median Do ilead=1,Nlead ist=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then ist=ist+1 Do it=1,Ntime ispl=0 ! 1->Nmodel*(Ncase_1-2)=4 Do imod=1,Nmodel Do imem=1,1 !!!!! Nmem <<<<<< Use only member 1 Do icas=2,Ncase_1 ispl=ispl+1 xens(ispl)=SPEI_All(ilon,ilat,imem,ilead,it,icas,imod) Enddo Enddo Enddo !print"(20F7.3)",xens(1:ispl) SPEI_ens (ilon,ilat,ilead,it)=Ave (xens,ispl,Miss) SPEI_ens1(ilon,ilat,ilead,it)=percentile(xens,ispl,0.5,Miss) !print*,ist,SPEI_ens (ilon,ilat,ilead,it),SPEI_ens1(ilon,ilat,ilead,it),ispl Enddo endif Enddo Enddo !print*,"SPEI_ens" !pause Enddo End Subroutine Calc_IDX_Ens !------------------------------------------------------------------------------ Subroutine Calc_IDX_Prob Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, latx Real :: xens(Nmodel*Ncase_1) ! xens(Nmodel*Nmem*Ncase_1) !if (Allocated(SPEI_ens)) Deallocate (SPEI_ens) !Allocate (SPEI_ens(Nlon,Nlat,Nlead,Ntime)); SPEI_ens=Miss ! Ens Mean !if (Allocated(SPEI_ens1)) Deallocate (SPEI_ens1) !Allocate (SPEI_ens1(Nlon,Nlat,Nlead,Ntime)); SPEI_ens1=Miss ! Ens Median ! Calc probability: if (Allocated(Prob_mild)) Deallocate (Prob_mild) if (Allocated(Prob_mode)) Deallocate (Prob_mode) if (Allocated(Prob_seve)) Deallocate (Prob_seve) Allocate (Prob_mild(Nlon,Nlat,Nlead,Ntime)); Prob_mild=Miss Allocate (Prob_mode(Nlon,Nlat,Nlead,Ntime)); Prob_mode=Miss Allocate (Prob_seve(Nlon,Nlat,Nlead,Ntime)); Prob_seve=Miss Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Ntime ispl=0 ! 1->Nmodel*Nmem*Ncase_1=42 Do imod=1,Nmodel Do imem=1,1 !!!!! Nmem <<<<<< Use only member 1 Do icas=1,Ncase_1 ispl=ispl+1 xens(ispl)=SPEI_All(ilon,ilat,imem,ilead,it,icas,imod) Enddo Enddo Enddo !SPEI_ens (ilon,ilat,ilead,it)=Ave(xens,ispl,Miss) !SPEI_ens1(ilon,ilat,ilead,it)=percentile(xens,ispl,0.5,Miss) Prob_mild(ilon,ilat,ilead,it)=Probab(xens,ispl,Dr_The(1),.True.,Miss) Prob_mode(ilon,ilat,ilead,it)=Probab(xens,ispl,Dr_The(2),.True.,Miss) Prob_seve(ilon,ilat,ilead,it)=Probab(xens,ispl,Dr_The(3),.True.,Miss) Enddo endif Enddo Enddo Enddo End Subroutine Calc_IDX_Prob ! Calc_IDXs !------------------------------------------------------------------------------ Real Function Probab (X,N,Thre,Below,Miss) ! Tinh xac suat xuat hien P(X<=Thre) neu Below=.true., P(X>=Thre) neu Below=.false. Integer :: N,i Real :: X(N), Thre, Miss, tmp, ntmp, ic Logical :: Below tmp=0.; ntmp=0; ic=0 Do i=1,N if (X(i)/=Miss) then ntmp=ntmp+1 ! Num of valid values if (Below) then if (X(i)<=Thre) ic=ic+1 else if (X(i)>=Thre) ic=ic+1 endif endif Enddo if (ntmp>0) then Probab=ic/ntmp else Probab=Miss endif End Function Probab !------------------------------------------------------------------------------ Subroutine Write_Add_DrIDXs (IDX_nm,ilen) Character (Len=*) :: IDX_nm 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) Real, Allocatable :: X4d(:,:,:,:),xtime(:) ! Write out for testing: SPEI_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase_1,Nmodel) ! SPEI_obs(Nlon,Nlat,Ntime); SPEI_ens(Nlon,Nlat,Nlead,Ntime) ! X4d(Nlon,Nlat,Nz,Nlead): Nz=Nmodel*Nmem*Ncase !Fname="Areal_mean_Pr_T2_"//trim(IDX_nm)//"_"//trim(Model_nm)//".nc" write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Ens_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(xtime)) Deallocate (xtime) Allocate (X4d(Nlon,Nlat,Nlead,Ntime),xtime(Ntime)) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="lead"; t_nm="time" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do i=1,Ntime; xtime(i)=i; Enddo Where (IsNaN(SPEI_sel)) SPEI_sel=Miss Where (IsNaN(SPEI_sel_ens)) SPEI_sel_ens=Miss i=0 !Do imod=1,Nmodel ! Do imem=1,Nmem Do icas=1,Ncase_2 Do ilead=1,Nlead Do it=1,Ntime X4d(:,:,ilead,it)=SPEI_sel(:,:,ilead,it,icas) Enddo Enddo Var_nm=Case_2(icas) write (Var_ln,"(A,I2,A,I2,A)") "Values of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " month(s) of Model "//trim(Model_nm(imod)) i=i+1 if (i==1) then Call wr_ncdf_4d (Fname,Lon,Lat,xlead,xtime, Nlon,Nlat,Nlead,Ntime,X4d,trim(Var_nm), & t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) else Call wr_ncdf_4d_apd (Fname,X4d,trim(Var_nm),Nlon,Nlat,Nlead,Ntime,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo !Enddo !Enddo !----------------------------------------------------------- write(Var_ln,"(A,I2.2,A)")"Values of model sensemble mean of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen) Var_nm="en_all" Call wr_ncdf_4d_apd (Fname,SPEI_sel_ens,trim(Var_nm),Nlon,Nlat,Nlead,Ntime,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) write(Var_ln,"(A,I2.2,A)")"Values of oserved "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen)," month." Var_nm="obs" Call wr_ncdf_3d_apd (Fname,SPEI_obs,trim(Var_nm),Nlon,Nlat,Ntime,Var_ln,"-",lon_nm,lat_nm,t_nm,Miss) Deallocate (X4d,xtime) End Subroutine Write_Add_DrIDXs !------------------------------------------------------------------------------ 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_All)) Deallocate (SPEI_All) Allocate (SPEI_All(Nlon,Nlat,Nmem,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 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_All(:,:,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_All)) SPEI_All=Miss Where (IsNaN(SPEI_obs)) SPEI_obs=Miss Where (SPEI_All>10.0) SPEI_All=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"