Use NetCDF Real, Parameter :: Miss=-99.0, amssng = Miss Integer, Parameter :: Nreg=7, Nmon=12, IMiss=-99, NLead=5, NMonth = Nmon Integer, Parameter :: Yr1=1961, Yr2=2018, Ntime=(Yr2-Yr1+1)*Nmon Integer, Parameter :: NYear=(Yr2-Yr1+1), maxyrs=NYear ! Need for Sub_SPI Integer, Parameter :: NLen=4 Integer, Parameter :: Len_SPI(NLen)=(/1,3,6,12/) Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255) :: Dir Character (Len=10) :: Model_Name Real, Allocatable :: lon(:), lat(:), Regs(:,:) !Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatPr(:,:), DatT2m(:,:), DatTx(:,:), DatTm(:,:), Rain(:) ! For calculating SPI, SPEI Real, Allocatable :: SPI(:,:,:) ! (Nst, Ntime, NLen) Real, Allocatable :: SPI_O(:,:,:), SPEI_OT(:,:,:), SPEI_OH(:,:,:) ! (Nst, Ntime, NLen) Real, Allocatable :: SPI_O1(:,:,:), SPEI_OT1(:,:,:) Real, Allocatable :: SPI_F(:,:,:,:), SPEI_FT(:,:,:,:), SPEI_FH(:,:,:,:) ! (Nst, Ntime, NLen, NLead) Real :: beta_spi(NMonth), gamm(NMonth), pzero(NMonth) Real, Allocatable :: beta_spi_store(:,:,:), gamm_spi_store(:,:,:), pzero_spi_store(:,:,:) ! (Nst,Nmon,Nlen) Real :: LogParams(0:NMonth,0:2) Real, Allocatable :: SPEI_Pars(:,:,:,:) ! (Nst,0:NMonth,0:2,NLen) !---------------------------------------------------------------------------------------------- ! Input data: Real, Allocatable, dimension(:,:,:,:) :: Pr, T2m, Tx, Tm ! (Nlon,Nlat,Nlead,Ntime) Real, Allocatable, dimension(:,:,:) :: Pr_O, T2m_O, Tx_O, Tm_O ! (Nlon,Nlat,Ntime) Real, Allocatable :: xlon(:), xlat(:), xtime(:), Reg(:) Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) Integer :: Nlon, Nlat, nlo, nla, iyr, imon, it, ist !---------------------------------------------------------------------------------------------- ! Model Climatology: Real, Allocatable, dimension(:,:,:,:) :: Pr_M_cli, T2m_M_cli, Tx_M_cli, Tm_M_cli ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:) :: Pr_O_cli, T2m_O_cli, Tx_O_cli, Tm_O_cli ! (Nlon,Nlat,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Pr_BS, T2m_BS, Tx_BS, Tm_BS ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- ! Differences between Model and Obs: Real, Allocatable, dimension(:,:,:,:) :: Pr_Df, T2_Df, Tx_Df, Tm_Df ! (Nlon,Nlat,Nlead,Ntime) !---------------------------------------------------------------------------------------------- ! Verification Scores: Real, Allocatable, dimension(:,:,:,:) :: Pr_ME, Pr_MAE, Pr_RME, Pr_RMAE, Pr_COR, Pr_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: T2_ME, T2_MAE, T2_COR, T2_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tx_ME, Tx_MAE, Tx_COR, Tx_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tm_ME, Tm_MAE, Tm_COR, Tm_pval ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:) ! (Nreg,Mon,LT,Mem) Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:) ! (Nreg,Mon,Nws,Mem) Integer :: Nst, ilon,ilat Character (Len=20),Allocatable :: Stname(:) !---------------------------------------------------------------------------------------------- Real, Allocatable :: beta_store(:,:,:), gamm_store(:,:,:), pzero_store(:,:,:) ! (Nst,Nmon,NLen) - Gamma distribution's parameters !------------------------------------------------ Call Read_OBS ! Read obs data of Pr, T2m, Tx, Tm. Return to Pr_O, T2m_O, Tx_O, Tm_O ! (Nlon,Nlat,Ntime) Allocate (beta_spi_store(Nst,Nmon,Nlen), gamm_spi_store(Nst,Nmon,Nlen), pzero_spi_store(Nst,Nmon,Nlen)) ! (Nst,Nmon,Nlen) Allocate (SPEI_Pars(Nst,0:NMonth,0:2,NLen)) Allocate (DatPr(Ntime,Nst), DatT2m(Ntime,Nst), DatTx(Ntime,Nst), DatTm(Ntime,Nst), Rain(Ntime)) Allocate (SPI(Nst,Ntime,NLen),SPI_O(Nst,Ntime,NLen),SPEI_OT(Nst,Ntime,NLen),SPEI_OH(Nst,Ntime,NLen)) ! T-Thornthwaite; H-Hargreaves Do ist=1,Nst it=0 Do iyr=1,Nyear Do imon=1,Nmon it=it+1 DatPr (it,ist)=Pr_O (ist,iyr,imon) DatT2m(it,ist)=T2m_O(ist,iyr,imon) Enddo Enddo Enddo SPI(:,:,:) = amssng; SPI_O(:,:,:) = Miss; SPEI_OT(:,:,:) = Miss; SPEI_OH(:,:,:) = Miss ! SPI: Call Cal_SPI; SPI_O(:,:,:) = SPI(:,:,:) ! SPEI: SPI(:,:,:) = amssng Call Cal_SPEI_thornthwaite; SPEI_OT(:,:,:) = SPI(:,:,:) ! Replace two years with missing data into origin data: !DatPr(241:288,:) = Miss; DatT2m(241:288,:) = Miss DatPr(Ntime-2*Nmon+1:Ntime,:) = Miss; DatT2m(Ntime-2*Nmon+1:Ntime,:) = Miss Allocate (SPI_O1(Nst,Ntime,NLen),SPEI_OT1(Nst,Ntime,NLen)) ! SPI: SPI(:,:,:) = amssng; SPI_O1(:,:,:) = Miss; SPEI_OT1(:,:,:) = Miss Call Cal_SPI; SPI_O1(:,:,:) = SPI(:,:,:) ! SPEI: SPI(:,:,:) = amssng Call Cal_SPEI_thornthwaite; SPEI_OT1(:,:,:) = SPI(:,:,:) !SPI(:,:,:) = amssng !Call Cal_SPEI__Hargreaves; SPEI_OH(:,:,:) = SPI(:,:,:) !------------------------------------------------ Call Write_to_file ("Test.nc") CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Write_Parameter_of_Distribution (Fname) Character (Len=*) :: Fname Character (Len=100) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm End Subroutine Write_Parameter_of_Distribution !---------------------------------------------------------------------------------------------- Subroutine Write_to_file (Fname) Integer :: imon, id Character (Len=*) :: Fname Character (Len=100) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Real :: StnID(Nst), Xtime(Ntime), Xlen(NLen) lon_nm="stn"; lat_nm="time"; t_nm="scale"; t_unit="months since 1961-01-01 " Do ist=1,Nst; StnID(ist)=ist; Enddo Do it =1,Ntime;Xtime(it)=it; Enddo Xlen(:) = Len_SPI(:) ! Obs cli: var_ln="SPI"; var_unit="-" Call wr_ncdf_3d(Fname,StnID,Xtime,Xlen,Nst,Ntime,NLen,SPI_O,"spi",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="SPEI"; var_unit="-" Call wr_ncdf_3d_apd (Fname,SPEI_OT,"spei",Nst,Ntime,NLen,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="SPI with missing data" Call wr_ncdf_3d_apd (Fname,SPI_O1,"spi1",Nst,Ntime,NLen,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="SPEI with missing data" Call wr_ncdf_3d_apd (Fname,SPEI_OT1,"spei1",Nst,Ntime,NLen,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) End Subroutine Write_to_file !---------------------------------------------------------------------------------------------- Subroutine Read_OBS ! Read observed gridded data of Pr, T2m, Tx, Tm Integer :: j, it, Nyear, ilon, ilat, Nt, Yr0, it1, it2, imon Real, Allocatable :: StnID(:), xyear(:), xmon(:) Character (Len=255) Fname Fname="/home/tanpv/Data/OBS/T_R_Monthly_for_Dr/R_T2m_1961-2018_mon_Tram_KT_61_18_3D.nc" Call rd_ncdf_3d (Fname,StnID,xyear,xmon,Nst,Nyear,Nt,Pr_O,"R") Deallocate (StnID,xyear,xmon) Call rd_ncdf_3d (Fname,StnID,xyear,xmon,Nst,Nyear,Nt,T2m_O,"T2m") Deallocate (StnID,xyear,xmon) Call rd_ncdf_1d (Fname,StnID,Nst,lon,"Lon") Deallocate (StnID) Call rd_ncdf_1d (Fname,StnID,Nst,lat,"Lat") Deallocate (StnID) Call rd_ncdf_1d_Char (Fname,StnID,Nst,Stname,"Stname") End Subroutine Read_OBS !---------------------------------------------------------------------------------------------- Subroutine Correlation(X,Y,N,Corr,prob,Miss) Integer :: N, N1, i, j Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0 Do i=1,N if (X(i)/=Miss.AND.Y(i)/=Miss) then n1=n1+1 Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i) Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i) Rxy=Rxy+X(i)*Y(i) endif Enddo if (N1>2) then Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1) Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb tmp=Rxy/Real(N1) Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy) t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) ) Prob = T_Dist (t, Real(N1-2)) else Rxy=Miss; prob=Miss endif Corr=Rxy End Subroutine Correlation !--------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/DayOfMon.f90" Include "/home/tanpv/Mylibs/Sub_SPI_2022.f90" Include "/home/tanpv/Mylibs/Sub_SPEI_2022.f90" END !--------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Student_Prob.f90"