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 :: Ntime=(Yr2-Yr1+1)*N365 ! Integer, Parameter :: Nlen=2 ! 3 ! Number of timascales ! Integer, Parameter :: TimeScale(Nlen)=(/31,91/) !7,15,31/) ! in day ! 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(:,:,:) :: Pr, Temp, Tmax, Tmin ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: SPEI_th, SPEI_ha ! Nlon,Nlat,Nlen,Ntime Real, Allocatable, Dimension(:,:,:) :: HWevt ! Nlon,Nlat,Ntime - HW events Real, Allocatable, Dimension(:,:,:) :: HWsev ! Nlon,Nlat,Ntime - HW severity = Tx-Qp 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 Integer :: Yr1, Yr2, Ntime, Nlen Integer, Allocatable :: TimeScale(:) !------------------------------------------------------------------------------ ! 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,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) Ntime=(Yr2-Yr1+1)*N365 print*, Yr1, Yr2, Ntime, Nlen print*, TimeScale !pause ! Inp_Dir = "/home/tanpv/Data/ERA5/SEA/" ! Out_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 !------------------------------------------------------------------------ Call Read_all_data ! Read all data of "mod","cli","ann","mlr" Call Invert_Lat (Pr) ! Invert Latitude to Mask file's regulation Call Invert_Lat (Temp) Call Invert_Lat (Tmax) Call Invert_Lat (Tmin) Allocate (SPEI_out(Ntime), Rain(Ntime), T2m(Ntime), Tx(Ntime), Tm(Ntime)) Allocate (SPEI_th(Nlon,Nlat,Nlen,Ntime), SPEI_ha(Nlon,Nlat,Nlen,Ntime)) Allocate (Qp_tx(Nlon,Nlat,N365),HWevt(Nlon,Nlat,Ntime)) Allocate (HWsev(Nlon,Nlat,Ntime)) SPEI_out=Miss; SPEI_th=Miss; SPEI_ha=Miss; Rain=Miss; T2m=Miss; Tx=Miss; Tm=Miss; Qp_tx=Miss; HWevt=Miss; HWsev=Miss Call Calc_spei_all Call Calc_percentile_Tx Call Calc_HW_events !Fname=trim(Out_Dir)//"Dr_HW_SEA.nc" Fname_th=trim(Out_Dir)//"Dr_Th_SEA.nc"; Fname_ha=trim(Out_Dir)//"Dr_Ha_SEA.nc"; Fname_hw=trim(Out_Dir)//"HW_Evt_SEA.nc"; Fname_qp=trim(Out_Dir)//"Q90_Tx_SEA.nc" Fname =trim(Out_Dir)//"HW_Sev_SEA.nc" Call Write_Daily_Dr_IDX_and_HWevt (Fname_th,Fname_ha,Fname_hw,Fname,Fname_qp) !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_HW_events Integer :: ilon,ilat,it,iday HWevt = 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 (Tmax(ilon,ilat,it)/=Miss.AND.Qp_tx(ilon,ilat,iday)/=Miss) then if (Tmax(ilon,ilat,it)>=Qp_tx(ilon,ilat,iday)) then HWevt(ilon,ilat,it)=1 HWsev(ilon,ilat,it)=Tmax(ilon,ilat,it)-Qp_tx(ilon,ilat,iday) else HWevt(ilon,ilat,it)=0.0 HWsev(ilon,ilat,it)=0.0 endif endif Enddo !print*,"HW evts: ",HWevt(ilon,ilat,100:200); pause endif Enddo Enddo End Subroutine Calc_HW_events !------------------------------------------------------------------------------ Subroutine Calc_percentile_Tx Integer :: ilon,ilat,it Real :: Qp(N365) Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Tx (1:Ntime) = Tmax(ilon,ilat,1:Ntime) Call Daily_per_HW (Tx,Ntime,p,Mwd,Qp,Miss) Qp_tx(ilon,ilat,1:N365)=Qp(1:N365) !print*,"Qp: ",Qp_tx(ilon,ilat,1:N365); pause endif Enddo Enddo End Subroutine Calc_percentile_Tx !------------------------------------------------------------------------------ Subroutine Calc_spei_all Integer :: ilon,ilat,it, ilen Real :: xlat0 ! Pr, Temp, Tmax, Tmin ! Nlon,Nlat,Ntime print*,"In Calc_spei_all ..." Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then xlat0 = Lat(ilat) Rain(1:Ntime) = Pr (ilon,ilat,1:Ntime) T2m (1:Ntime) = Temp(ilon,ilat,1:Ntime) Call Calc_SPEI_Thorn_Daily_One_Stn (T2m, Rain, xlat0, Ntime, TimeScale(ilen), SPEI_out, Miss) SPEI_th(ilon,ilat,ilen,1:Ntime) = SPEI_out(1:Ntime) Tx (1:Ntime) = Tmax(ilon,ilat,1:Ntime) Tm (1:Ntime) = Tmin(ilon,ilat,1:Ntime) Call Calc_SPEI_Harg_Daily_One_Stn (Tx, Tm, Rain, xlat0, Ntime, TimeScale(ilen), SPEI_out, Miss) SPEI_ha(ilon,ilat,ilen,1:Ntime) = SPEI_out(1:Ntime) ! print to test: !Do it=100,200 ! print"(4F6.2)",SPEI_th(ilon,ilat,ilen,it),SPEI_ha(ilon,ilat,ilen,it) !Enddo endif Enddo Enddo Enddo End Subroutine Calc_spei_all !------------------------------------------------------------------------------ Subroutine Write_Daily_Dr_IDX_and_HWevt (Fname_th,Fname_ha,Fname_hw,Fname_sev,Fname_qp) Character (Len=*) :: Fname_th,Fname_ha,Fname_hw,Fname_qp,Fname_sev 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(Ntime)) t_unit="days since 1961-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="day" Do i=1,Ntime; time(i)=i; Enddo Var_ln="Daily drought SPEI index using Thornthwaite PET. Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d (Fname_th,Lon,Lat,xscale,time, Nlon,Nlat,Nlen,Ntime,SPEI_th,"spei_th", & t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Daily drought SPEI index using Hargreaves PET. Z=Timescales of 7,15 and 31 days" Call wr_ncdf_4d (Fname_ha,Lon,Lat,xscale,time, Nlon,Nlat,Nlen,Ntime,SPEI_ha,"spei_ha", & t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) !Call wr_ncdf_4d_apd (Fname,SPEI_ha,"spei_ha",Nlon,Nlat,Nlen,Ntime,Var_ln,"-", & ! lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Heatwave events (daily maximum temperature >= q90)" !Call wr_ncdf_3d_apd (Fname,HWevt,"hw",Nlon,Nlat,Ntime, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Call wr_ncdf_3d (Fname_hw, Lon,Lat,time,Nlon,Nlat,Ntime,HWevt,"hw", t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="Heatwave severity = Tx - Qp" Call wr_ncdf_3d (Fname_sev,Lon,Lat,time,Nlon,Nlat,Ntime,HWsev,"sev",t_unit, var_ln, "C deg", lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_3d_apd (Fname_sev,HWsev,"sev",Nlon,Nlat,Ntime, var_ln, "C deg", lon_nm, lat_nm, t_nm, Miss) Var_ln="90%-percentile of daily Tmax (Q90) - time = 1..365" Deallocate (time) Allocate (time(N365)) Do i=1,N365; time(i)=i; Enddo !Call wr_ncdf_3d_apd (Fname,Qp_tx,"q90",Nlon,Nlat,Ntime, var_ln,"C", lon_nm, lat_nm, t_nm, Miss) Call wr_ncdf_3d (Fname_qp,Lon,Lat,time,Nlon,Nlat,N365,Qp_tx,"q90",t_unit, var_ln,"C", lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_Daily_Dr_IDX_and_HWevt !------------------------------------------------------------------------------ Subroutine Invert_Lat (X) Integer :: ilat, j Real :: X(Nlon,Nlat,Ntime), Y(Nlon,Nlat,Ntime) print*,Ntime,Nlon,Nlat Do ilat=1,Nlat j = Nlat-ilat+1 Y(:,ilat,:) = X(:,j,:) Enddo X = Y End Subroutine Invert_Lat !------------------------------------------------------------------------------ Subroutine Read_all_data ! (Model_nm) !Character (Len=*) :: Model_nm Integer :: it,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X3d(:,:,:), xtime(:), xlon(:),xlat(:),xmem(:) !,xtime(:) Real, parameter :: Mis_val = -32767.0 ! Pr, Temp, Tmax, Tmin if (Allocated(X3d)) Deallocate (X4d); if (Allocated(xlon)) Deallocate (xlon); if (Allocated(xlat)) Deallocate (xlat) if (Allocated(xtime)) Deallocate (xtime) it1 = 0 Do iyr=Yr1,Yr2 write (Fname,"(A,I4.4,A)") trim(Inp_Dir)//"pr_",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tp") if (iyr==Yr1) then !Nlon=Nlo; Nlat=Nla !Allocate (Lon(Nlon),Lat(Nlat)) Allocate (Pr (Nlon,Nlat,Ntime),Temp(Nlon,Nlat,Ntime)) Allocate (Tmax(Nlon,Nlat,Ntime),Tmin(Nlon,Nlat,Ntime)) Pr = Miss; Temp = Miss; Tmax = Miss; Tmin = Miss !Lon = xlon; Lat = xlat endif it1 = (iyr-Yr1)*N365+1 it2 = it1 - 1 + N365 !print*,it1,it2; pause Where (X3d==Mis_val) X3d = Miss Where (X3d/=Miss) X3d = X3d*1000.0 ! convert m/day to mm/day Pr(:,:,it1:it2) = X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) !print*,Pr(:,:,it1:it2); pause write (Fname,"(A,I4.4,A)") trim(Inp_Dir)//"t2m_",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"t2m") Where (X3d==Mis_val) X3d = Miss Where (X3d/=Miss) X3d = X3d - 273.15 ! convert K to C Temp(:,:,it1:it2) = X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) !print*,Temp(:,:,it1:it2); pause write (Fname,"(A,I4.4,A)") trim(Inp_Dir)//"tx_",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"mx2t") Where (X3d==Mis_val) X3d = Miss Where (X3d/=Miss) X3d = X3d - 273.15 ! convert K to C Tmax(:,:,it1:it2) = X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) !print*,Tmax(:,:,it1:it2); pause write (Fname,"(A,I4.4,A)") trim(Inp_Dir)//"tm_",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"mn2t") Where (X3d==Mis_val) X3d = Miss Where (X3d/=Miss) X3d = X3d - 273.15 ! convert K to C Tmin(:,:,it1:it2) = X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) !print*,Tmin(:,:,it1:it2); pause Enddo Print*,"Reading data ... Done" End Subroutine Read_all_data !------------------------------------------------------------------------------ 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