Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI Use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2, Nreg = 7 Integer, Parameter :: NMonth = 12, NLen=4 ! 3,6,12 Integer :: Yr1, Yr2 Integer, Parameter :: Len_SPI(NLen)=(/1,3,6,12/) ! Integer :: Len_SPI(NLen) !!! Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI Real, Allocatable :: DatPr(:,:), DatTx(:,:), DatTm(:,:), DatT2m(:,:), Rain(:) ! For calculating SPI, SPEI 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, Allocatable :: SPI(:,:,:), SPEI_th(:,:,:), SPEI_ha(:,:,:) !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 (:) :: Xt !------------------------------------------------------------------------------ Fname="/home/tanpv/OBS_Stn/Monthly/NC/R_Tx_Tm_T2m_1980-2020_mon_3D_VN.nc" Call Read_T_R (Fname) ! Assign drought timescales: !Do ilen=1,NLen ! Len_SPI(ilen)=ilen !Enddo Allocate (Rain(Ntime)) Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth)) Allocate (SPI(Nst,Ntime,NLen), SPEI_th(Nst,Ntime,NLen), SPEI_ha(Nst,Ntime,NLen)) ! Calculating SPEI: !------------------------------------------------------------------------------ SPI(:,:,:) = amssng; SPEI_th(:,:,:) = amssng ! Call Cal_SPEI ! <<<<<<<< Should check again if SPE is SPI or SPEI ???? 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_stn_1980-2020_4timescales.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=1,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=1,3,6,12 months, PET-thornthwaite" Call wr_ncdf_3d_apd (Fname,SPEI_th,"spei_th",Nst,Ntime,NLen, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="SPEI for period 1980-2020 with drought timescale=1,3,6,12 months, PET-Hargreaves" Call wr_ncdf_3d_apd (Fname,SPEI_ha,"spei_ha",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 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)) 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 ! 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_2022.f90" ! Chuong trinh tinh chi so han SPI Include "/home/tanpv/Mylibs/Sub_SPEI_2022.f90" ! Chuong trinh tinh chi so han SPEI !Include "/home/tanpv/Mylibs/Percentile.f90" !-------------------------------------------------------------------- END ! Program