Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI Use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2, Nreg = 7, Nidx = 2 Integer, Parameter :: NMonth = 12, Nmon = NMonth, NLen=3 ! 6,12 Character (Len=5),Parameter :: IDXname(Nidx) = (/"spi","spei"/) Integer :: Yr1, Yr2 Integer, Parameter :: Len_SPI(NLen)=(/3,6,12/) ! Integer :: Len_SPI(NLen) !!! Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI Real, Allocatable :: DatPr(:,:), DatTx(:,:), DatTm(:,:), DatT2m(:,:) ! For calculating SPI, SPEI Real, Allocatable :: Rain(:), T2m(:) Real, Allocatable :: RainReg(:,:), TxReg(:,:), TmReg(:,:), T2mReg(:,:) ! Regional mean Character (Len=20), Allocatable :: Stname(:) Real, Allocatable :: Lon(:), Lat(:), Reg(:) Real, Allocatable :: beta_spi(:), gamm(:), pzero(:) Real :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2) Real, Allocatable :: SPI(:,:,:), SPEI_th(:,:,:), SPEI_ha(:,:,:) !------------------------------------------------------------------------------ ! Statistics: Real, Allocatable, Dimension (:,:) :: Pr_Stat, Tx_Stat, Tm_Stat, T2m_Stat ! (Nst,NStat) - Nstat=4: Max,Min,Mean,Std !------------------------------------------------------------------------------ !Real, Allocatable :: SPI_reg(:,:,:), SPEI_reg_th(:,:,:), SPEI_reg_ha(:,:,:) ! Notes: ! 1) Dat_Mon(ivar,ist,iyr,imw), where ivar=1:Temp (F), ivar=2:Rainfall (inch); ist=1..Nst; iyr=YBas1..YBas2; imw=1..Nmon (or Nweek) ! Dat_Mon used to calculating PDSI ! 2) Temp_Inp (Temp in C deg), Rain_Inp (Rainfall in mm/month) (1:Nlon1,1:Nlat1,1:Ntime), Ntime=(YBas2-YBas1+1)*Nmon ! Temp_Inp and Rain_Inp to be used to write out as input data ! 3) Prec (mm/month), Temp (C) (1:Ntime,1:Nst) ! 4) DatPr(1:Ntime,Nst), DatTx(1:Ntime,Nst), DatTm(1:Ntime,Nst), DatT2m(1:Ntime,Nst) (1:Ntime,1:Nst): Used to calculating SPI, SPEI !------------------------------------------------------------------------------ Character (Len=255) :: Fname, R_Fname, T_Fname, BL_Dir, PDSI_Dir, PDSI_WDir, Data_dir, MaskFname, AWC_Fname, Out_Dir Character (Len=3) :: TVar, RVar Character (Len=50) :: Pr_nm, T_nm Logical :: status, Empty Integer :: YBas1, YBas2, NYear, maxyrs Integer :: Nst, Nlon, Nlat, IO, Ntime, Nt, it Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6 Real, allocatable, Dimension (:) :: X1d, Xt, xspi !------------------------------------------------------------------------------ ! Data input: Fname="/home/tanpv/OBS_Stn/Monthly/NC/R_Tx_Tm_T2m_1980-2020_mon_3D_VN.nc" !Fname="/home/tanpv/OBS_Stn/Monthly/NC/R_T2m_1980-2018_mon_Tram_KT_61_18_3D.nc" Call Read_T_R (Fname) Fname="/home/tanpv/Progs/Dr_4Books/Res/Pr_T_Statistics" Call Calc_Stat_Pr_T (Fname) Allocate (Rain(Ntime),T2m(Ntime)) Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth)) Allocate (SPI(Nst,Ntime,NLen), SPEI_th(Nst,Ntime,NLen)) !, SPEI_ha(Nst,NLen,Ntime)) Do iidx=1,Nidx print*,"Calculating ",IDXname(iidx) Call Calc_IDXs (IDXname(iidx)) Enddo ! Output: Fname="/home/tanpv/Progs/Dr_4Books/NC/spi_spei_stn_1980-2018.nc" Call Write_to_NC (Fname) !------------------------------------------------------------------------------ ! Calc regional mean: ! Deallocate (SPI,SPEI_th,SPEI_ha,DatPr, DatTx, DatTm, DatT2m) ! Nst=Nreg+1 ! Allocate (SPI(Nst,Ntime,NLen), SPEI_th(Nst,Ntime,NLen), SPEI_ha(Nst,Ntime,NLen)) ! Allocate (DatPr(Ntime,Nst),DatTx(Ntime,Nst),DatTm(Ntime,Nst),DatT2m(Ntime,Nst)) ! DatPr(:,:)=RainReg(:,:); DatTx(:,:)=TxReg(:,:); DatTm(:,:)=TmReg(:,:); DatT2m(:,:)=T2mReg(:,:) ! Calculating SPEI: ! SPI(:,:,:) = amssng; SPEI_th(:,:,:) = amssng ! Call Cal_SPEI_thornthwaite ! SPEI_th(:,:,:) = SPI(:,:,:) ! SPI(:,:,:) = amssng; SPEI_ha(:,:,:) = amssng ! Call Cal_SPEI_Hargreaves ! SPEI_ha(:,:,:) = SPI(:,:,:) ! Calculating SPI: ! SPI(:,:,:) = amssng ! Call Cal_SPI ! Fname="Results/spi_spei_Reg_1980-2020_4timescales.nc" ! Call Write_to_NC (Fname) CONTAINS !------------------------------------------------------------------------------ Subroutine Write_to_NC (Fname) Character (Len=*) :: Fname Real, Allocatable :: Stn(:), time(:), LenDr(:), X2D(:,:) Integer :: it, Nm, ist, imon, iyr, Nstn Character (Len=255) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Allocate (Stn(Nst),time(Ntime),LenDr(NLen)) Do ist=1,Nst; Stn(ist)=ist; Enddo Do it=1,Ntime; time(it)=it; Enddo Do it=1,NLen; LenDr(it)=Len_SPI(it); Enddo ! SPI t_unit="months since 1980-01-01 "; lon_nm="stn"; lat_nm="mon"; t_nm="len"; var_unit="-" var_ln="SPI for period 1980-2020 with drought timescale=3,6,12 months" Call wr_ncdf_3d (Fname,Stn,time,LenDr,Nst,Ntime,Nlen,SPI,"spi",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! SPEI var_ln="SPEI for period 1980-2020 with drought timescale=3,6,12 months, PET-thornthwaite" Call wr_ncdf_3d_apd (Fname,SPEI_th,"spei",Nst,Ntime,Nlen, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Precipitation Allocate (X2D(Nst,Ntime)) Do ist=1,Nst; X2D(ist,:)=DatPr(:,ist); Enddo var_ln="Monthly rainfall of all stations"; var_unit="mm/month" Call wr_ncdf_2d_apd (Fname,X2D,"pr",Nst,Ntime, var_ln, var_unit, lon_nm, lat_nm, Miss) ! Temperature Do ist=1,Nst; X2D(ist,:)=DatT2m(:,ist); Enddo var_ln="Monthly mean temperature of all stations"; var_unit="C" Call wr_ncdf_2d_apd (Fname,X2D,"t2m",Nst,Ntime, var_ln, var_unit, lon_nm, lat_nm, Miss) Do ist=1,Nst; X2D(ist,:)=DatTx(:,ist); Enddo var_ln="Monthly maximum temperature of all stations"; var_unit="C" Call wr_ncdf_2d_apd (Fname,X2D,"tx",Nst,Ntime, var_ln, var_unit, lon_nm, lat_nm, Miss) Do ist=1,Nst; X2D(ist,:)=DatTm(:,ist); Enddo var_ln="Monthly minimum temperature of all stations"; var_unit="C" Call wr_ncdf_2d_apd (Fname,X2D,"tm",Nst,Ntime, var_ln, var_unit, lon_nm, lat_nm, Miss) ! Lon, Lat, Reg and Station name var_ln="Longitude in East degree"; var_unit="deg-E" Call wr_ncdf_1d_apd (Fname,Lon,"lon", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Latitude in North degree"; var_unit="deg-N" Call wr_ncdf_1d_apd (Fname,Lat,"lat", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Climatic sub-regions"; var_unit="-" Call wr_ncdf_1d_apd (Fname,Reg,"reg", Nst, var_ln, var_unit, lon_nm, Miss) var_ln="Name of Stations"; var_unit="st" Call wr_ncdf_1d_apd_char (Fname,Stname,"stname", Nst, var_ln, var_unit, lon_nm, Miss) End Subroutine Write_to_NC !------------------------------------------------------------------------------ Subroutine Calc_IDXs (IDX_nm) ! 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, Allocatable :: xspi(:) Logical :: use_baseline if (Allocated (xspi)) Deallocate (xspi) Allocate (xspi(Ntime)) Do ist=1,Nst latx=Lat(ilat) Do ilen=1,Nlen use_baseline=.false.; xspi=Miss; Rain(1:Ntime)=DatPr (1:Ntime,ist) T2m (1:Ntime)=DatT2m(1:Ntime,ist) if (IDX_nm=="spi") then Call Cal_SPI_One_Stn (Rain,Ntime,Len_SPI(ilen),alpha,beta_spi,gamm,pzero,xspi,use_baseline,Miss) SPI(ist,1:Ntime,ilen)=xspi(1:Ntime) else if (IDX_nm=="spei") then Call Calc_SPEI_One_Stn (T2m,Rain,latx,Ntime,Nmon,Len_SPI(ilen),logLogisticParams,use_baseline,xspi,Miss) SPEI_th(ist,1:Ntime,ilen)=xspi(1:Ntime) endif Enddo Enddo End Subroutine Calc_IDXs !------------------------------------------------------------------------------ Subroutine Calc_Stat_Pr_T (Fname) ! Real, Allocatable, Dimension (:,:) :: Pr_Stat, Tx_Stat, Tm_Stat, T2m_Stat ! (Nst,NStat) - Nstat=4: Max,Min,Mean,Std ! Allocate (DatPr(Ntime,Nst),DatTx(Ntime,Nst),DatTm(Ntime,Nst),DatT2m(Ntime,Nst)) Character (Len=*) :: Fname Integer :: ist,i,j,k Integer, parameter :: Nstat=4 Character (Len=5) :: Statname(Nstat)=(/"Max","Min","Mean","Stdev"/) Real :: X(Ntime) Integer, Parameter :: Nvar=4 Character (Len=4) :: vname(Nvar)=(/"pr","tx","tm","t2m"/) Real, Allocatable :: Y(:,:,:) Allocate (Pr_Stat(Nst,NStat),Tx_Stat(Nst,NStat),Tm_Stat(Nst,NStat),T2m_Stat(Nst,NStat)) Do ist=1,Nst X(1:Ntime) = DatPr(1:Ntime,ist) Pr_Stat(ist,1)=Xmax(X,Ntime,Miss); Pr_Stat(ist,2)=Xmin(X,Ntime,Miss); Pr_Stat(ist,3)=Ave(X,Ntime,Miss); Pr_Stat(ist,4)=Stdev(X,Ntime,Miss) X(1:Ntime) = DatTx(1:Ntime,ist) Tx_Stat(ist,1)=Xmax(X,Ntime,Miss); Tx_Stat(ist,2)=Xmin(X,Ntime,Miss); Tx_Stat(ist,3)=Ave(X,Ntime,Miss); Tx_Stat(ist,4)=Stdev(X,Ntime,Miss) X(1:Ntime) = DatTm(1:Ntime,ist) Tm_Stat(ist,1)=Xmax(X,Ntime,Miss); Tm_Stat(ist,2)=Xmin(X,Ntime,Miss); Tm_Stat(ist,3)=Ave(X,Ntime,Miss); Tm_Stat(ist,4)=Stdev(X,Ntime,Miss) X(1:Ntime) = DatT2m(1:Ntime,ist) T2m_Stat(ist,1)=Xmax(X,Ntime,Miss); T2m_Stat(ist,2)=Xmin(X,Ntime,Miss); T2m_Stat(ist,3)=Ave(X,Ntime,Miss); T2m_Stat(ist,4)=Stdev(X,Ntime,Miss) Enddo Open (3,file=trim(Fname)//".txt", status="Unknown") Write (3,*) "Climatology of Rainfall, Maximun, Minimum and Mean Temperature" write (3,"(35x,4(A24,5x))") "Precipitation","Maximum Temp.", "Minimum Temp.", "Mean Temp." write (3,"(A4,A16,A8,A7,A3,4(4A7,A1))") "No ","Stn mane ","Lon ","Lat ","Reg",("Max","Min","Mean","Std","|",j=1,Nstat) Do ist=1,Nst write (3,"(I3,1x,A16,F8.3,F7.3,F3.0,4(4F7.1,1x))") ist,stname(ist),Lon(ist),Lat(ist),Reg(ist), & (Pr_Stat(ist,j),j=1,Nstat),(Tx_Stat(ist,j),j=1,Nstat),(Tm_Stat(ist,j),j=1,Nstat), (T2m_Stat(ist,j),j=1,Nstat) Enddo Close (1) ! To GrADS Allocate (Y(Nst,Nvar,Nstat)) Y(:,1,:) = Pr_Stat(:,:); Y(:,2,:) = Tx_Stat(:,:); Y(:,3,:) = Tm_Stat(:,:); Y(:,4,:) = T2m_Stat(:,:); Call Write_GrADS_Stn_3D (vname,Y,Nst,Nvar,Nstat,Lon,Lat,Stname,Fname) End Subroutine Calc_Stat_Pr_T !------------------------------------------------------------------------------ Subroutine Read_T_R (Fname) Character (Len=*) :: Fname Integer, parameter :: Nst_sea=9 Real, Allocatable :: X3D1(:,:,:),X3D2(:,:,:),X3D3(:,:,:), Y3D(:,:,:), Stn(:), Year(:), Mon(:) Real, Allocatable :: Lon1(:), Lat1(:), Reg1(:) Character (Len=20), Allocatable :: Stnm(:) Character (Len=20), parameter :: IslandStn(Nst_sea)=(/"COTO","BACHLONGVI","HONDAU", & "CONCO","LYSON","PHUQUY","CONDAO","TRUONGSA","PHUQUOC"/) Integer :: ireg, i, it, Nm, ist, ist1, ist2, imon, iyr, Nstn Real :: Nc1, Nc2, Nc3, Nc4, tmp1, tmp2, tmp3, tmp4 Call rd_ncdf_3d (Fname,Stn,Year,Mon,Nstn,Nyear,Nm, Y3D,"R") Deallocate (Stn,Year,Mon) Call rd_ncdf_3d (Fname,Stn,Year,Mon,Nstn,Nyear,Nm, X3D1,"Tx") Deallocate (Stn,Year,Mon) Call rd_ncdf_3d (Fname,Stn,Year,Mon,Nstn,Nyear,Nm, X3D2,"Tm") Deallocate (Stn,Year,Mon) Call rd_ncdf_3d (Fname,Stn,Year,Mon,Nstn,Nyear,Nm, X3D3,"T2m") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nstn, Lon1, "Lon") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nstn, Lat1, "Lat") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nstn, Reg1, "Reg") Deallocate (Stn) Call rd_ncdf_1d_Char (Fname,Stn,Nstn,Stnm,"Stname") !Yr1 = Int(Year(1)); Yr2 = Int(Year(Nyear)) ! 1980-2020 Yr1 = 1980; Yr2 = 2018 Ntime = Nyear*NMonth; maxyrs = Nyear Nst = Nstn - Nst_sea ! So tram tren dat lien Allocate (DatPr(Ntime,Nst),DatTx(Ntime,Nst),DatTm(Ntime,Nst),DatT2m(Ntime,Nst)) Allocate (Lon(Nst),Lat(Nst),Reg(Nst),Stname(Nst)) ist1=0 Loop1: Do ist=1,Nstn if (trim(Stnm(ist))=="COTO".OR.trim(Stnm(ist))=="BACHLONGVI".OR.trim(Stnm(ist))=="HONDAU".OR. & trim(Stnm(ist))=="CONCO".OR.trim(Stnm(ist))=="LYSON".OR.trim(Stnm(ist))=="PHUQUY".OR. & trim(Stnm(ist))=="CONDAO".OR.trim(Stnm(ist))=="TRUONGSA".OR.trim(Stnm(ist))=="PHUQUOC") Cycle Loop1 it=0; ist1=ist1+1 Do iyr=1,Nyear Do imon=1,Nm it=it+1 DatPr (it,ist1) = Y3D (ist,iyr,imon) DatTx (it,ist1) = X3D1(ist,iyr,imon) DatTm (it,ist1) = X3D2(ist,iyr,imon) DatT2m(it,ist1) = X3D3(ist,iyr,imon) Enddo Enddo Lon(ist1)=Lon1(ist); Lat(ist1)=Lat1(ist) Reg(ist1)=Reg1(ist); Stname(ist1)=Stnm(ist) Enddo Loop1 Where (IsNaN(DatPr)) DatPr = Miss Where (IsNaN(DatTx)) DatTx = Miss Where (IsNaN(DatTm)) DatTm = Miss Where (IsNaN(DatT2m)) DatT2m = Miss Where (DatTx>40.) DatTx = Miss Where (DatTm>40.) DatTm = Miss Where (DatT2m>40.) DatT2m = Miss ! Regional mean: Allocate (RainReg(Ntime,Nreg+1), TxReg(Ntime,Nreg+1), TmReg(Ntime,Nreg+1), T2mReg(Ntime,Nreg+1)) Do it=1,Ntime Do ireg=1,Nreg Nc1=0; Nc2=0; Nc3=0; Nc4=0; tmp1=0; tmp2=0; tmp3=0; tmp4=0 Do ist=1, Nst if (Reg(ist)==ireg) then if (DatPr(it,ist)/=Miss) then Nc1=Nc1+1; tmp1=tmp1+DatPr(it,ist) endif if (DatTx(it,ist)/=Miss) then Nc2=Nc2+1; tmp2=tmp2+DatTx(it,ist) endif if (DatTm(it,ist)/=Miss) then Nc3=Nc3+1; tmp3=tmp3+DatTm(it,ist) endif if (DatT2m(it,ist)/=Miss) then Nc4=Nc4+1; tmp4=tmp4+DatT2m(it,ist) endif endif Enddo RainReg(it,ireg) = tmp1/Nc1 TxReg(it,ireg) = tmp2/Nc2; TmReg(it,ireg) = tmp3/Nc3; T2mReg(it,ireg) = tmp4/Nc4 Enddo Nc1=0; Nc2=0; Nc3=0; Nc4=0; tmp1=0; tmp2=0; tmp3=0; tmp4=0 Do ist=1, Nst ! Entire country if (DatPr(it,ist)/=Miss) then Nc1=Nc1+1; tmp1=tmp1+DatPr(it,ist) endif if (DatTx(it,ist)/=Miss) then Nc2=Nc2+1; tmp2=tmp2+DatTx(it,ist) endif if (DatTm(it,ist)/=Miss) then Nc3=Nc3+1; tmp3=tmp3+DatTm(it,ist) endif if (DatT2m(it,ist)/=Miss) then Nc4=Nc4+1; tmp4=tmp4+DatT2m(it,ist) endif Enddo RainReg(it,Nreg+1) = tmp1/Nc1 TxReg (it,Nreg+1) = tmp2/Nc2 TmReg (it,Nreg+1) = tmp3/Nc3 T2mReg (it,Nreg+1) = tmp4/Nc4 !print"(8F10.1)",RainReg(it,:) !print"(8F10.1)",TempReg(it,:) Enddo; ! pause Deallocate (Stn,Year,Mon,X3D1,X3D2,X3D3,Y3D,Lon1,Lat1,Reg1,Stnm) End Subroutine Read_T_R !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" !Include "/home/tanpv/Mylibs/Sub_RW_Station_nc.f90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_SPI_One_Stn_2023.f90" ! Chuong trinh tinh chi so han SPI Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_2023.f90" ! Chuong trinh tinh chi so han SPEI Include "/home/tanpv/Mylibs/Xmax.f90" Include "/home/tanpv/Mylibs/Xmin.f90" Include "/home/tanpv/Mylibs/Stdev.f90" Include "/home/tanpv/Mylibs/Ave.f90" Include "Sub_GrADS.f90" !Include "/home/tanpv/Mylibs/Percentile.f90" !-------------------------------------------------------------------- END ! Program