Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, N365=365 !, Yr1=1961, Yr2=2020 ! 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(:), time(:), xmod(:), Mask(:,:) Integer :: Nlon, Nlat, it, ilon, ilat !------------------------------------------------------------------------------ ! Real, Allocatable, Dimension(:,:,:) :: Pr, Temp, Tmax, Tmin ! Nlon,Nlat,Ntime !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: Tx_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Tx_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Tm_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Tm_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Pr_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Pr_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Qp_tx_obs ! Nlon,Nlat,N365 (Percentile for HW Real, Allocatable, Dimension(:,:,:,:) :: Qp_tx_mod ! Nlon,Nlat,Nmodel,N365 !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: HWevt_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: HWevt_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: HWsev_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: HWsev_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: SPEI_obs ! Nlon,Nlat,Ntime - SPEI with a fixed timescale in days (e.g. 31 days, 91 days) Real, Allocatable, Dimension(:,:,:,:) :: SPEI_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: DrEvt_obs ! Nlon,Nlat,Ntime - Drought events with a fixed threshold Real, Allocatable, Dimension(:,:,:,:) :: DrEvt_mod ! Nlon,Nlat,Nmodel,Ntime !------------------------------------------------------------------------------ ! Real, Allocatable, Dimension(:,:,:) :: Qp_tx ! Nlon,Nlat,N365 (Percentile for HW ! Real, parameter :: p=0.9 ! p=P(X>Qp) ! Integer, Parameter :: Mwd=15 ! data window for calculating percentile Real :: p ! Read from Namelist Integer :: Mwd ! Read from Namelist !------------------------------------------------------------------------------ Character (Len=255) :: Fname,Fname_evt,Fname_sev,Fname_hw,Fname_qp Character (Len=255) :: Obs_Dir, Inp_Dir, Out_Dir, Mask_File Character (Len=255) :: Fevt_obs,Fevt_mod,Fsev_obs,Fsev_mod Integer :: Yr1, Yr2, Nmodel, Timescale, Ntime Character (Len=20), Allocatable :: Modelnm(:) Character (Len=6) :: Tx_nm, Tm_nm, Pr_nm ! Varname in files Logical :: status, Empty, calc_qp ! Integer, Allocatable :: TimeScale(:) Integer :: Nlo,Nla,Nt,Nz !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, SPEI_out, Rain, T2m, Tx, Tm !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Mask_File,Obs_Dir,Inp_Dir,Out_Dir Namelist /Data_Info/ Yr1, Yr2, Timescale, calc_qp, p, Mwd Namelist /Modelname/ Tx_nm, Tm_nm, Pr_nm !------------------------------------------------------------------------------ Open (1, file="Namelist_CDH_Vn_GRD.txt",Status="old") Read (1, File_Dir_Names) Read (1, Data_Info) Read (1, Modelname) Ntime=(Yr2-Yr1+1)*N365 ! Read landmask data for SEA region with resolution of CORDEX SEA 2024 (~0.25 deg) Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Mask,"reg") ! Land points !------------------------------------------------------------------------ ! Read Tx, Tm and Pr obs Call Read_Tx_Tm_Pr ! Read obs data !------------------------------------------------------------------------ ! Since the calculation of Qp is time-consuming, it is necessary to store the results in files for future use: Fname = trim(Out_Dir)//"Qp_obs.nc" if (calc_qp) then ! Need to calculate Qp again Call Calc_percentile_Tx_obs Call Write_Qp (Fname) else ! Read from file Call rd_ncdf_3d (Fname,xlon,xlat,time,nlo,nla,nt,Qp_tx_obs,"qp") Deallocate (xlon,xlat,time) endif !------------------------------------------------------------------------ Call Calc_HW_events Fsev_obs=trim(Out_Dir)//"HW_sev_obs.nc" Call Write_HWevt_and_HWsev (Fsev_obs) Call Calc_spei_Ha (TimeScale) Fname = trim(Out_Dir)//"SPEI_obs.nc" Call Write_Drought_IDX (Fname) !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_HW_events Integer :: imod,ilon,ilat,it,iday !-------- Obs: ----------- !Allocate (HWevt_obs(Nlon,Nlat,Ntime)); HWevt_obs=Miss Allocate (HWsev_obs(Nlon,Nlat,Ntime)); HWsev_obs=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do it=1,Ntime iday=mod(it,N365); if (iday==0) iday=N365 if (Tx_obs(ilon,ilat,it)/=Miss.AND.Qp_tx_obs(ilon,ilat,iday)/=Miss) then if (Tx_obs(ilon,ilat,it)>Qp_tx_obs(ilon,ilat,iday)) then !HWevt_obs(ilon,ilat,it)=Tx_obs(ilon,ilat,it)-Qp_tx_obs(ilon,ilat,iday) ! Severity HWsev_obs(ilon,ilat,it)=Tx_obs(ilon,ilat,it)-Qp_tx_obs(ilon,ilat,iday) ! Severity else !HWevt_obs(ilon,ilat,it)=0.0 HWsev_obs(ilon,ilat,it)=0.0 endif endif Enddo endif Enddo Enddo End Subroutine Calc_HW_events !------------------------------------------------------------------------------ Subroutine Calc_spei_Ha (Timescale) ! Using Hargreaves method Integer :: TimeScale, ilon,ilat,it,imod, ilen Real :: xlat0 print*,"In Calc_spei_all ..." Allocate (Tx(Ntime), Tm(Ntime), Rain(Ntime), SPEI_out(Ntime)) Allocate (SPEI_obs(Nlon,Nlat,Ntime)) !, SPEI_mod(Nlon,Nlat,Nmodel,Ntime)) SPEI_obs=Miss; !SPEI_mod=Miss ! OBS: Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then xlat0 = Lat(ilat) Tx (1:Ntime) = Tx_obs(ilon,ilat,1:Ntime) Tm (1:Ntime) = Tm_obs(ilon,ilat,1:Ntime) Rain(1:Ntime) = Pr_obs(ilon,ilat,1:Ntime) Call Calc_SPEI_Harg_Daily_One_Stn (Tx, Tm, Rain, xlat0, Ntime, Timescale, SPEI_out, Miss) !print"(20F7.2)",Tx (1:Ntime)-Tm (1:Ntime) ; pause !print"(20F7.2)",Tm (1:Ntime); pause !print"(20F7.2)",Rain(1:Ntime); pause !print"(20F7.2)",SPEI_out(1:Ntime); pause SPEI_obs(ilon,ilat,1:Ntime) = SPEI_out(1:Ntime) endif Enddo Enddo End Subroutine Calc_spei_Ha !------------------------------------------------------------------------------ Subroutine Write_Drought_IDX (Fobs) Character (Len=*) :: Fobs 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,Nz,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg if (Allocated(time)) Deallocate (time) Allocate (time(Ntime)) t_unit="days since 1980-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,Ntime; time(i)=i-1; Enddo Write(Var_ln,"(A,I2.2,A)") "The obs drought index (SPEI) with timescale of ",Timescale, " days calculated from VN OBS grid data" Call wr_ncdf_3d (Fobs,Lon,Lat,time,Nlon,Nlat,Ntime,SPEI_obs,"spei",t_unit, var_ln,"-", lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_3d_apd (Fevt,SPEI_obs,"obs",Nlon,Nlat,Ntime,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_4d_apd (Fevt,X4d,"obs", Nlon,Nlat,Nz,Ntime, & ! Var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) End Subroutine Write_Drought_IDX !------------------------------------------------------------------------------ Subroutine Write_HWevt_and_HWsev (Fsev_obs) Character (Len=*) :: Fsev_obs 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,Nz,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg if (Allocated(time)) Deallocate (time) Allocate (time(Ntime)) t_unit="days since 1980-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,Ntime; time(i)=i-1; Enddo Var_ln="The obs HW severity calculated from Tx" Call wr_ncdf_3d (Fsev_obs,Lon,Lat,time,Nlon,Nlat,Ntime,HWsev_obs,"sev",t_unit, var_ln,"C deg", lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_HWevt_and_HWsev !------------------------------------------------------------------------------ Subroutine Write_Qp (Fobs) !,Fmod) ! write percentiles of obs and models: Character (Len=*) :: Fobs !,Fmod 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 :: xmod(Nmodel) if (Allocated(time)) Deallocate (time) Allocate (time(N365)) t_unit="days since 1980-01-01 "; var_unit="C deg"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,N365; time(i)=i-1; Enddo !Do i=1,Nmodel; xmod(i)=i; Enddo !Var_ln="The model percentiles. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" !Call wr_ncdf_4d (Fmod,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,N365,Qp_tx_mod,"qp", & ! t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The obs Tx percentiles." Call wr_ncdf_3d (Fobs,Lon,Lat,time,Nlon,Nlat,N365,Qp_tx_obs,"qp",t_unit, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_Qp !------------------------------------------------------------------------------ Subroutine Calc_percentile_Tx_obs Integer :: ilon,ilat,it Real :: Qp(N365) Allocate(Tx(Ntime),Qp_tx_obs(Nlon,Nlat,N365)); Qp_tx_obs=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Tx (1:Ntime) = Tx_obs(ilon,ilat,1:Ntime) Call Daily_per_HW (Tx,Ntime,p,Mwd,Qp,Miss) ! p and Mwd are needed to be defined Qp_tx_obs(ilon,ilat,1:N365)=Qp(1:N365) endif Enddo Enddo Deallocate (Tx) End Subroutine Calc_percentile_Tx_obs !------------------------------------------------------------------------------ Subroutine Read_Tx_Tm_Pr !Real, parameter :: FillValue=-9999.0 Integer :: iyr,it,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X3d(:,:,:), xtime(:), xlon(:),xlat(:) !,xtime(:) Fname=trim(Obs_Dir)//"tx_daily_1980-2021_0.25.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,Tx_obs,"tx") Deallocate (xlon,xlat,xtime) Fname=trim(Obs_Dir)//"tm_daily_1980-2021_0.25.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,Tm_obs,"tm") Deallocate (xlon,xlat,xtime) Fname=trim(Obs_Dir)//"pr_daily_1980-2021_0.25.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,Pr_obs,"pr") Deallocate (xlon,xlat,xtime) if (Nt/=Ntime) then print*,"Something wrong during reading obs data." stop endif End Subroutine Read_Tx_Tm_Pr !------------------------------------------------------------------------------ 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_SPEI_Daily_One_Stn_2024.f90" Include "/home/tanpv/Mylibs/Sub_Daily_Percentile.f90" !-------------------------------------------------------------------- END ! Program