Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Yr1=1976, Yr2=2005, Nmon=12, Nmem=7 Character (Len=100), parameter :: fmem(Nmem)=(/"CNRM","CSIRO","ECEA","GFDL","HadG","MPI","ENS"/) Character (Len=4) :: BC ! = "ORG" or "QM" !BL_CNRM_ORIG.nc BL_CSIRO_ORIG.nc BL_ECEA_ORIG.nc BL_Ens2_ORIG.nc !BL_Ens_ORIG.nc BL_GFDL_ORIG.nc BL_HadG_ORIG.nc BL_MPI_ORIG.nc !BL_CNRM_QM.nc BL_CSIRO_QM.nc BL_ECEA_QM.nc BL_Ens2_QM.nc BL_Ens_QM.nc !BL_GFDL_QM.nc BL_HadG_QM.nc BL_MPI_QM.nc ! 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 ! 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 :: NumEvt(:,:,:) ! (Nlon,Nlat,NDr_idx) - Number of Events 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 :: 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 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) Integer :: Nst, iidx, Nyear ! Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 ! Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & ! (/"pdsi","spi12","spe12"/) Integer, Parameter :: NDr_idx=2 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"spe12","spi12"/) Character (Len=100) :: Fname, Dir, Fch1, Fch2, Fch3, Fch4 !--------------------------------------------------------------------------------- !Dir="/work/users/hmo/remoclic/Phuong/PHD/BC_2023/Tan/RES_Dr_IDX/ORIG/"; BC="ORIG" Dir="/work/users/hmo/remoclic/Phuong/PHD/BC_2023/Tan/RES_Dr_IDX/QM/"; BC="QM" Do imem=1,Nmem Fname=trim(Dir)//"BL_"//trim(fmem(imem))//"_"//trim(BC)//".nc" ! "Dr_PDSI_SPEI_SPI.nc" Call Read_Drough_IDX (Fname) Nyear = Ntime/Nmon !print*,Dr_idx(:,:,20:30,3) 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)) 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 (NumEvt(Nlon,Nlat,NDr_idx)) !Allocate (Mean(Nlon,Nlat,NDr_idx), q25(Nlon,Nlat,NDr_idx), q75(Nlon,Nlat,NDr_idx), IQR(Nlon,Nlat,NDr_idx)) Allocate (Fr_max(Nlon,Nlat,NDr_idx), Fr_Std(Nlon,Nlat,NDr_idx)) ! Set all to Miss: Freq=Miss; FrMo=Miss; Geo=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 NumEvt=Miss Call Calc_Chars Fch1="NC/Freq_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Fch2="NC/Freq_mon_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Call Write_Fre_to_NC (Fch1,Fch2) !"OUT_CHAR/Freq_Ens2_QM_BL.nc","OUT_CHAR/Fre_mon_Ens2_QM_BL.nc") Fch3="NC/Chars_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Fch4="NC/Geo_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Call Write_Char_to_NC(Fch3,Fch4) ! "OUT_CHAR/Chars_Ens2_QM_BL.nc","OUT_CHAR/Geo_Ens2_QM_BL.nc") Deallocate (Mask_all,Dr_IDX,Qp,Events,Freq,FrMo,Geo,Dur_M,Dur_X) Deallocate (Sev_M,Sev_X,Inten_M,Inten_X,Inter_M,Inter_X,NumEvt) Deallocate (Fr_max,Fr_Std,Mask) Deallocate (SPEI12,SPI12) !,PDSI) Deallocate (xlon,xlat,xtime) Enddo CONTAINS !--------------------------------------------------------------------------------- Subroutine Write_Char_to_NC(Fname,Fn_geo) Character (Len=*) :: Fname, Fn_geo Character (Len=150) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Real :: x_idx(NDr_idx), xmon(Ntime) lon_nm="lon"; lat_nm="lat"; t_nm="idx"; t_unit="years since 1977-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:SPEI12, 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:SPEI12, 2:SPI12" !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:SPEI12, 2:SPI12" !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:SPEI12, 2:SPI12" !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 intensity: var_ln="Mean intensity of drought. IDXs is 1:SPEI12, 2:SPI12" !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:SPEI12, 2:SPI12" !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 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:SPEI12, 2:SPI12" !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) ! Number of drought events: var_ln="Total number of drought events. IDXs is 1:SPEI12, 2:SPI12" !1:PDSI, 2:SPI12, 3:SPEI12" var_unit="event" Call wr_ncdf_3d_apd (Fname,NumEvt,"evt",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: Do iyr=1, Ntime; xmon(iyr)=iyr; Enddo var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1977. IDXs is 1:SPEI12, 2:SPI12" !1:PDSI, 2:SPI12, 3:SPEI12" var_unit="%"; lon_nm="mon"; lat_nm="idx" Call wr_ncdf_2d (Fn_geo,xmon,x_idx,Ntime,NDr_idx, Geo,"geo", var_ln, var_unit, lon_nm, lat_nm, Miss) 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:SPEI12, 2:SPI12" !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:SPEI12, 2:SPI12" !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 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 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 NumEvt(ilon,ilat,iidx)=N_Spell_1 ! So dot trong toan chuoi 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=0.; In_mean=0; In_max=0. Do i=1, N_Spell ista=Int(I_Start(i)); iend=Int(I_End(i)) Amean=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 Amean = Amean + ABS(Dr_IDX(ilon,ilat,k,iidx)) ! Lay tri tuyet doi Enddo Smean = Smean+Amean ! Tong Dr_IDX trong tat ca cac dot (Tong Seve) if (Smax0) 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 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 endif ! Seve & Inten ! Inter-arrival: if (N_Spell>1) then 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: 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 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) Qp(ilon,ilat,iidx) = -1.0 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 Character (Len=*) :: Fname Real :: tmp, xmean(12), ic(12) Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:) ! 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) Fname="/work/users/hmo/tanpv/CORDEX_SEACLID/Land_mask_sea.nc" Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land") !Phuong chinh sua 1976-1995 Nyear=30 ! 1976-2005 ! Nyear =20 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 (SPEI12,SPI12) !PDSI) Allocate (SPEI12(nlon,nlat,Ntime),SPI12(nlon,nlat,Ntime)) !,PDSI(nlon,nlat,Ntime)) ! PDSI=X1; SPEI12=X2; SPI12=X3 !Deallocate (X1,X2,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. & if (.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,:,1)=SPEI12(ilon,ilat,:) else Mask_all(ilon,ilat)=Miss endif endif Enddo Enddo Nst=ist ! Number of gridpoints !Deallocate (PDSI,SPEI12,SPI12) !Deallocate (X1,X2,X3) 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 "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/sub_rd_ncdf_Tan.F90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/sub_wr_ncdf_Tan.F90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sub_Percentile.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sub_NSpell_and_Leng.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Ave.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Stdev.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sum.f90" END Program