Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, N365=365, Nreg=7 Integer, Parameter :: Max_evt=1000 ! Maximum number of events over all period Yr1-Yr2 Real, Allocatable :: Lon(:), Lat(:), xlon(:), xlat(:), time(:), xmod(:), Mask(:,:) Integer :: Nlon, Nlat, it, ilon, ilat !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:) :: Frq_obs, Dm_obs, Dx_obs ! Nreg,Nyear Real, Allocatable, Dimension(:,:) :: Sev_obs, HWday_obs Real, Allocatable, Dimension(:,:) :: Svx_obs Real, Allocatable, Dimension(:) :: Frq_obs_all, Dm_obs_all, Dx_obs_all ! Nreg Real, Allocatable, Dimension(:) :: Sev_obs_all, HWday_obs_all Real, Allocatable, Dimension(:) :: Svx_obs_all !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:) :: Reg_point ! Number of grid points in each subregion Real, Allocatable, Dimension(:,:) :: HWevt_obs ! Nreg,Ntime - for each subregion Real, Allocatable, Dimension(:,:) :: HWsev_reg ! Nreg,Ntime - mean values over subreg of HWsev_obs Real, Allocatable, Dimension(:,:,:) :: HWsev_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:) :: SPEI_obs ! Nlon,Nlat,Ntime !------------------------------------------------------------------------------ Character (Len=255) :: Fname,Fname1,Fname2,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, Ntime, Nyear, Mwd Real :: p, Dr_Thre Character (Len=6) :: Varnm ! tasmax Logical :: status, Empty, calc_qp !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, SPEI_out, Rain, T2m, Tx, Tm !------------------------------------------------------------------------------ Integer :: Nlo,Nla,Nt,Nz,imod !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Mask_File,Obs_Dir,Inp_Dir,Out_Dir Namelist /Data_Info/ Yr1, Yr2, Dr_Thre, p !------------------------------------------------------------------------------ Open (1, file="Namelist_CDH_Dr_HwTx_Reg.txt",Status="old") Read (1, File_Dir_Names) Read (1, Data_Info) Ntime=(Yr2-Yr1+1)*N365 Nyear=Yr2-Yr1+1 Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Mask,"reg") ! Land points Call Calc_reg_point ! Land point in each subregion !------------------------------------------------------------------------ if (Allocated(xlon)) Deallocate (xlon) if (Allocated(xlat)) Deallocate (xlat) if (Allocated(time)) Deallocate (time) !------------------------------------------------------------------------ Fsev_obs=trim(Out_Dir)//"HW_VN_sev_obs.nc" Call rd_ncdf_3d (Fsev_obs,xlon,xlat,time,Nlo,Nla,Nt,HWsev_obs,"sev") Deallocate (xlon,xlat,time) Fsev_obs=trim(Out_Dir)//"SPEI_obs.nc" Call rd_ncdf_3d (Fsev_obs,xlon,xlat,time,Nlo,Nla,Nt,SPEI_obs,"spei") Deallocate (xlon,xlat,time) Call Calc_HW_Evnts_obs_for_subreg Call Calc_Chars_obs Fname =trim(Out_Dir)//"CDH_HwTx_Evnt_VN_Reg.nc" Fname1=trim(Out_Dir)//"Chars_CDH_HwTx_yearly_obs_VN_Reg.nc" Fname2=trim(Out_Dir)//"Chars_CDH_HwTx_Clim_obs_VN_Reg.nc" Call Write_all_Chars_obs (Fname,Fname1,Fname2) print*,"Normal terminal of Chars_HW_Tx.f90" !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_reg_point Integer :: ilon,ilat,ic,ireg Allocate (Reg_point(Nreg)) Do ireg=1,Nreg ic=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Mask(ilon,ilat))==ireg) ic=ic+1 Enddo Enddo Reg_point(ireg)=ic Enddo End Subroutine Calc_reg_point !------------------------------------------------------------------------------ Subroutine Calc_Chars_obs Integer :: it,it1,it2,ilon,ilat,imod Real, Allocatable :: X1d(:), Y1d(:), M1(:), Spells(:), I_Start(:), I_End(:) Integer :: Nspls,Nspl3,istar,iend,Max_evt Real :: IMAX1, N_Spell_1, duration Real, dimension (Nyear) :: Frq, Dm, Dx, Nd, Sev, Svx Integer :: i,j,iyr Max_evt = Ntime Allocate (X1d(Ntime), Y1d(Ntime)) Allocate (Spells(Max_evt), I_Start(Max_evt), I_End(Max_evt)) Allocate (M1(Max_evt)) !--------------------------------------------------------------------- Allocate (Frq_obs(Nreg,Nyear),Dm_obs (Nreg,Nyear),Dx_obs(Nreg,Nyear)) Allocate (Sev_obs(Nreg,Nyear),HWday_obs(Nreg,Nyear)) Allocate (Svx_obs(Nreg,Nyear)) Frq_obs=Miss; Dm_obs=Miss; Dx_obs=Miss; Sev_obs=Miss; HWday_obs=Miss; Svx_obs=Miss; !--------------------------------------------------------------------- Allocate (Frq_obs_all(Nreg),Dm_obs_all (Nreg),Dx_obs_all (Nreg)) Allocate (Sev_obs_all(Nreg),Svx_obs_all(Nreg),HWday_obs_all(Nreg)) Frq_obs_all=Miss; Dm_obs_all=Miss; Dx_obs_all=Miss; Sev_obs_all=Miss; HWday_obs_all=Miss; Svx_obs_all=Miss; !--------------------------------------------------------------------- print*, "Calling Calc_Chars_obs..." ! Obs: !Do ilon=1,Nlon Do ireg=1,Nreg !if (Mask(ilon,ilat)/=Miss) then X1d(1:Ntime) = HWevt_obs(ireg,1:Ntime) ! Timeseries of Event days Y1d(1:Ntime) = HWsev_reg(ireg,1:Ntime) ! Timeseries of severity Call NSpell_and_Leng (X1d, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss) Nspls = N_Spell_1 ! So dot Frq=0; Dm=0; Dx=0; Nd=0; Sev=0; Svx=0 ! Freq, Mean Dur, Max Dur, NHWday, mean Sev Max Sev Do i=1,Nspls duration = Spells(i) ! I_End(i) - I_Start(i) + 1 if (duration >= 3.) then istar = Int(I_Start(i)) ! Thu tu Ngay bat mot HW iend = Int(I_End(i)) ! Thu tu Ngay ket thuc mot HW iyr = istar/N365 + 1 !mod(istar,N365) ! Thu tu Nam co HW Frq(iyr) = Frq(iyr) + 1 ! Number of HW event in year iyr Dm (iyr) = Dm (iyr) + duration ! Total duration of HW in year iyr Nd (iyr) = Dm (iyr) ! Number of HW days if (Dx(iyr)0) tmp = tmp + Y1d(j) ! Sum of positive difference between Tx and Qp Enddo Sev(iyr) = Sev(iyr) + tmp ! Total severity in year iyr if (Svx(iyr)0) then ! Trung binh nam Dm (iyr) = Dm (iyr)/Frq(iyr) Sev(iyr) = Sev(iyr)/Frq(iyr) endif Enddo Frq_obs(ireg,1:Nyear) = Frq(1:Nyear) ! Freq (evt) Dm_obs (ireg,1:Nyear) = Dm (1:Nyear) ! Mean duration (day) Dx_obs (ireg,1:Nyear) = Dx (1:Nyear) ! Max duration (day) HWday_obs(ireg,1:Nyear) = Nd(1:Nyear) ! Num of HW days (day) Sev_obs(ireg,1:Nyear) = Sev(1:Nyear) ! Mean sev (C) Svx_obs(ireg,1:Nyear) = Svx(1:Nyear) ! Max sev (C) ! Clim mean: Frq_obs_all(ireg) = Ave_No_Zero(Frq,Nyear,Miss) ! Clim mean Dm_obs_all (ireg) = Ave_No_Zero(Dm, Nyear,Miss) Dx_obs_all (ireg) = Ave_No_Zero(Dx, Nyear,Miss) Sev_obs_all(ireg) = Ave_No_Zero(Sev,Nyear,Miss) HWday_obs_all(ireg)= Ave_No_Zero(Nd, Nyear,Miss) Svx_obs_all(ireg) = Ave_No_Zero(Svx,Nyear,Miss) !endif Enddo !Enddo print*, "Calling Calc_Chars_obs Done" End Subroutine Calc_Chars_obs !------------------------------------------------------------------------------ Subroutine Calc_HW_Evnts_obs_for_subreg Integer :: imod,ilon,ilat,it,iday,iday1, ireg Real :: DrIdx, sev, ic, tmp, tmp1, sevreg if (Allocated (HWevt_obs)) Deallocate (HWevt_obs) if (Allocated (HWsev_reg)) Deallocate (HWsev_reg) Allocate (HWevt_obs(Nreg,Ntime), HWsev_reg(Nreg,Ntime)) HWevt_obs=Miss; HWsev_reg=Miss print*,"Calling Calc_HW_Evnts_obs_for_subreg ..." ! Obs: Do it=1,Ntime Do ireg=1,Nreg ic=0; sevreg=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Mask(ilon,ilat))==ireg) then ! if gridpoint is in the subregion ireg sev = HWsev_obs (ilon,ilat,it) ! HW severity = Tx - Qp DrIdx = SPEI_obs (ilon,ilat,it) if (sev>0.AND.DrIdx/=Miss.AND.DrIdx<=Dr_Thre) then ! HW occurs ic=ic+1; sevreg=sevreg-sev*DrIdx ! Total of product of (Tx - Qp)*(-DrIdx) endif endif Enddo Enddo tmp=ic/Reg_point(ireg); ! Ratio of area experienced CDH if (ic>0) then tmp1=sevreg/ic else tmp1=0 endif if (tmp>=p) then ! CDH occurs in ireg-th subreg at it if Ratio of its area exceeds p HWevt_obs(ireg,it)=1 HWsev_reg(ireg,it)=tmp1 else HWevt_obs(ireg,it)=0 HWsev_reg(ireg,it)=0 endif Enddo Enddo End Subroutine Calc_HW_Evnts_obs_for_subreg !------------------------------------------------------------------------------ Subroutine Write_all_Chars_obs (Fevnt,Fannual,Fclim) Character (Len=*) :: Fevnt,Fannual,Fclim Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str Integer :: ireg,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 :: xreg(Nreg) Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo if (Allocated(time)) Deallocate (time) ! Regional Event Allocate (time(Ntime)) Do it=1,Ntime; time(it)=it; Enddo lon_nm="reg"; lat_nm="day"; Var_ln="HW events for each climatic sub-regions of VN" Call wr_ncdf_2d (Fevnt,xreg,time,Nreg,Ntime,HWevt_obs,"evt", var_ln, "-", lon_nm, lat_nm, Miss) Var_ln="Mean severity of HW events for each climatic sub-regions of VN" Call wr_ncdf_2d_apd (Fevnt,HWsev_reg,"sev",Nreg,Ntime, var_ln,"C deg",lon_nm, lat_nm, Miss) Deallocate (time) Allocate (time(Nyear)) t_unit="years since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="year" Do i=1,Nyear; time(i)=i-1; Enddo ! Annual Obs: Var_ln="Annual number of the observed heatwave events (Frequency)" Call wr_ncdf_2d (Fannual,xreg,time,Nreg,Nyear,Frq_obs,"frq", var_ln, "-", lon_nm, lat_nm, Miss) Var_ln="Annual mean duration of the observed HW" Call wr_ncdf_2d_apd (Fannual,Dm_obs,"dm",Nreg,Nyear, var_ln,"day",lon_nm, lat_nm, Miss) Var_ln="Annual maximum duration of the observed HW" Call wr_ncdf_2d_apd (Fannual,Dx_obs,"dx",Nreg,Nyear, var_ln,"day",lon_nm, lat_nm, Miss) Var_ln="Annual mean severity of the observed HW" Call wr_ncdf_2d_apd (Fannual,Sev_obs,"sev",Nreg,Nyear, var_ln,"C deg",lon_nm, lat_nm, Miss) Var_ln="Annual maximum severity of the observed HW" Call wr_ncdf_2d_apd (Fannual,Svx_obs,"svx",Nreg,Nyear, var_ln,"C deg",lon_nm, lat_nm, Miss) Var_ln="Annual number of the observed HW days" Call wr_ncdf_2d_apd (Fannual,HWday_obs,"nday",Nreg,Nyear, var_ln,"day",lon_nm, lat_nm, Miss) ! Clim obs: Var_ln="Climatological mean of the observed heatwave events" Call wr_ncdf_1d (Fclim,xreg,Nreg,Frq_obs_all,"frq", var_ln,"-", lon_nm,Miss) Var_ln="Climatological mean of the observed mean HW duration" Call wr_ncdf_1d_apd (Fclim,Dm_obs_all,"dm",Nreg, var_ln,"day", lon_nm, flag) Var_ln="Climatological mean of the observed maximum HW duration" Call wr_ncdf_1d_apd (Fclim,Dx_obs_all,"dx",Nreg, var_ln,"day", lon_nm, flag) Var_ln="Climatological mean of the observed mean HW severity" Call wr_ncdf_1d_apd (Fclim,Sev_obs_all,"sev",Nreg, var_ln,"C", lon_nm, flag) Var_ln="Climatological mean of the observed maximum HW severity" Call wr_ncdf_1d_apd (Fclim,Svx_obs_all,"svx",Nreg, var_ln,"C", lon_nm, flag) Var_ln="Climatological mean of the observed number of HW days" Call wr_ncdf_1d_apd (Fclim,HWday_obs_all,"nday",Nreg, var_ln,"day", lon_nm, flag) End Subroutine Write_all_Chars_obs !------------------------------------------------------------------------------ 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_NSpell_and_Leng.f90" Include "/home/tanpv/Mylibs/Sub_Trend_MK.f90" Include "/home/tanpv/Mylibs/Ave_No_Zero.f90" !Include "/home/tanpv/Mylibs/Sub_SPEI_Daily_One_Stn_2024.f90" !Include "/home/tanpv/Mylibs/Sub_Daily_Percentile.f90" !-------------------------------------------------------------------- END ! Program