Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI ! use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99 Integer, Parameter :: Nmon=12, N365=365, Yr1=1961, Yr2=2011, Ntime=(Yr2-Yr1+1)*N365 Integer, Parameter :: Nday=31, Nyear=(Yr2-Yr1+1) Integer, Parameter :: TimeScale=7 Integer, parameter :: DayOfMn(Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! Distribution function's parameters: Real, Allocatable :: beta_spi(:), gamm(:), pzero(:) Real :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2) Real, Allocatable :: Lon(:), Lat(:), xlon(:), xlat(:) Integer :: Nlon, Nlat, it, ilon, ilat Real :: xlat0 !------------------------------------------------------------------------------ ! Real, Allocatable, Dimension(:,:,:) :: Pr, Temp, Tmax, Tmin ! Nlon,Nlat,Ntime ! Real, Allocatable, Dimension(:,:,:) :: SPEI ! Nlon,Nlat,Ntime !------------------------------------------------------------------------------ Real, Dimension (Nyear,Nmon,Nday) :: EV, Pr, T2m, Tx, Tm, X3d Real, Dimension (Ntime) :: pet0, pr0, t2m0, tx0,tm0, SPEI_out, spei_th, spei_ha Real, Dimension (Ntime) :: Y ! for testing Integer :: Ny,Nyr ! for testing Real, Allocatable :: Mean(:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File Logical :: status, Empty !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, xspi Real :: evth, evha, HI !------------------------------------------------------------------------------ Inp_Dir = "/home/tanpv/OBS_Stn/Daily/T0_Huong_2011/" Out_Dir = "/home/tanpv/Progs/Dr_Daily/Res/" !------------------------------------------------------------------------ xlat0 = 21.028511 Fname=trim(Inp_Dir)//"R_HANOI.txt" Call Read_One_Sta (Fname, Pr) Fname=trim(Inp_Dir)//"T2m_HANOI.txt" Call Read_One_Sta (Fname, T2m) Fname=trim(Inp_Dir)//"Tx_HANOI.txt" Call Read_One_Sta (Fname, Tx) Fname=trim(Inp_Dir)//"Tm_HANOI.txt" Call Read_One_Sta (Fname, Tm) Fname=trim(Inp_Dir)//"BH_HANOI.txt" Call Read_One_Sta (Fname, EV) ! Call Reshape_3d_to_1d (EV,pet0) Call Reshape_3d_to_1d (Pr,pr0) Call Reshape_3d_to_1d (T2m,t2m0) Call Reshape_3d_to_1d (Tx,tx0) Call Reshape_3d_to_1d (Tm,tm0) HI = Heat_Index (T2m0, Ntime, Miss) !print*,"HI=", HI; pause !Do it=1,Ntime ! iday = mod(it,N365); if (iday==0) iday=N365 ! evth = PET_daily_Thorn (t2m0(it), HI, xlat0, iday) ! evha = PET_daily_Harg (tx0(it), tm0(it), xlat0, iday) ! print*,it,iday,evth,evha,pr0(it) !Enddo !pause !-------------------------------------------- ! Call Calc_SPEI_Thorn_Daily_One_Stn (T2m0, pr0, xlat0, Ntime, TimeScale, SPEI_out, Miss) ! spei_th(:)=SPEI_out(:) ! Call Calc_SPEI_Harg_Daily_One_Stn (Tx0, Tm0, pr0, xlat0, Ntime, TimeScale, SPEI_out, Miss) ! spei_ha(:)=SPEI_out(:) ! Open (3,file="HN_SPEI.txt",status="unknown") ! write (3,"(A8,3A10)") "Date","pet","spei_th","spei_ha" ! it=0 ! Do iyr=1,Nyear ! Do imon=1,Nmon ! Do iday=1,DayOfMn(imon) ! it=it+1 ! write(3,"(I4.4,2I2.2,10F8.3)") 1960+iyr,imon,iday,pet0(it),spei_th(it),spei_ha(it) ! Enddo ! Enddo ! Enddo ! close(3) !-------------------------------------------- ! Test for extract daily to monthly data: Do imon=1,Nmon Call Monthly_Data_Filter (pr0,Ntime,Nyr,Y,Ny,imon,Mean,Miss) xtmp=0 Do j=1,Ny xtmp=xtmp+Y(j) print"(2I6,2F10.1)", imon,j,Y(j),xtmp Enddo print*,"===" print"(I6,F10.1)",Ny,Sumx(Y,Ny,Miss) pause Do j=1,Nyr print*,j,Mean(j) Enddo pause Enddo !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Reshape_3d_to_1d (X3d,X1d) Integer :: iyr,imon,iday,it Real, Dimension (Nyear,Nmon,Nday) :: X3d Real :: X1d(Ntime) it=0 Do iyr=1,Nyear Do imon=1,Nmon Do iday=1,DayOfMn(imon) it=it+1 X1d(it)=X3d(iyr,imon,iday) Enddo Enddo Enddo End Subroutine Reshape_3d_to_1d !------------------------------------------------------------------------------ Subroutine Read_One_Sta (Fname, X3d) Character (Len=*) :: Fname Integer :: imon, iday, day, IO, iyr, yy Real, Dimension (Nyear,Nmon,Nday) :: X3d Real :: w1(Nmon) Open (1,file=Fname,Status="Old") IO = 0; Read (1,*); Read (1,*) ! Bo 2 dong dau Do iyr=1,Nyear !if (IO /= 0) exit read (1,*, IOStat=IO) yy !if (iyr > Yr2) Exit Do iday=1,Nday read (1,*,IOStat=IO) day, (w1(imon), imon=1,Nmon) Do imon=1,Nmon X3d(iyr,imon,iday) = w1(imon) Enddo Enddo Enddo Close (1) End Subroutine Read_One_Sta !------------------------------------------------------------------------------ Subroutine Write_GeoExt (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str Integer :: it,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),xreg(Nreg),xlev(Dr_Lev),xtime(Ntime) Logical :: append write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Geo_Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="reg"; lat_nm="dr_lev"; z_nm="lead"; t_nm="time" ! Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo ! Do i=1,Dr_Lev; xlev(i)=i; Enddo ! Do i=1,Nreg; xreg(i)=i; Enddo ! Do i=1,Ntime; xtime(i)=i; Enddo Var_ln="Model drought geographic extent. reg=1..7 subregs; dr_lev=1..3: mild,moder,extreme; lead=1..6 months; time=1/1983-12/2020" Var_nm="geof" ! Call wr_ncdf_4d (Fname,xreg,xlev,xlead,xtime, Nreg,Dr_Lev,Nlead,Ntime,Geo_F,trim(Var_nm), & ! t_unit,Var_ln,"%",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Observed drought geographic extent. reg=1..7 subregs; dr_lev=1..3: mild,moder,extreme; lead=1..6 months; time=1/1983-12/2020" Var_nm="geob" ! Call wr_ncdf_3d_apd (Fname,Geo_O,trim(Var_nm),Nreg,Dr_Lev,Ntime,Var_ln,"%",lon_nm,lat_nm,t_nm,Miss) End Subroutine Write_GeoExt !------------------------------------------------------------------------------ ! 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/Sum.f90" ! Include "/home/tanpv/Mylibs/Sub_Monthly_Data_Filter.f90" ! Include "/home/tanpv/Mylibs/Sub_SPEI_Daily_One_Stn_2024.f90" Include "Sub_SPEI_Daily_One_Stn_test_2025.f90" !-------------------------------------------------------------------- END ! Program