Use NetCDF Real, Parameter :: Miss=-99. !Integer, Parameter :: Yr1=1986, Yr2=2005, Integer, Parameter :: Nmon=12, Nreg=7 ! D0 (abnormally dry, 20-30% percentile), D1 (moderate drought, 10-20%), D2 (severe drought, 5-10%), ! D3 (extreme drought, 2-5%) and D4 (exceptional drought, <2%) Real, Parameter :: pd0=0.30, pd1=0.20, pd2=0.10, pd3=0.05, pd4=0.02 ! Probability to calculate percentiles Integer, parameter :: NGeo_bin=14 Real, Parameter :: Geo_bin(0:NGeo_bin)=(/ 0, 2, 4, 6, 10, 14, 18, 24, 30, 36, 42, 50, 60, 76, 101/) Real, Parameter :: XGeo_bin(NGeo_bin) =(/ 1, 3, 5, 8, 12, 16, 21, 27, 33, 39, 46, 55, 68, 80/) ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime Logical :: Status Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:) ! (Nlon,Nlat,Ntime) Real, allocatable :: Mask(:,:), Mask_all(:,:) ! (Nlon,Nlat) Real, allocatable :: Dr_IDX(:,:,:,:) ! (Nlon,Nlat,Ntime,NDr_idx) Real, allocatable :: Qp(:,:,:) ! (Nlon,Nlat,NDr_idx) !--------------------------------------------------------------------------------- Real, allocatable :: Events(:,:,:,:) ! (Nlon,Nlat,Ntime,NDr_idx): =1 if Dr_IDX < Qp Real, allocatable :: FrMo(:,:,:,:) ! (Nlon,Nlat,Nmon,NDr_idx)-% of Event/mon Real, allocatable :: Freq(:,:,:,:) ! (Nlon,Nlat,Nyear,NDr_idx)-Num of Event/Year Real, allocatable :: Geo(:,:) ! (Ntime,NDr_idx) - Geo Ext in each mon Real, allocatable :: GeoReg(:,:,:) ! (Ntime,NDr_idx,Nreg) - Geo Ext in each mon Real, allocatable :: Dur_M(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Duration Real, allocatable :: Dur_X(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Duration Real, allocatable :: Sev_M(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Severity Real, allocatable :: Sev_X(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Severity Real, allocatable :: Inten_M(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Intensity Real, allocatable :: Inten_X(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Intensity Real, allocatable :: Inter_M(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Inter-arrival Real, allocatable :: Inter_X(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Inter-arrival ! New define for Severity and Intensity: ! Severity = Sum of DrIDX-Threshold; Intensity = Severity/Duration: Real, allocatable :: SEVm(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Severity Real, allocatable :: SEVx(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Severity Real, allocatable :: INTEm(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Mean Intensity Real, allocatable :: INTEx(:,:,:) ! ((Nlon,Nlat,NDr_idx) - Max Intensity Real, allocatable :: Fr_Max(:,:,:) ! (Nlon,Nlat,NDr_idx) - Maximum annual Freq (mon/year) Real, allocatable :: Fr_Std(:,:,:) ! (Nlon,Nlat,NDr_idx) - Sdt annual Freq (mon/year) !Real, allocatable :: q25(:,:,:) ! (Nlon,Nlat,NDr_idx) - Freq q25 (mon/year) !Real, allocatable :: q75(:,:,:) ! (Nlon,Nlat,NDr_idx) - Freq q75 (mon/year) !Real, allocatable :: IQR(:,:,:) ! (Nlon,Nlat,NDr_idx) - Freq IQR (mon/year) Real, allocatable :: Geo_PDF(:,:,:) ! (NDr_idx,Nreg,NGeo_bin) PDF Real, allocatable :: GeoRegMon(:,:,:) ! (Ntime,NDr_idx,Nreg+1) - Geo Ext in each mon Integer :: Nst, iidx, Yr1, Yr2, Nyear, iper,isce,imod, ipe, isc Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"pdsi","spi12","spe12"/) !(/"PDSI","SPI12","SPE12"/) Integer, Parameter :: NPeriod=3, NSce=2, Nmodels=6 ! Chi tinh cho 6-models de tao ra Ens Integer, Parameter :: NChars=12, NPeSc=5 ! 12 dac trung, 5=BL +4 SCe (MF45,MF85,FF45,FF85) Character (Len=100) :: Fname, LandMaskFile Character (Len=100) :: Data_Dir, Model_nm(Nmodels), Model_Obs, Out_Dir Integer :: YSta(NPeriod), YEnd(NPeriod) Character (Len=5) :: Sce_nm(NSce), Period(NPeriod) !--------------------------------------------------------------------------------- ! For Ensemble: Real, allocatable, Dimension (:,:,:,:,:,:) :: & ! (NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels DuM, Dux, SVm1, SVx1, SVm2, SVx2, ITm1, ITx1, ITm2, ITx2, IRm, IRx Real, allocatable, Dimension (:,:,:,:,:,:,:) :: FrM ! (NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmon,Nmodels) Real, allocatable, Dimension (:,:,:,:,:,:,:) :: FrY ! (NDr_idx,NPeriod,NSce,Nlon,Nlat,Nyear,Nmodels) Real, allocatable, Dimension (:,:,:,:,:) :: GeoVN ! (NDr_idx,NPeriod,NSce,Ntime,Nmodels) Real, allocatable, Dimension (:,:,:,:,:,:) :: GeoRg ! (NDr_idx,NPeriod,NSce,Nreg,Ntime,Nmodels) Real, allocatable :: XChars(:,:,:,:,:) ! (NDr_idx,Nlon,Nlat,NChars,5); 5=BL+2 Sce*2 Per Real, allocatable :: XFrm(:,:,:,:,:) ! NDr_idx,Nlon,Nlat,Nmon,5 Real, allocatable :: XFry(:,:,:,:,:) ! NDr_idx,Nlon,Nlat,Nyear,5 Real, allocatable :: GeoRgVN(:,:,:,:) ! NDr_idx,Nreg+1,Ntime,5 Real :: X1D(Nmodels) ! Working array !--------------------------------------------------------------------------------- Namelist /Infor/ Data_Dir, Period, YSta, YEnd, Model_nm, Sce_nm, LandMaskFile, Out_Dir Open (1, file="Namelist_Chars_Ens.txt", Status="Old") Read (1, Infor) Close (1) !---------------------------------------------- Call Read_Land_Mask Do iper=1, NPeriod Yr1=YSta(iper); Yr2=YEnd(iper) if (iper==1) then ! BL: 6 models+Obs, No Scenarios Do imod=1, Nmodels Fname=trim(Period(iper))//"_"//trim(Model_nm(imod))//".nc" Call Calc_Chars_one_Case (iper,1,imod) ! isce=1, iper=1: BL Enddo else Do isce=1,NSce Do imod=1, Nmodels Fname="SCE_"//trim(Period(iper))//"_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//".nc" Call Calc_Chars_one_Case (iper,isce,imod) !(trim(Data_Dir)//Fname) Enddo Enddo endif Enddo !---------------------------------------------- ! Calculate Ensemble: NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels Allocate (XChars(NDr_idx,Nlon,Nlat,NChars,NPeSc)) ! NChars=12, NPeSc=5 Allocate (XFrm(NDr_idx,Nlon,Nlat,Nmon,NPeSc)) Allocate (XFry(NDr_idx,Nlon,Nlat,Nyear,NPeSc)) Allocate (GeoRgVN(NDr_idx,Nreg+1,Ntime,NPeSc)) Do iidx=1,NDr_idx ipe=0 Do iper=1,NPeriod if (iper==1) then ! iper=1, isce=1: BL ipe=ipe+1 Do ilon=1,Nlon Do ilat=1,Nlat X1D(1:Nmodels)=DuM(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,1,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=DuX(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,2,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVm1(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,3,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVx1(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,4,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVm2(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,5,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVx2(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,6,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITm1(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,7,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITx1(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,8,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITm2(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,9,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITx2(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,10,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=IRm(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,11,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=IRx(iidx,iper,1,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,12,ipe)=Ave(X1D,Nmodels,Miss) Do imon=1,Nmon X1D(1:Nmodels)=FrM(iidx,iper,1,ilon,ilat,imon,1:Nmodels); XFrm(iidx,ilon,ilat,imon,ipe)=Ave(X1D,Nmodels,Miss) Enddo Do iyr=1,Nyear X1D(1:Nmodels)=FrY(iidx,iper,1,ilon,ilat,iyr,1:Nmodels); XFry(iidx,ilon,ilat,iyr,ipe)=Ave(X1D,Nmodels,Miss) Enddo Enddo Enddo Do it=1,Ntime X1D(1:Nmodels)=GeoVN(iidx,iper,1,it,1:Nmodels); GeoRgVN(iidx,8,it,ipe)=Ave(X1D,Nmodels,Miss) Do ireg=1,Nreg X1D(1:Nmodels)=GeoRg(iidx,iper,1,ireg,it,1:Nmodels); GeoRgVN(iidx,ireg,it,ipe)=Ave(X1D,Nmodels,Miss) Enddo Enddo else ! iper =2,3; isce=1,2: Scenarios Do isce=1,NSce ipe=ipe+1 Do ilon=1,Nlon; Do ilat=1,Nlat X1D(1:Nmodels)=DuM(iidx,iper,isce,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,1,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=DuX(iidx,iper,isce,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,2,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVm1(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,3,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVx1(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,4,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVm2(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,5,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=SVx2(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,6,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITm1(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,7,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITx1(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,8,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITm2(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,9,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=ITx2(iidx,iper,isce,ilon,ilat,1:Nmodels);XChars(iidx,ilon,ilat,10,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=IRm(iidx,iper,isce,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,11,ipe)=Ave(X1D,Nmodels,Miss) X1D(1:Nmodels)=IRx(iidx,iper,isce,ilon,ilat,1:Nmodels); XChars(iidx,ilon,ilat,12,ipe)=Ave(X1D,Nmodels,Miss) Do imon=1,Nmon X1D(1:Nmodels)=FrM(iidx,iper,isce,ilon,ilat,imon,1:Nmodels);XFrm(iidx,ilon,ilat,imon,ipe)=Ave(X1D,Nmodels,Miss) Enddo Do iyr=1,Nyear X1D(1:Nmodels)=FrY(iidx,iper,isce,ilon,ilat,iyr,1:Nmodels);XFry(iidx,ilon,ilat,iyr,ipe)=Ave(X1D,Nmodels,Miss) Enddo Enddo; Enddo Do it=1,Ntime X1D(1:Nmodels)=GeoVN(iidx,iper,isce,it,1:Nmodels); GeoRgVN(iidx,8,it,ipe)=Ave(X1D,Nmodels,Miss) Do ireg=1,Nreg X1D(1:Nmodels)=GeoRg(iidx,iper,isce,ireg,it,1:Nmodels);GeoRgVN(iidx,ireg,it,ipe)=Ave(X1D,Nmodels,Miss) Enddo Enddo Enddo endif Enddo Enddo Call Write_Char_Ens_to_NC !---------------------------------------------- CONTAINS !--------------------------------------------------------------------------------- Subroutine Read_Land_Mask LandMaskFile="/home/tanpv/Data/AWC/Land_Mask_VN.nc" Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"all") ! all=2..8=7 regs Allocate (Mask_all(nlon,nlat)) Mask_all=Mask End Subroutine Read_Land_Mask !--------------------------------------------------------------------------------- Subroutine Calc_Chars_one_Case (iper,isce,imod) !(Fname) Character (Len=100) :: Fout Integer :: iper,isce,imod Integer :: iidx,ilon,ilat,it Call Read_Drough_IDX (trim(Data_Dir)//trim(Fname)) Nyear = Ntime/Nmon ! New: For Ensemble: if (.not.Allocated(DuM)) Allocate (DuM(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(DuX)) Allocate (DuX(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(SVm1)) Allocate (SVm1(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(SVx1)) Allocate (SVx1(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(SVm2)) Allocate (SVm2(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(SVx2)) Allocate (SVx2(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(ITm1)) Allocate (ITm1(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(ITx1)) Allocate (ITx1(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(ITm2)) Allocate (ITm2(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(ITx2)) Allocate (ITx2(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(IRm)) Allocate (IRm(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(IRx)) Allocate (IRx(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmodels)) if (.not.Allocated(FrM)) Allocate (FrM(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nmon,Nmodels)) if (.not.Allocated(FrY)) Allocate (FrY(NDr_idx,NPeriod,NSce,Nlon,Nlat,Nyear,Nmodels)) if (.not.Allocated(GeoVN)) Allocate (GeoVN(NDr_idx,NPeriod,NSce,Ntime,Nmodels)) if (.not.Allocated(GeoRg)) Allocate (GeoRg(NDr_idx,NPeriod,NSce,Nreg,Ntime,Nmodels)) !----------------------------------------- Allocate (Qp(Nlon,Nlat,NDr_idx), Events(Nlon,Nlat,Ntime,NDr_idx)) Events (:,:,:,:) = Miss Call Calc_Per_Events (pd3) ! 0.05 Allocate (Freq(Nlon,Nlat,Nyear,NDr_idx),FrMo(Nlon,Nlat,Nmon,NDr_idx)) Allocate (Geo(Ntime,NDr_idx), GeoReg(Ntime,NDr_idx,Nreg)) Allocate (Dur_M(Nlon,Nlat,NDr_idx), Dur_X(Nlon,Nlat,NDr_idx)) Allocate (Sev_M(Nlon,Nlat,NDr_idx), Sev_X(Nlon,Nlat,NDr_idx)) Allocate (Inten_M(Nlon,Nlat,NDr_idx), Inten_X(Nlon,Nlat,NDr_idx)) Allocate (Inter_M(Nlon,Nlat,NDr_idx), Inter_X(Nlon,Nlat,NDr_idx)) Allocate (SEVm(Nlon,Nlat,NDr_idx),SEVx(Nlon,Nlat,NDr_idx)) Allocate (INTEm(Nlon,Nlat,NDr_idx), INTEx(Nlon,Nlat,NDr_idx)) Allocate (Fr_max(Nlon,Nlat,NDr_idx), Fr_Std(Nlon,Nlat,NDr_idx)) Allocate (Geo_PDF(NDr_idx,Nreg,NGeo_bin)) ! PDF of Geo Allocate (GeoRegMon(Nmon,NDr_idx,Nreg+1)) ! Including VN ! Set all to Miss: Freq=Miss; FrMo=Miss; Geo=Miss; GeoReg=Miss; Dur_M=Miss; Dur_X=Miss Sev_M=Miss; Sev_X=Miss; Inten_M=Miss; Inten_X=Miss; Inter_M=Miss; Inter_X=Miss Fr_Max=Miss; Fr_Std=Miss !q25=Miss; q75=Miss; IQR=Miss Geo_PDF=Miss; GeoRegMon=Miss SEVm=Miss; SEVx=Miss; INTEm=Miss; INTEx=Miss Call Calc_Chars Call Write_Fre_to_NC (trim(Out_Dir)//"Freq_"//trim(Fname),trim(Out_Dir)//"Fre_mon_"//trim(Fname)) Call Write_Char_to_NC(trim(Out_Dir)//"Chars_"//trim(Fname),trim(Out_Dir)//"Geo_"//trim(Fname), & trim(Out_Dir)//"Geo_PDF_"//trim(Fname), & trim(Out_Dir)//"Geo_Mon_"//trim(Fname)) ! For Ensemble: !--------------------------------------- iper,isce,imod Do iidx=1,NDr_idx Do ilon=1,Nlon Do ilat=1,Nlat DuM(iidx,iper,isce,ilon,ilat,imod) =Dur_M(ilon,ilat,iidx) DuX(iidx,iper,isce,ilon,ilat,imod) =Dur_X(ilon,ilat,iidx) SVm1(iidx,iper,isce,ilon,ilat,imod)=Sev_M(ilon,ilat,iidx) SVx1(iidx,iper,isce,ilon,ilat,imod)=Sev_X(ilon,ilat,iidx) SVm2(iidx,iper,isce,ilon,ilat,imod)=SEVm (ilon,ilat,iidx) SVx2(iidx,iper,isce,ilon,ilat,imod)=SEVx (ilon,ilat,iidx) ITm1(iidx,iper,isce,ilon,ilat,imod)=Inten_M(ilon,ilat,iidx) ITx1(iidx,iper,isce,ilon,ilat,imod)=Inten_X(ilon,ilat,iidx) ITm2(iidx,iper,isce,ilon,ilat,imod)=INTEm(ilon,ilat,iidx) ITx2(iidx,iper,isce,ilon,ilat,imod)=INTEx(ilon,ilat,iidx) IRm (iidx,iper,isce,ilon,ilat,imod)=Inter_M(ilon,ilat,iidx) IRx (iidx,iper,isce,ilon,ilat,imod)=Inter_X(ilon,ilat,iidx) FrM (iidx,iper,isce,ilon,ilat,1:Nmon,imod)=FrMo(ilon,ilat,1:Nmon,iidx) FrY (iidx,iper,isce,ilon,ilat,1:Nyear,imod)=Freq(ilon,ilat,1:Nyear,iidx) Enddo Enddo Do it=1,Ntime GeoVN(iidx,iper,isce,it,imod)=Geo(it,iidx) GeoRg(iidx,iper,isce,1:Nreg,it,imod)=GeoReg(it,iidx,1:Nreg) Enddo Enddo Deallocate (Qp,Events) Deallocate (Freq,FrMo,Geo,GeoReg,Dur_M,Dur_X,Sev_M,Sev_X,Inten_M,Inten_X,Fr_max,Fr_Std) Deallocate (SEVm,SEVx,INTEm,INTEx) Deallocate (Geo_PDF,GeoRegMon) Deallocate (Dr_IDX, Inter_M,Inter_X) End Subroutine Calc_Chars_one_Case !--------------------------------------------------------------------------------- Subroutine Write_Char_Ens_to_NC !(Fname,Fn_geo,Fn_geo_pdf,Fn_geo_mon) Character (Len=255) :: Fname, Fn_geo,Fn_geo_pdf,Fn_geo_mon Character (Len=50) :: var_n, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=250) :: var_ln Real :: xtime(Ntime), xreg(Nreg+1), xmon(Nmon), xyear(Nyear), xsce(NPeSc), xcha(NChars) Real, Allocatable :: X4D(:,:,:,:), X3D(:,:,:) Integer :: i, j, ipe, imon, it, ireg, iidx, iyr lon_nm="lon"; lat_nm="lat"; t_unit="years since 1987-01-01 " Do i=1,NChars; xcha(i)=i; Enddo Do i=1,NPeSc; xsce(i)=i; Enddo Do i=1,Nmon; xmon(i)=i; Enddo Do i=1,Nyear; xyear(i)=i; Enddo Do i=1,Ntime; xtime(i)=i; Enddo ! Chars files: Allocate (X4D(Nlon,Nlat,NPeSc,NChars)) var_ln="Ensemble mean of Drought Characteristics (Z: BL,MF45,MF85,FF45,FF85); (t=12: "// & "1:DuM,2:DuX,3:SVm1,4:SVx1,5:SVm2,6:SVx2,7:ITm1,8:ITx1,9:ITm2,10:ITx2,11:IRm,12:IRx)" z_nm="sce"; t_nm="cha"; Do iidx=1,NDr_idx Fname=trim(Out_Dir)//"Chars_"//trim(IDX_nm(iidx))//"_All_Ens.nc" Do ilon=1,Nlon Do ilat=1,Nlat Do ipe=1,NPeSc Do j=1,NChars X4D(ilon,ilat,ipe,j)=XChars(iidx,ilon,ilat,j,ipe) Enddo Enddo Enddo Enddo Call wr_ncdf_4d (Fname,xlon,xlat,xsce,xcha,nlon,nlat,NPeSc,NChars,X4D,"drcha", & t_unit, var_ln,"-", lon_nm, lat_nm, z_nm, t_nm, Miss) Enddo Deallocate (X4D) ! Monthly Frequency: Allocate (X4D(Nlon,Nlat,NPeSc,Nmon)) var_ln="Ensemble mean of Monthly Drought frequency (Z: BL,MF45,MF85,FF45,FF85)" z_nm="sce"; t_nm="mon"; Do iidx=1,NDr_idx Fname=trim(Out_Dir)//"Freq_Mon_"//trim(IDX_nm(iidx))//"_All_Ens.nc" Do ilon=1,Nlon Do ilat=1,Nlat Do ipe=1,NPeSc X4D(ilon,ilat,ipe,1:Nmon)=XFrm(iidx,ilon,ilat,1:Nmon,ipe) Enddo Enddo Enddo Call wr_ncdf_4d (Fname,xlon,xlat,xsce,xmon,nlon,nlat,NPeSc,Nmon,X4D,"fr_m", & t_unit, var_ln,"-", lon_nm, lat_nm, z_nm, t_nm, Miss) Enddo Deallocate (X4D) ! Yearly Frequency: Allocate (X4D(Nlon,Nlat,NPeSc,Nyear)) var_ln="Ensemble mean of annual Drought frequency (Z: BL,MF45,MF85,FF45,FF85)" t_nm="year"; Do iidx=1,NDr_idx Fname=trim(Out_Dir)//"Freq_Yrs_"//trim(IDX_nm(iidx))//"_All_Ens.nc" Do ilon=1,Nlon Do ilat=1,Nlat Do ipe=1,NPeSc X4D(ilon,ilat,ipe,1:Nyear)=XFry(iidx,ilon,ilat,1:Nyear,ipe) Enddo Enddo Enddo Call wr_ncdf_4d (Fname,xlon,xlat,xsce,xyear,nlon,nlat,NPeSc,Nyear,X4D,"fr_y", & t_unit, var_ln,"-", lon_nm, lat_nm, z_nm, t_nm, Miss) Enddo Deallocate (X4D) ! GeoRgVN: Allocate (X3D(Nreg+1,NPeSc,Ntime)) var_ln="Ensemble mean of Geographic Extent by sub-regions and VN (Z: BL,MF45,MF85,FF45,FF85)" lon_nm="reg"; lat_nm="sce"; t_nm="mon"; Do iidx=1,NDr_idx Fname=trim(Out_Dir)//"Geo_Reg_VN"//trim(IDX_nm(iidx))//"_All_Ens.nc" Do ireg=1,Nreg+1 Do ipe=1,NPeSc X3D(ireg,ipe,1:Ntime)=GeoRgVN(iidx,ireg,1:Ntime,ipe) Enddo Enddo Call wr_ncdf_3d (Fname,xreg,xsce,xtime,Nreg+1,NPeSc,Ntime, X3D, "geo",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Enddo Deallocate (X3D) End Subroutine Write_Char_Ens_to_NC !--------------------------------------------------------------------------------- Subroutine Write_Char_to_NC(Fname,Fn_geo,Fn_geo_pdf,Fn_geo_mon) Character (Len=*) :: Fname, Fn_geo,Fn_geo_pdf,Fn_geo_mon Character (Len=50) :: var_n, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=250) :: var_ln Real :: x_idx(NDr_idx), xmon(Ntime), xreg(Nreg), xmon1(Nmon), xreg1(Nreg+1) Real, Allocatable, Dimension (:,:,:) :: Xtmp Integer :: ireg, iidx, iyr lon_nm="lon"; lat_nm="lat"; t_nm="idx"; t_unit="years since 1987-01-01 " x_idx(1)=1; x_idx(2)=2; x_idx(3)=3 ! Drought duration: var_ln="Mean duration of drought (month). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="mon" Call wr_ncdf_3d (Fname,xlon,xlat,x_idx,nlon,nlat,NDr_idx, Dur_M, "dur_m",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum duration of drought (month). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, Dur_X,"dur_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought severity: var_ln="Mean severity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="-" Call wr_ncdf_3d_apd (Fname, Sev_M,"sev_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum severity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, Sev_X,"sev_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought severity (New definition): var_ln="Mean severity of drought (New: Sum of DrIDX minus Threshold). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="-" Call wr_ncdf_3d_apd (Fname, SEVm,"svm",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum severity of drought(New: Sum of DrIDX minus Threshold). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, SEVx,"svx",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought intensity: var_ln="Mean intensity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, Inten_M,"inten_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum intensity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, Inten_X,"inten_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought intensity (New definition): var_ln="Mean intensity of drought (New: Sum of DrIDX minus Threshold divide to Duration). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, INTEm,"itm",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum intensity of drought (New: Sum of DrIDX minus Threshold divide to Duration). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, INTEx,"itx",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought inter-arrival: var_ln="Mean inter-arrival of drought"; var_unit="mon" Call wr_ncdf_3d_apd (Fname, Inter_M,"inter_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Maximum inter-arrival of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d_apd (Fname, Inter_X,"inter_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Drought characteristics' frequency (Mean, q25, q75, IQR of Freq): var_ln="Maximum Annual drought frequency (mon/year)" var_unit="mon/year" Call wr_ncdf_3d_apd (Fname, Fr_Max,"fr_max",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Standard deviation of Annual drought frequency (mon/year)" Call wr_ncdf_3d_apd (Fname, Fr_Std,"fr_std",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !var_ln="First quartile (q25) of annual drought frequency (mon/year)" !Call wr_ncdf_3d_apd (Fname, q25,"fr_q25",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !var_ln="Third quartile (q75) of annual drought frequency (mon/year)" !Call wr_ncdf_3d_apd (Fname, q75,"fr_q75",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !var_ln="Interquartile (IQR) of annual drought frequency (mon/year)" !Call wr_ncdf_3d_apd (Fname, IQR,"fr_iqr",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Geographic Extent: Geo(:,:) ! (Ntime,NDr_idx) Do iyr=1, Ntime; xmon(iyr)=iyr; Enddo var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1987. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="%"; lon_nm="mon"; lat_nm="idx"; t_nm="reg" Allocate (Xtmp(Ntime,NDr_idx,Nreg)) Xtmp(1:Ntime,1:NDr_idx,1:Nreg)=GeoReg(1:Ntime,1:NDr_idx,1:Nreg) Do ireg=1, Nreg; xreg(ireg)=ireg; Enddo var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1987 for each subregions (1-7). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_3d (Fn_geo,xmon,x_idx,xreg,Ntime,NDr_idx,Nreg, Xtmp,"georeg",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1987 for entire country. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" Call wr_ncdf_2d_apd (Fn_geo,Geo,"geovn",Ntime,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, Miss) ! PDF of Geographic Extent: Deallocate (Xtmp); Allocate (Xtmp(NGeo_bin,NDr_idx,Nreg)) Do iidx=1,NDr_idx; Do ireg=1,Nreg; Xtmp(:,iidx,ireg)=Geo_PDF(iidx,ireg,:); Enddo; Enddo var_ln="PDF of geographic extent (%) calculated over entire period" var_unit="%"; lon_nm="bin"; lat_nm="idx"; t_nm="reg" Call wr_ncdf_3d (Fn_geo_pdf, XGeo_bin,x_idx,xreg, NGeo_bin,NDr_idx,Nreg, Xtmp,"geopdf",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) ! Monthly Geo Extent for sub-reg & VN: Do imon=1,Nmon; xmon1(imon)=imon; Enddo Do ireg=1,Nreg+1; xreg1(ireg)=ireg; Enddo var_ln="Monthly geographic extent (%) for sub-regs (1-7) and VN (8)" var_unit="%/mon"; lon_nm="mon"; lat_nm="idx"; t_nm="reg" Call wr_ncdf_3d (Fn_geo_mon, xmon1,x_idx,xreg1, Nmon,NDr_idx,Nreg+1,GeoRegMon,"geomon",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Deallocate (Xtmp) End Subroutine Write_Char_to_NC !--------------------------------------------------------------------------------- Subroutine Write_Fre_to_NC (Fname1, Fname2) Character (Len=*) :: Fname1, Fname2 Character (Len=50) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Real :: x_idx(NDr_idx), xyear(Nyear), xmon(Nmon) Integer :: iyr, iidx Real, Allocatable :: X(:,:,:,:) !Nlon,Nlat,NDr_idx,Nyear) ! Freq: lon_nm="lon"; lat_nm="lat"; t_nm="year"; z_nm="idx" t_unit="years since 1977-01-01 " x_idx(1)=1; x_idx(2)=2; x_idx(3)=3 Do iyr=1, Nyear; xyear(iyr)=Yr1+iyr; Enddo Allocate (X(Nlon,Nlat,NDr_idx,Nyear)) Do ilon=1, Nlon Do ilat=1, Nlat Do iidx=1, NDr_idx Do iyr=1, Nyear X(ilon,ilat,iidx,iyr) = Freq(ilon,ilat,iyr,iidx) Enddo Enddo Enddo Enddo var_ln="Number of drought month in the year. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="mon" print*, "Writing to file ..." Call wr_ncdf_4d (Fname1,xlon,xlat,x_idx,xyear,nlon,nlat,NDr_idx,Nyear, X, "freq", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) !FrMo: t_nm="mon"; t_unit="months since 1977-01-01 " Do iyr=1, Nmon; xmon(iyr)=iyr; Enddo ! months in the year Deallocate(X); Allocate (X(Nlon,Nlat,NDr_idx,Nmon)) Do ilon=1, Nlon Do ilat=1, Nlat Do iidx=1, NDr_idx X(ilon,ilat,iidx,:) = FrMo(ilon,ilat,:,iidx) Enddo Enddo Enddo var_ln="Monthly Drought frequency (%): Drought frequency in each month. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12" var_unit="%" Call wr_ncdf_4d (Fname2,xlon,xlat,x_idx,xmon,nlon,nlat,NDr_idx,Nmon, X, "freq", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) Deallocate(X) End Subroutine Write_Fre_to_NC !--------------------------------------------------------------------------------- Subroutine Calc_Chars Integer :: ilon, ilat, iidx, it, iyr, imon, imon1 Real :: X(Ntime), M1(Ntime), IMAX1, Spells(Ntime), N_Spell_1, I_Start(Ntime), I_End(Ntime) Integer :: N, i, k, N_Spell, Imax, ista, iend Real :: tmp, ntmp, Smean, Smax, In_mean, In_max, Amean, Se0, Sem, Sex, Inm, Inx Real :: XX(Nyear) M1=Miss; Spells=Miss; I_Start=Miss;I_End=Miss Do iidx=1,NDr_idx Do ilon=1, Nlon Do ilat=1, Nlat if (Mask_all(ilon,ilat) /= Miss) then X(:) = Events(ilon,ilat,:,iidx) Call NSpell_and_Leng (X, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss) print"(A,20I6)","iidx,Imax1,Maxval(Spells), N_Spell: ", iidx,int(Imax1),int(Maxval(Spells)), Int(N_Spell_1) N_Spell=Int(N_Spell_1); Imax=Int(Imax1) ! So dot va Do dai dot dai nhat tmp=0 Do i=1,N_Spell tmp=tmp+Spells(i) ! Tong do dai tat ca cac dot Enddo if (N_Spell>0) then Dur_M(ilon,ilat,iidx)=tmp/N_Spell_1; Dur_X(ilon,ilat,iidx)=Imax1 else Dur_M(ilon,ilat,iidx)=Miss; Dur_X(ilon,ilat,iidx)=Miss endif ! Do dai trung binh & Do dai dot dai nhat ! Severity and Intensity: Smean=0; Smax=999.; In_mean=0; In_max=999. Sem=0; Sex=999; Inm=0; Inx=999 ! New definitions of severity and intensity Do i=1, N_Spell ista=Int(I_Start(i)); iend=Int(I_End(i)) Amean=0; Se0=0 Do k=ista, iend ! Thoi diem bat dau/ket thuc dot thu i Amean = Amean + Dr_IDX(ilon,ilat,k,iidx) ! Seve dot thu i = Tong Dr_IDX trong dot thu i Se0 = Se0 + (Dr_IDX(ilon,ilat,k,iidx)-Qp(ilon,ilat,iidx)) ! Sum of (Dr_IDX-Qp) in a spell Enddo Smean = Smean+Amean ! Tong Dr_IDX trong tat ca cac dot (Tong Seve) Sem = Sem + Se0 ! Sum of (Dr_IDX-Qp) in all spells if (Smax>Amean) Smax = Amean ! Max Severity cua cac dot if (Sex >Sem ) Sex = Sem ! Max of Seve Amean = Amean/Spells(i) ! Intensity dot thu i Se0 = Se0 /Spells(i) ! Intensity dot thu i (New definition) In_mean = In_mean+Amean ! Tong Intensity cua cac dot Inm = Inm +Se0 ! Tong Intensity cua cac dot (New definition) if (In_max>Amean) In_max = Amean ! Max Intensity cua cac dot if (Inx >Se0 ) Inx = Se0 ! Max Intensity cua cac dot (New definition) Enddo if (N_Spell>0) then Sev_M(ilon,ilat,iidx) = Smean/N_Spell_1; Sev_X(ilon,ilat,iidx) = Smax Inten_M(ilon,ilat,iidx) = In_mean/N_Spell_1; Inten_X(ilon,ilat,iidx) = In_max ! (New definition): SEVm (ilon,ilat,iidx) = Sem/N_Spell_1; SEVx (ilon,ilat,iidx) = Sex INTEm(ilon,ilat,iidx) = Inm/N_Spell_1; INTEx(ilon,ilat,iidx) = Inx else Sev_M(ilon,ilat,iidx) = Miss; Sev_X(ilon,ilat,iidx) = Miss Inten_M(ilon,ilat,iidx) = Miss; Inten_X(ilon,ilat,iidx) = Miss ! (New definition): SEVm (ilon,ilat,iidx) = Miss; SEVx (ilon,ilat,iidx) = Miss INTEm(ilon,ilat,iidx) = Miss; INTEx(ilon,ilat,iidx) = Miss endif ! Seve & Inten ! Inter-arrival: if (N_Spell>3) then ! It nhat phai co 4 dot Amean=0; tmp=0 Do i=1, N_Spell-1 ntmp = I_Start(i+1)-I_Start(i) ! Khoang giua hai lan Onset lien tiep Amean = Amean+ntmp ! Tong cac Inter-Arrival if (ntmp>tmp) tmp = ntmp ! Max Inter-Arrival Enddo Inter_M(ilon,ilat,iidx) = Amean/(N_Spell_1-1); Inter_X(ilon,ilat,iidx) = tmp else Inter_M(ilon,ilat,iidx) = Miss; Inter_X(ilon,ilat,iidx) = Miss endif ! Freq: So events/nam (So thang han/nam) it=0 Do iyr=1, Nyear Freq(ilon,ilat,iyr,iidx) = 0 Do imon=1, Nmon it=it+1; Freq(ilon,ilat,iyr,iidx) = Freq(ilon,ilat,iyr,iidx) + Events(ilon,ilat,it,iidx) Enddo Enddo ! Chars of Freq: Mean, q25, q75, IQR !XX(1:Nyear) = Freq(ilon,ilat,1:Nyear,iidx) !X(1:Ntime) = Events(ilon,ilat,1:Ntime,iidx) !print*,xlon(ilon),xlat(ilat),Int(Sum(X,mask=X>=0)) !Mean(ilon,ilat,iidx) = Sumx(X, Ntime, Miss)!/Nyear !Sum(X,mask=X>=0)/Nyear ! Trung binh Mon/Year !print*,Mean(ilon,ilat,iidx) XX(1:Nyear) = Freq(ilon,ilat,1:Nyear,iidx) Fr_Max(ilon,ilat,iidx) = Maxval(XX) ! Nam co so thang han nhieu nhat Fr_Std(ilon,ilat,iidx) = Stdev (XX,Nyear,Miss) !q25(ilon,ilat,iidx) = Stdev (X,Ntime,Miss) ! !q75(ilon,ilat,iidx) = Maxval (XX(1:Nyear)) !if (q25(ilon,ilat,iidx)/=0) then ! IQR(ilon,ilat,iidx)=q25(ilon,ilat,iidx)/Mean(ilon,ilat,iidx) !else ! IQR(ilon,ilat,iidx)=Miss !endif ! FrMo: Tan suat xuat hien trong tung thang (%/mon) Do imon=1, Nmon FrMo(ilon,ilat,imon,iidx) = 0 Do it=1, Ntime imon1=mod(it,Nmon); if (imon1==0) imon1 = Nmon if (imon1==imon) then FrMo(ilon,ilat,imon,iidx) = FrMo(ilon,ilat,imon,iidx) + Events(ilon,ilat,it,iidx) endif Enddo FrMo(ilon,ilat,imon,iidx) = FrMo(ilon,ilat,imon,iidx)/Real(Nyear)*100.0 Enddo endif Enddo Enddo Enddo ! Geo Ext for entire country: Do iidx=1,NDr_idx Do it=1, Ntime Geo(it,iidx) = 0 Do ilon=1, Nlon Do ilat=1, Nlat if (Mask_all(ilon,ilat) /= Miss) then Geo(it,iidx) = Geo(it,iidx) + Events(ilon,ilat,it,iidx) endif Enddo Enddo Geo(it,iidx)=Geo(it,iidx)/Real(Nst)*100.0 Enddo !print*,Geo(:,iidx);pause Enddo ! Geo Ext for 7 subregions: GeoReg(:,:,:) ! (Ntime,NDr_idx,Nreg) Do iidx=1,NDr_idx Do it=1, Ntime Do ireg=1,Nreg GeoReg(it,iidx,ireg)=0; igrd=0 Do ilon=1, Nlon Do ilat=1, Nlat if (Mask_all(ilon,ilat) == ireg+1) then igrd=igrd+1 GeoReg(it,iidx,ireg) = GeoReg(it,iidx,ireg) + Events(ilon,ilat,it,iidx) endif Enddo Enddo if (igrd>0) then GeoReg(it,iidx,ireg)=GeoReg(it,iidx,ireg)/Real(igrd)*100.0 else GeoReg(it,iidx,ireg)=Miss endif Enddo Enddo Enddo ! Geo Ext Monthly for sub-reg and for entire country: Do iidx=1,NDr_idx Do ireg=1,Nreg ! 7 sub-reg Do imon=1,Nmon tmp=0; ntmp=0 Do it=1, Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon.AND.GeoReg(it,iidx,ireg)/=Miss) then ntmp=ntmp+1; tmp=tmp+GeoReg(it,iidx,ireg) endif Enddo GeoRegMon(imon,iidx,ireg)=tmp/ntmp Enddo Enddo Do imon=1,Nmon ! Entire VN tmp=0; ntmp=0 Do it=1, Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon.AND.Geo(it,iidx)/=Miss) then ntmp=ntmp+1; tmp=tmp+Geo(it,iidx) endif Enddo GeoRegMon(imon,iidx,Nreg+1)=tmp/ntmp Enddo Enddo ! PDF of Geo Ext: Do iidx=1,NDr_idx Do ireg=1,Nreg Do i=1,NGeo_bin ! Tan so/tan suat khoang tmp=0; ntmp=0 Do it=1,Ntime Smean=GeoReg(it,iidx,ireg) if (Smean/=Miss) then ntmp=ntmp+1 ! Dem so phan tu if (Smean>=Geo_bin(i-1).AND.SmeanNtime/2) then ! Dung luong du lon Geo_PDF(iidx,ireg,i) = tmp/ntmp*100.0 ! (%) else Geo_PDF(iidx,ireg,i) = Miss endif Enddo Enddo Enddo End Subroutine Calc_Chars !--------------------------------------------------------------------------------- Subroutine Calc_Per_Events (p) ! Return: Qp and Events Integer :: ilon, ilat, iidx, it Real :: p Real :: X(Ntime) Do iidx=1,NDr_idx Do ilon=1, Nlon Do ilat=1, Nlat if (Mask_all(ilon,ilat) /= Miss) then X(:) = Dr_IDX(ilon,ilat,:,iidx) Qp(ilon,ilat,iidx) = percentile (X,Ntime,p,Miss) Do it=1, Ntime if (Dr_IDX(ilon,ilat,it,iidx) < Qp(ilon,ilat,iidx)) then Events(ilon,ilat,it,iidx) = 1 else Events(ilon,ilat,it,iidx) = 0 endif Enddo !print "(30F6.2)",X(:);print* !print "(30F6.2)", Qp(ilon,ilat,iidx);print* !print "(30I6)", Int(Events(ilon,ilat,:,iidx)); print* !X(:)=Events(ilon,ilat,:,iidx) !Do it=1,Ntime; if (X(it)>0) print*,it," ",Dr_IDX(ilon,ilat,it,iidx),Qp(ilon,ilat,iidx),X(it); Enddo !print*,Int(Sum(X,mask=X>0)) !print*,Int(sum(Events(ilon,ilat,1:Ntime,iidx))); pause endif Enddo Enddo Enddo End Subroutine Calc_Per_Events !--------------------------------------------------------------------------------- Subroutine Read_Drough_IDX (Fname) Integer :: ist, ilon, ilat, it, Nt, imon, imon1, iyr Integer :: Nyear, Nmon, nlat1, nlat0, nlon, nlat Character (Len=*) :: Fname Real :: tmp, xmean(12), ic(12) Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) Real, Allocatable :: xlon(:), xlat(:) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "ori_pdsi") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "spei_12") Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "spi_12") Deallocate (xlon, xlat, xtime) !LandMaskFile="/work/projects/G09/SEACLID/Land_Mask_VN.nc" !Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"all") ! all=2..8=7 regs Nyear=20 ! 1986-2005 Nmon =12 Ntime = Nt - Nmon !Allocate (Mask_all(nlon,nlat)) !Mask_all=Mask Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime)) !X1(:,1:nlat,1:Ntime)=PDSI(:,Nlat0+1:Nlat1,Nmon+1:Nt) !X2(:,1:nlat,1:Ntime)=SPEI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) !X3(:,1:nlat,1:Ntime)=SPI12(:,Nlat0+1:Nlat1,Nmon+1:Nt) X1(:,:,1:Ntime)=PDSI (:,:,Nmon+1:Nt) X2(:,:,1:Ntime)=SPEI12(:,:,Nmon+1:Nt) X3(:,:,1:Ntime)=SPI12 (:,:,Nmon+1:Nt) Deallocate (PDSI,SPEI12,SPI12) Allocate (PDSI(nlon,nlat,Ntime),SPEI12(nlon,nlat,Ntime),SPI12(nlon,nlat,Ntime)) PDSI=X1; SPEI12=X2; SPI12=X3 Allocate (Dr_IDX(nlon,nlat,Ntime,NDr_idx)) Dr_IDX(:,:,:,:) = Miss !Dr_IDX(:,:,:,1)=PDSI(:,:,:); Dr_IDX(:,:,:,2)=SPI12(:,:,:); Dr_IDX(:,:,:,3)=SPEI12(:,:,:) Deallocate (X1, X2, X3) ist=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat)>0) then if (.not.Empty_Ar1(PDSI (ilon,ilat,1:Ntime), Ntime).AND. & .not.Empty_Ar1(SPI12 (ilon,ilat,1:Ntime), Ntime).AND. & .not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then ist=ist+1 Dr_IDX(ilon,ilat,:,1)=PDSI(ilon,ilat,:) Dr_IDX(ilon,ilat,:,2)=SPI12(ilon,ilat,:) Dr_IDX(ilon,ilat,:,3)=SPEI12(ilon,ilat,:) else Mask_all(ilon,ilat)=Miss endif endif Enddo Enddo Nst=ist ! Number of gridpoints Deallocate (PDSI,SPEI12,SPI12) print*,"Reading data ... Done !" End Subroutine Read_Drough_IDX !--------------------------------------------------------------------------------- Logical Function Empty_Ar1(X,N) Integer :: N,i Real :: X(N) Logical :: Empty Empty=.False. Do i=1,N if (X(i)==Miss.OR.IsNaN(X(i))) then Empty=.True.; Exit endif Enddo Empty_Ar1=Empty End Function Empty_Ar1 !--------------------------------------------------------------------------------- !--------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_Percentile.f90" Include "/home/tanpv/Mylibs/Sub_NSpell_and_Leng.f90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Stdev.f90" !Include "/home/tanpv/Mylibs/Sum.f90" END Program