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" Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"mod","cli","ann","mlr"/) 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 !------------------------------------------------------------------------------ ! 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 :: 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 Calc_IDXs (trim(IDXname(idx)),ilen) ! Outputs: SPEI_obs, SPEI_All,SPEI_ens,SPEI_ens Call Calc_IDX_Ens ! Outputs: SPEI_ens,SPEI_ens1 Call Calc_HQTT ! Outputs: SPEI_mlr. Based on SPEI_All Call Write_DrIDXs (trim(IDXname(idx)),ilen) ! Out: SPEI_All,SPEI_ens,SPEI_ens1,SPEI_mlr !!! ,Prob_mild,Prob_mode,Prob_seve ! 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 Test Integer :: Nspl, ilead,imod,ilon,ilat,imem,icas,ic Real, Allocatable :: X(:),Y(:) Real :: pobs,pens,tobs,tens,pper,tper ! Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel) Do imod=2,Nmodel Do ilead=1,Nlead-1 Do imem=1,Nmem print*,Pr_All(:,:,imem,ilead,:,1,imod)-Pr_All(:,:,imem,ilead+1,:,1,imod) print*,ilead,imem,imod pause Enddo Enddo Enddo !------- Nspl = Nmem*3*Nmodel Allocate (X(Nspl),Y(Nspl)) Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do it=1,Ntime if (ilead==1) then pobs=Pr_All(ilon,ilat,1,1,it,1,1); tobs=T2_All(ilon,ilat,1,1,it,1,1) endif ic=0 Do imod=1,Nmodel Do imem=1,Nmem-1 Do icas=3,Ncase ic=ic+1 X(ic)=Pr_All(ilon,ilat,imem,ilead,it,icas,imod) Y(ic)=T2_All(ilon,ilat,imem,ilead,it,icas,imod) Enddo Enddo Enddo pens=Ave(X,ic,Miss); tens=Ave(Y,ic,Miss) pper=percentile(X,ic,0.5,Miss); tper=percentile(Y,ic,0.5,Miss) ! if (pobs/=Miss.AND.pens/=Miss.AND.tobs/=Miss.AND.tens/=Miss) & ! print"(10F6.1)",pens-pobs,pens/pobs*100,tens-tobs,pper-pobs,pper/pobs*100,tper-tobs Enddo !pause endif Enddo Enddo Enddo End Subroutine Test !------------------------------------------------------------------------------ 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) 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, 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),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*4=8 (RegCM,clWRF, "mod","cli","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=1,Ncase_1 j=j+1 X_inp(it,j) = SPEI_All(ilon,ilat,imem,ilead,it,icas,imod) !X_inp(it,2) = SPEI_ens1(ilon,ilat,ilead,it) 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 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) End Subroutine HQTT !------------------------------------------------------------------------------ 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_mlr >-0.5) SPEI_mlr = 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_mlr <=-1.5.AND.SPEI_mlr /=Miss) SPEI_mlr = 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_mlr <=-1.0.AND.SPEI_mlr /=Miss) SPEI_mlr = 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_mlr <=-0.5.AND.SPEI_mlr /=Miss) SPEI_mlr = 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_mlr(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_mlr(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_ens(ilon,ilat,ilead,it)/=Miss) then igrdf=igrdf+1 if (SPEI_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 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),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_All_Reg: Do imod=1,Nmodel Do ilead=1,Nlead Do imem=1,Nmem Do icas=1,Ncase_1 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_All(ilon,ilat,imem,ilead,it,icas,imod) endif Enddo Enddo SPEI_All_Reg(ireg,imem,ilead,it,icas,imod)=percentile(xtmp,ispl,0.5,Miss) Enddo Enddo Enddo Enddo Enddo Enddo ! 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; xtmp1=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==Real(ireg)) then ispl=ispl+1 xtmp (ispl)=SPEI_ens (ilon,ilat,ilead,it) xtmp1(ispl)=SPEI_ens1(ilon,ilat,ilead,it) xtmp2(ispl)=SPEI_mlr(ilon,ilat,ilead,it) endif Enddo Enddo SPEI_ens_Reg (ireg,ilead,it) = percentile(xtmp, ispl,0.5,Miss) SPEI_ens1_Reg(ireg,ilead,it) = percentile(xtmp1,ispl,0.5,Miss) SPEI_mlr_Reg (ireg,ilead,it) = percentile(xtmp2,ispl,0.5,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) 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 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) Enddo endif Enddo Enddo 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_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 !Character (Len=3) :: Modelnm(2)=("reg","clw") 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)//"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_All)) SPEI_All=Miss Where (IsNaN(SPEI_obs)) SPEI_obs=Miss Where (IsNaN(SPEI_ens)) SPEI_ens=Miss i=0 Do imod=1,Nmodel Do imem=1,1 ! Nmem Do icas=1,Ncase_1 i=i+1 Do ilead=1,Nlead Do it=1,Ntime X4d(:,:,ilead,it)=SPEI_All(:,:,imem,ilead,it,icas,imod) Enddo Enddo !write (Var_nm,"(A,I2.2)") "idx_",i Var_nm = trim(Model_nm(imod))//"_"//Case_1(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))//", for z=1 to ",Nlead," months leadtime, and BC correction using "//Case_1(icas) 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), & " month. (Ensemble over Nmodel*Ncase_1)" Var_nm="ens" Call wr_ncdf_4d_apd (Fname,SPEI_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 model sensemble Median of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " month. (Ensemble over Nmodel*Ncase_1)" Var_nm="emed" Call wr_ncdf_4d_apd (Fname,SPEI_ens1,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 MLR based on all dr_idxs of mem1 of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " month. (Ensemble over Nmodel*Ncase_1)" Var_nm="mlr" Call wr_ncdf_4d_apd (Fname,SPEI_mlr,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)") "Probability of occurrence of Mild drought based on "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen)," months" !Var_nm="mild" !Call wr_ncdf_4d_apd (Fname,Prob_mild,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)") "Probability of occurrence of Moderate drought based on "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen)," months" !Var_nm="moder" !Call wr_ncdf_4d_apd (Fname,Prob_mode,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)") "Probability of occurrence of Seve drought based on "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen)," months" !Var_nm="seve" !Call wr_ncdf_4d_apd (Fname,Prob_seve,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) !Call wr_ncdf_4d (Fname,Lon,Lat,xcase,xlead, Nlon,Nlat,Nz,Nlead,X4d,"dridx", & ! t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) !Call wr_ncdf_4d_apd (Fname,Pet_Reg,"pet",Nreg,Nlead,Ncase_1,Ntime, & ! Var_ln,"mm/month",lon_nm, lat_nm, z_nm, t_nm,Miss) !----------------------------------------- Deallocate (X4d,xtime) End Subroutine Write_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 Read_all_data ! (Model_nm) !Character (Len=*) :: Model_nm Integer :: imod,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:) !,xtime(:) it1=Nmon+1 ! Tinh tu 1/1983 pr_clWRF_lead_01.nc ! Case_1(Ncase_1)=(/"mod","cli","ann","mlr"/) ! Ntime=(Yr2-Yr1)*Nmon - bo nam dau tien 1982 Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase_1,Nmodel)) Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase_1,Nmodel)) Allocate (Pr_obs(Nlon,Nlat,Ntime),T2_obs(Nlon,Nlat,Ntime)) T2_All=Miss; Pr_All=Miss; T2_obs=Miss; Pr_obs=Miss if (Allocated(X4d)) Deallocate (X4d); if (Allocated(xmem)) Deallocate (xmem) if (Allocated(xlon)) Deallocate (xlon); if (Allocated(xlat)) Deallocate (xlat) if (Allocated(time)) Deallocate (time); if (Allocated(xtime)) Deallocate (xtime) Do imod=1,Nmodel Do ilead=1,Nlead write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" print*,"Reading ",trim(Fname) if (imod==1.AND.ilead==1) then Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,"obs") Do it=1,Ntime T2_obs(:,:,it)=X4d(:,:,1,Nmon+it) Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) endif Do icas=1,Ncase_1 Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Case_1(icas))) Do imem=1,Nmem Do it=1,Ntime T2_All(:,:,imem,ilead,it,icas,imod)=X4d(:,:,imem,Nmon+it) ! Luu vao T2_All voi time tu Nmon+1 = it1 (Tinh tu 1/1983) Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(RVar)//"_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" if (imod==1.AND.ilead==1) then Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,"obs") Do it=1,Ntime Pr_obs(:,:,it)=X4d(:,:,1,Nmon+it) Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) endif Do icas=1,Ncase_1 Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Case_1(icas))) print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt Do imem=1,Nmem Do it=1,Ntime Pr_All(:,:,imem,ilead,it,icas,imod)=X4d(:,:,imem,Nmon+it) ! Luu vao Pr_All voi time tu Nmon+1 = it1 (Tinh tu 1/1983) Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo Print*,"Reading data ...LT=",ilead Enddo Enddo Print*,"Nt=",Nt Print*,"Reading data ... Done" End Subroutine Read_all_data !------------------------------------------------------------------------------ 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" !-------------------------------------------------------------------- END ! Program