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 = 4 ! "obs","cli","ann","mlr" Integer, Parameter :: Ncase_1 = Ncase-1 ! "cli","ann","mlr" Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"cli","ann","mlr"/) Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","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 :: Agre_ens (:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Dr_Lev ! Agreement between model and obs for each month Real, Allocatable :: Agre_ens1(:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Dr_Lev Real, Allocatable :: Agre_ens_All(:,:,:,:) ! Nlon,Nlat,Nlead,Dr_Lev ! For all timeseries Real, Allocatable :: Agre_ens1_All(:,:,:,:) ! Nlon,Nlat,Nlead,Dr_Lev !------------------------------------------------------------------------------ ! 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) 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 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) Call Write_Pr_Temp_DrIDXs (trim(IDXname(idx)),ilen) ! (Nlon,Nlat,Nz,Nlead), Nz=2 models x 7 members x 4 cases = 56 Enddo Enddo !Call Write_Pr_Temp_DrIDXs ("spei") !------------------------------------------------------------------------------ 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 Agreement Integer :: ispl,ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: nct ! Real, Allocatable :: Agre_ens (:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Dr_Lev ! Agreement between model and obs for each month ! Real, Allocatable :: Agre_ens1(:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Dr_Lev ! Real, Allocatable :: AgAgre_ens_All(:,:,:,:) ! Nlon,Nlat,Nlead,Dr_Lev ! For all timeseries ! Real, Allocatable :: AgAgre_ens1_All(:,:,:,:) ! Nlon,Nlat,Nlead,Dr_Lev ! 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 ! Integer, Parameter :: Dr_Lev = 3 ! Real, parameter :: Dr_The(Dr_Lev)=(/-0.5, -1.0, -1.5/) ! 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.5) SPEI_ens = 0 Where (SPEI_ens1>-0.5) SPEI_ens1 = 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 <=-1.5.AND.SPEI_ens/=Miss) SPEI_ens = 3 Where (SPEI_ens1<=-1.5.AND.SPEI_ens1/=Miss) SPEI_ens1 = 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 <=-1.0.AND.SPEI_ens/=Miss) SPEI_ens = 2 Where (SPEI_ens1<=-1.0.AND.SPEI_ens1/=Miss) SPEI_ens1 = 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.5.AND.SPEI_ens/=Miss) SPEI_ens = 1 Where (SPEI_ens1<=-0.5.AND.SPEI_ens1/=Miss) SPEI_ens1 = 1 if (Allocated (Agre_ens)) Deallocate (Agre_ens) if (Allocated (Agre_ens1)) Deallocate (Agre_ens1) if (Allocated (Agre_ens_All)) Deallocate (Agre_ens_All) if (Allocated (Agre_ens1_All)) Deallocate (Agre_ens1_All) Allocate (Agre_ens(Nlon,Nlat,Nlead,Nmon,Dr_Lev), Agre_ens1(Nlon,Nlat,Nlead,Nmon,Dr_Lev)) Allocate (Agre_ens_All(Nlon,Nlat,Nlead,Dr_Lev), Agre_ens1_All(Nlon,Nlat,Nlead,Dr_Lev)) Agre_ens=Miss; Agre_ens1=Miss; Agre_ens_All=Miss; Agre_ens1_All=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ilead=1,Nlead Do iidx=1,Dr_Lev Do imon=1,Nmon nct=0; ict=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then endif Enddo Enddo Enddo Enddo endif Enddo Enddo End Subroutine Agreement !------------------------------------------------------------------------------ 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 ! Ens: SPEI_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime ! SPEI_ens1 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,Nmem 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_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_Pr_Temp_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)//"Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" Nz=Nmem*Ncase_1*Nmodel ! 7*3*2 = 42 ("cli","ann","mlr") 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="case"; t_nm="lead" 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,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 write (Var_ln,"(A,I2,A,I2,A)") "Values of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " of Model "//trim(Model_nm(imod))//" leadtime ",ilead," months, 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*Nmem*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*Nmem*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) ! 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 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_Pr_Temp_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 !------------------------------------------------------------------------------ Subroutine Write_Dr_IDXs_to_file (Fname,ilead) !,Model_nm) Character (Len=*) :: Fname !,Model_nm Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: it, imem, ilead, ireg, icas, ilen, imod Real :: zlev(Nlen) Real, Allocatable :: X3d(:,:,:) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="time" Do ilen=1,Nlen; zlev(ilen)=Len_SPI(ilen); Enddo Allocate (X4d(Nlon,Nlat,Nlen,Ntime)) ! SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen) ! Pr_ens (Nlon,Nlat,Nlead,Ntime,Ncase_1) Do icas=1,Ncase_1 Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat Do it=1,Ntime !X4d(ilon,ilat,ilen,it)=SPEI_ens(ilon,ilat,it,icas,ilen) Enddo Enddo Enddo Enddo write(Var_ln,"(A,I2.2,A,I1)") "SPEI for "//trim(Case_1(icas))//" with timescale in z=1,2,3 corresponding to 1,3,6 month(s) for Leadtime=",ilead write(Var_nm,"(A,I2.2)") trim(Case_1(icas)) !//"_",Len_SPI(ilen) if (icas==1) then Call wr_ncdf_4d (Fname, xlon,xlat,zlev,xtime, Nlon,Nlat,Nlen,Ntime,X4d,Var_nm, & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) else Call wr_ncdf_4d_apd (Fname,X4d,Var_nm,Nlon,Nlat,Nlen,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo Deallocate (X4d) ! Write out Pr and T2m: Allocate (X3d(Nlon,Nlat,Ntime)) X3d(:,:,1:Ntime)=Pr_ens (:,:,ilead,1:Ntime,1) ! Pr_obs Call wr_ncdf_3d_apd (Fname,X3d,"pr_obs",Nlon,Nlat,Ntime,"Pr obs","mm/month",lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Ntime)=Pr_ens (:,:,ilead,1:Ntime,2) ! Pr_model_ens Call wr_ncdf_3d_apd (Fname,X3d,"pr_ens",Nlon,Nlat,Ntime,"Pr model ens","mm/month",lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Ntime)=T2_ens (:,:,ilead,1:Ntime,1) ! T2_obs Call wr_ncdf_3d_apd (Fname,X3d,"t2_obs",Nlon,Nlat,Ntime,"T2 obs","C",lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Ntime)=T2_ens (:,:,ilead,1:Ntime,2) ! Pr_model_ens Call wr_ncdf_3d_apd (Fname,X3d,"t2_ens",Nlon,Nlat,Ntime,"T2 model ens","C",lon_nm,lat_nm,t_nm,Miss) Deallocate (X3d) End Subroutine Write_Dr_IDXs_to_file !------------------------------------------------------------------------------ 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)=(/"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; 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(Cases(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" !-------------------------------------------------------------------- END ! Program