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 !Integer, Parameter :: Nyear=Yr2-Yr1+1 !Integer, Parameter :: Nlen=1 ! 3 ! Number of timascales !Integer, Parameter :: TimeScale(Nlen)=(/31/) !7,15,31/) ! in day Real, parameter :: Dr_Thre = -1.0 ! Moderate Drought ! 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(:), Mask(:,:) Integer :: Nlon, Nlat, it, ilon, ilat !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:,:) :: SPEI ! Nlon,Nlat,Nlen,Ntime Real, Allocatable, Dimension(:,:,:) :: HWevt ! Nlon,Nlat,Ntime HW events Real, Allocatable, Dimension(:,:,:,:) :: CDH ! Nlon,Nlat,Nlen,Ntime (Evts) Real, Allocatable, Dimension(:,:,:,:) :: CDH_Frq ! Nlon,Nlat,Nlen,Nyear (Frequency - Num of evts in each year) Real, Allocatable, Dimension(:,:,:,:) :: CDH_Dur ! Nlon,Nlat,Nlen,Nyear (Mean duration - Average the num of days) Real, Allocatable, Dimension(:,:,:,:) :: CDH_Dmx ! Nlon,Nlat,Nlen,Nyear (Max duration - Number of days of longest duration)) Real, Allocatable, Dimension(:,:,:,:) :: CDH_Ndy ! Nlon,Nlat,Nlen,Nyear (Number of CDH days in the year) !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: Sen_Frq, Sen_Dur, Sen_Dmx, Sen_Ndy ! Nlon,Nlat,Nlen - Sen's Trends Real, Allocatable, Dimension(:,:,:) :: Prb_Frq, Prb_Dur, Prb_Dmx, Prb_Ndy ! Nlon,Nlat,Nlen - MK !------------------------------------------------------------------------------ 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 !------------------------------------------------------------------------------ Character (Len=255) :: Fname,Fname_th,Fname_ha,Fname_hw,Fname_qp 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, SPEI_out, Rain, T2m, Tx, Tm !------------------------------------------------------------------------------ Integer :: Yr1, Yr2, Ntime, Nyear, Nlen Integer, Allocatable :: TimeScale(:) !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Mask_File,Inp_Dir,Out_Dir Namelist /Data_Info/ Yr1, Yr2, Nlen Namelist /Dr_Timescale/ TimeScale !------------------------------------------------------------------------------ Open (1, file="Namelist_Dr_HW_SEA.txt",Status="old") Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (TimeScale(Nlen)) Read (1, Dr_Timescale) Nyear=Yr2-Yr1+1 Inp_Dir = trim(Out_Dir) !"/home/tanpv/Progs/Dr_Daily/Res/" !Out_Dir = trim(Inp_Dir) !"/home/tanpv/Progs/Dr_Daily/Res/" !Fname = "/home/tanpv/Data/AWC/SEA_landmask_ERA5.nc" ! Land mask file created from APH with preci>=0 Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Mask,"land") ! Land points !------------------------------------------------------------------------ Fname_th=trim(Inp_Dir)//"Dr_Th_SEA.nc"; Fname_hw=trim(Inp_Dir)//"HW_Evt_SEA.nc"; Call Read_Dr_HW (Fname_th,Fname_hw,"spei_th") Allocate (CDH(Nlon,Nlat,Nlen,Ntime),CDH_Frq(Nlon,Nlat,Nlen,Nyear)) Allocate (CDH_Dur(Nlon,Nlat,Nlen,Nyear),CDH_Dmx(Nlon,Nlat,Nlen,Nyear)) Allocate (CDH_Ndy(Nlon,Nlat,Nlen,Nyear)) Allocate (Sen_Frq(Nlon,Nlat,Nlen), Sen_Dur(Nlon,Nlat,Nlen), Sen_Dmx(Nlon,Nlat,Nlen)) Allocate (Prb_Frq(Nlon,Nlat,Nlen), Prb_Dur(Nlon,Nlat,Nlen), Prb_Dmx(Nlon,Nlat,Nlen)) Allocate (Sen_Ndy(Nlon,Nlat,Nlen), Prb_Ndy(Nlon,Nlat,Nlen)) CDH = Miss; CDH_Frq = Miss; CDH_Dur = Miss; CDH_Dmx = Miss; CDH_Ndy = Miss Sen_Frq = Miss; Sen_Dur = Miss; Sen_Dmx = Miss; Sen_Ndy = Miss; Prb_Frq = Miss; Prb_Dur = Miss; Prb_Dmx = Miss; Prb_Ndy = Miss; Call Calc_CDH_evts ! Determine CDH events Call Calc_CDH_Chars ! Calculate CDH Characteristics Call Calc_Trends ! Calculate CDH Characteristics' trends Fname = trim(Out_Dir)//"CDH_Chars_1.nc" Fname_hw = trim(Out_Dir)//"CDH_Trend_Chars_1.nc" Call Write_CDH_Chars (Fname,Fname_hw) !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_Trends Integer :: it,it1,it2,i,ilon,ilat,ilen Real :: Sen, Prb Real :: X1d(Nyear) Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then X1d(1:Nyear) = CDH_Frq(ilon,ilat,ilen,1:Nyear) Call MK_Sensl(X1d,Nyear,Sen,prb,Miss) Sen_Frq(ilon,ilat,ilen) = Sen; Prb_Frq(ilon,ilat,ilen) = prb X1d(1:Nyear) = CDH_Dur(ilon,ilat,ilen,1:Nyear) Call MK_Sensl(X1d,Nyear,Sen,prb,Miss) Sen_Dur(ilon,ilat,ilen) = Sen; Prb_Dur(ilon,ilat,ilen) = prb X1d(1:Nyear) = CDH_Dmx(ilon,ilat,ilen,1:Nyear) Call MK_Sensl(X1d,Nyear,Sen,prb,Miss) Sen_Dmx(ilon,ilat,ilen) = Sen; Prb_Dmx(ilon,ilat,ilen) = prb X1d(1:Nyear) = CDH_Ndy(ilon,ilat,ilen,1:Nyear) Sen_Ndy(ilon,ilat,ilen) = Sen; Prb_Ndy(ilon,ilat,ilen) = prb endif Enddo Enddo Enddo End Subroutine Calc_Trends !------------------------------------------------------------------------------ Subroutine Calc_CDH_Chars Integer :: it,it1,it2,i,ilon,ilat,ilen Real :: X1d(N365) Real :: M1(N365), IMAX1, Spells(N365), N_Spell_1, I_Start(N365), I_End(N365) Real :: tmp,tmpx,duration,numday Integer :: Nspls,Nspl3 Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do iyr=1,Nyear it1=(iyr-1)*N365+1; it2=iyr*365 X1d(1:N365) = CDH(ilon,ilat,ilen,it1:it2) Call NSpell_and_Leng (X1d, N365, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss) Nspls = N_Spell_1 ! So dot tmp = 0 ! Do dai trung binh tmpx = 0 ! Do dai Max Nspl3 = 0 ! So dot co duration >=3 days numday = 0 ! Number of CDH days Do i=1,Nspls duration = Spells(i) ! I_End(i) - I_Start(i) + 1 if (duration >= 3.) then Nspl3 = Nspl3 +1 ! Dem so dot co duration >= 3 days tmp = tmp + duration ! Tong duration if (tmpx < duration) tmpx = duration ! Max duration endif Enddo numday = tmp if (Nspl3>0) then ! CDH_Frq(ilon,ilat,ilen,iyr) = Nspl3 CDH_Dur(ilon,ilat,ilen,iyr) = tmp/Nspl3 CDH_Dmx(ilon,ilat,ilen,iyr) = tmpx else CDH_Frq(ilon,ilat,ilen,iyr) = 0 CDH_Dur(ilon,ilat,ilen,iyr) = 0 CDH_Dmx(ilon,ilat,ilen,iyr) = 0 endif CDH_Ndy(ilon,ilat,ilen,iyr) = numday Enddo endif Enddo Enddo Enddo End Subroutine Calc_CDH_Chars !------------------------------------------------------------------------------ Subroutine Calc_CDH_evts Integer :: it,i,ilon,ilat,ilen Logical :: Drou, HW, Undef, NoDr, NoHW Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do it=1,Ntime Drou = SPEI(ilon,ilat,ilen,it)/=Miss.AND.SPEI(ilon,ilat,ilen,it)<=Dr_Thre HW = HWevt(ilon,ilat,it)==1.0 Undef = SPEI(ilon,ilat,ilen,it)==Miss.OR.HWevt(ilon,ilat,it)==Miss !NoDr = SPEI(ilon,ilat,ilen,it)>Dr_Thre !NoHW = HWevt(ilon,ilat,it)==0.0 ! if (Drou.AND.HW) then ! CDH(ilon,ilat,ilen,it) = 1 ! else if (Undef) then ! CDH(ilon,ilat,ilen,it) = Miss ! else ! CDH(ilon,ilat,ilen,it) = 0.0 ! endif if (SPEI(ilon,ilat,ilen,it)/=Miss.AND.SPEI(ilon,ilat,ilen,it)<=Dr_Thre.AND.HWevt(ilon,ilat,it)==1.0) then CDH(ilon,ilat,ilen,it) = 1 else if (SPEI(ilon,ilat,ilen,it)>Dr_Thre.OR.HWevt(ilon,ilat,it)==0) then CDH(ilon,ilat,ilen,it) = 0 endif !print*,Int(CDH(ilon,ilat,ilen,it)) Enddo endif Enddo Enddo Enddo End Subroutine Calc_CDH_evts !------------------------------------------------------------------------------ Subroutine Write_CDH_Chars (Fname,Ftrend) Character (Len=*) :: Fname,Ftrend 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 :: xscale(Nlen) !=(/31/) ! !7,15,31/) xscale(:) = TimeScale(:) if (Allocated(time)) Deallocate (time) Allocate (time(Nyear)) t_unit="years since 1961-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="year" Do i=1,Nyear; time(i)=i; Enddo Var_ln="Coumpound Drought-Heatwave frequency (number of events/year. Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d (Fname,Lon,Lat,xscale,time, Nlon,Nlat,Nlen,Nyear,CDH_Frq,"frq", & t_unit,Var_ln,"events/year",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Coumpound Drought-Heatwave mean duration (number of days). Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d_apd (Fname,CDH_Dur,"dur",Nlon,Nlat,Nlen,Nyear,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Coumpound Drought-Heatwave Maximum duration (number of days). Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d_apd (Fname,CDH_Dmx,"dmx",Nlon,Nlat,Nlen,Nyear,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Number of days of the Coumpound Drought-Heatwave. Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d_apd (Fname,CDH_Ndy,"ndy",Nlon,Nlat,Nlen,Nyear,Var_ln,"-", & lon_nm, lat_nm, z_nm, t_nm,Miss) !------------------------------------------------------------------------------ t_nm="timescale" Var_ln="Trend of changes in CDH_Frq" ! Sen_Frq, Sen_Dur, Sen_Dmx Call wr_ncdf_3d (Ftrend,Lon,Lat,xscale, Nlon,Nlat,Nlen,Sen_Frq,"sen_frq",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Trend of changes in CDH_Dur" Call wr_ncdf_3d_apd (Ftrend,Sen_Dur,"sen_dur",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Trend of changes in CDH_Dmx" Call wr_ncdf_3d_apd (Ftrend,Sen_Dmx,"sen_dmx",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Trend of changes in CDH_Ndy" Call wr_ncdf_3d_apd (Ftrend,Sen_Ndy,"sen_ndy",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Statistical significance of trend of CDH_Frq" Call wr_ncdf_3d_apd (Ftrend,Prb_Frq,"prb_frq",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Statistical significance of trend of CDH_Dur" Call wr_ncdf_3d_apd (Ftrend,Prb_Dur,"prb_dur",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Statistical significance of trend of CDH_Dmx" Call wr_ncdf_3d_apd (Ftrend,Prb_Dmx,"prb_dmx",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Statistical significance of trend of CDH_Ndy" Call wr_ncdf_3d_apd (Ftrend,Prb_Ndy,"prb_ndy",Nlon,Nlat,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Deallocate (time) End Subroutine Write_CDH_Chars !------------------------------------------------------------------------------ Subroutine Read_Dr_HW (Fname_dr,Fname_hw,dr_idx) Character (Len=*) :: Fname_dr,Fname_hw,dr_idx Integer :: it,ilon,ilat,ilead,Nlo,Nla,Nlev,Nt,icas, it1, it2 Real, Allocatable :: X4d(:,:,:,:), xtime(:), xlon(:), xlat(:), xlev(:) !,xtime(:) Call rd_ncdf_4d (Fname_dr,xlon,xlat,xlev,xtime,nlo,nla,nlev,nt,SPEI,trim(dr_idx)) Deallocate (xlon,xlat,xlev,xtime) Call rd_ncdf_3d (Fname_hw,xlon,xlat,xtime,Nlo,Nla,Nt,HWevt,"hw") Deallocate (xlon,xlat,xtime) Ntime = Nt End Subroutine Read_Dr_HW !------------------------------------------------------------------------------ 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_NSpell_and_Leng.f90" Include "/home/tanpv/Mylibs/Sub_Trend_MK.f90" !Include "/home/tanpv/Mylibs/Sub_SPEI_Daily_One_Stn_2024.f90" !Include "/home/tanpv/Mylibs/Sub_Daily_Percentile.f90" !-------------------------------------------------------------------- END ! Program