Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6 Integer, Parameter :: Ncase = 3 ! mod,cli,ann Integer, Parameter :: Ncase_1 = 2 ! obs, ens Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"obs","ens"/) Character (Len=3), Parameter :: Cases(Ncase)=(/"mod","cli","ann"/) !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_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_Out(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel,Nlen !------------------------------------------------------------------------------ ! 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, Yr_Ope, Mon_Ope, Nmodel, NLen, Nidx Namelist /Model_name/ Model_nm, IDXname Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_Dr_IDXs.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 = (Yr_Ope - Yr1 + 1) Call Read_Regs print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime Call Reinit Call Read_all_data Allocate (SPEI_Out(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel,Nlen)); SPEI_Out=Miss Do idx=1,Nidx Do ilen=1,Nlen Call Calc_IDXs (trim(IDXname(idx)),ilen) Call Write_Pr_Temp_DrIDXs (trim(IDXname(idx)),ilen) Enddo Enddo !Call Write_Pr_Temp_DrIDXs ("spei") CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_IDXs (IDX_nm,ilen) ! IDX_nm="spi", "spei" Character (Len=*) :: IDX_nm Integer :: ilen,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, latx !Real :: xspi(Ntime) ! 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,Nmodel)) SPEI_All=Miss Do imod=1,Nmodel Do ilead=1,Nlead Nt=(Yr_Ope - Yr1)*Nmon + Mon_Ope + ilead ! Do dai chuoi SL thuc te (khong tinh Miss) if (Allocated (Rain)) Deallocate (Rain) if (Allocated (T2m)) Deallocate (T2m) Allocate (Rain(Nt),T2m(Nt)) Do imem=1,Nmem !Do ilen=1,Nlen Do icas=1,Ncase Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then use_baseline=.false.; xspi=Miss; latx=Lat(ilat) !Rain(1:Ntime)=Pr_All(ilon,ilat,imem,ilead,1:Ntime,icas,imod) !T2m (1:Ntime)=T2_All(ilon,ilat,imem,ilead,1:Ntime,icas,imod) 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,Ntime,Miss), Length_series (T2m,Ntime,Miss)) 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) !print*,Length_series (xspi,Nt1,Miss),xspi(Nt1),Rain(Nt1),T2m(Nt1) endif SPEI_Out(ilon,ilat,imem,ilead,icas,imod,ilen)=SPEI_All(ilon,ilat,imem,ilead,Nt,icas,imod) !print"(7I5,F8.4)",Len_SPI(ilen),imod,ilead,imem,icas,Nt1,Nt,SPEI_Out(ilon,ilat,imem,ilead,icas,imod,ilen) !print"(A,10F8.4)","SPEI_All: ",SPEI_All(ilon,ilat,imem,ilead,Nt1-4:Nt1,icas,imod) endif endif Enddo Enddo Enddo !Enddo Enddo Enddo !write (Fname,"(A,I1,A)") trim(Out_Dir)//trim(IDX_nm)//"_ens_lead_",ilead,".nc" !Call Write_Dr_IDXs_to_file (Fname,ilead) !,Model_nm(imod)) Enddo End Subroutine Calc_IDXs !------------------------------------------------------------------------------ Subroutine Write_Pr_Temp_DrIDXs (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: 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(:,:,:,:),xcase(:) ! Write out for testing: SPEI_Out(ilon,ilat,imem,ilead,icas,imod,ilen) ! 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_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" Nz=Nmem*Ncase*Nmodel Allocate (X4d(Nlon,Nlat,Nz,Nlead),xcase(Nz)) 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 icas=1,Nz; xcase(icas)=icas; Enddo Where (IsNaN(SPEI_Out)) SPEI_Out=Miss i=0 Do imod=1,Nmodel Do imem=1,Nmem Do icas=1,Ncase i=i+1 Do ilead=1,Nlead X4d(:,:,i,ilead)=SPEI_Out(:,:,imem,ilead,icas,imod,ilen) Enddo Enddo Enddo Enddo write(Var_ln,"(A,I2.2,A)")"Values of "//trim(IDX_nm)//" drought index for timescale ",Len_SPI(ilen), & " month. The order is mod, cli, ann for each member (1..7) then for each model (1..2)" 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,xcase) 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 !it2=it1+Ntime-1 ! Tinh den 12/2021 ! Ntime=(Yr2-Yr1)*Nmon - bo nam dau tien 1982 Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)) Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)) 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) Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))) !print*,"Nt,Ntime=",Nt,Ntime,Nt-Nmon; pause !if (imod==1.AND.icas==1.AND.ilead==1) then ! Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel),time(Ntime)); ! T2_All=Miss; time(1:Ntime)=xtime(1:Ntime) ! Luu lai de write out !endif 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" Do icas=1,Ncase 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 !if (imod==1.AND.icas==1.AND.ilead==1) then; Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)); Pr_All=Miss; endif 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 !if (Allocated(xtime)) Deallocate (xtime) !Allocate (xtime(Ntime)) !pause;print*,time !xtime(:)=time(:) ! restore xtime !Deallocate (time) 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