! Sarah E. Perkins, 2015: A review on the scientific understanding of heatwaves—Their measurement, ! driving mechanisms, and changes at the global scale. Atmospheric Research 164–165 (2015) 242–267. ! There are 4 discrete events (HWN); ! the length of the longest event is 10 days (HWD); ! the number of heatwave days in the sum of the duration of all four events (HWF); ! the average heatwave magnitude is the average temperature across all for events (HWM); ! and the heatwave amplitude is the hottest day of the event with the hottest average (HWA) Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Nreg=7, Mon1=4, Mon2=9 ! Start and End months of Summer season Integer :: Nst, Nyrs, Nmon, Nday, Yr1, Yr2 Character (Len=20), Allocatable :: Stn(:), Stname(:,:) Real, Allocatable :: Dat(:,:,:,:), Station(:), Year(:), Mon(:), Day(:), Lon(:), Lat(:), H(:), Reg(:) Integer :: ist, iyr, imon, iday, ilon, ilat, ireg Character (Len=100) :: Fname !--------------------------------------------------------------------- ! Heat wave: Real, Allocatable :: Tx(:,:,:,:), HD(:,:,:,:) ! Tx and Hotdays (Nst, Nyear, Ntime, Nreg); Ntime=1/4-30/9 Real, Allocatable :: Hevent(:,:,:) ! (Nyear, Ntime, Nreg) ! Characteristics: Real, Allocatable :: HWN(:,:) ! Freq (Nyear, Nreg) - Number of HW events Real, Allocatable :: HWDx(:,:) ! Max Duration: the length of the longest event (day) Real, Allocatable :: HWDm(:,:) ! Mean Duration (day) Real, Allocatable :: HDays(:,:) ! (Nyear, Nreg) - Number of hot days (day) Real, Allocatable :: HWF(:,:) ! (??) Tong so ngay trong tat ca cac dot HW (day) ! HWF: the number of heatwave days in the sum of the duration of all events Real, Allocatable :: HWM(:,:) ! (Nyear, Nreg) - Magnitude (C): Trung binh Tx trong tat ca cac ngay co HW !the average heatwave magnitude is the average temperature across all for events (HWM) Real, Allocatable :: HWA(:,:) ! (Nyear, Nreg) - Amplitude (C): Ngay nong nhat cua HW co Tx trung binh cao nhat !the heatwave amplitude is the hottest day of the event with the hottest average (HWA) !--------------------------------------------------------------------- ! Hot spell: !********************************************************************* Real, Allocatable :: q(:,:) ! Percentile (ist,ireg) Real, Allocatable :: X(:) ! Working array Integer :: Ntime, Nst1, Nyr1, it, DayofYear, iyr1 Integer :: Nstn(Nreg) Real, Parameter :: p=0.95 ! probability Real, Parameter :: HW_len_min=3 ! It nhat 3 ngay lien tiep co Hevent=1 Real, Parameter :: H_Area=0.20 ! Phan tram so tram co Tx>Threshold !--------------------------------------------------------------------- Fname="/home/tanpv/OBS_Stn/Daily_NC/Tx_Daily.nc" Call Read_Tx(Fname) Yr1=1980; Yr2=2014; Nyr1=(Yr2-Yr1+1); DayofYear=0 Do imon=Mon1, Mon2 ! Tong so ngay tu 1/4-30/9 DayofYear=DayofYear+DayOfMon(imon,2019) Enddo Ntime=DayofYear ! Tong so ngay cac thang tu Mon1-Mon2 cua cac nam tu Yr1 den Yr2 Allocate (Tx(Nst,Nyr1,DayofYear,Nreg),Stname(Nst,Nreg),HD(Nst,Nyr1,DayofYear,Nreg)) Call Extract_Tx print*,"Ntime=",Ntime Allocate (q(Nst,Nreg),X(Nyr1*Ntime),Hevent(Nyr1,Ntime,Nreg)) Call Calc_percentile_and_HD (p) !Do ist=1, Nst ! print*,Stname(ist), q(ist) ! Do iyr=1, Nyr1 !print*,iyr+Yr1-1 !print "(60I2)", (Int(HD(ist,iyr,:))) ! Enddo !pause !Enddo Call Calc_Hot_Event !Do ireg=1, Nreg !print*, "Subregion: ", ireg !Do iyr=1, Nyr1 ! print*,iyr+Yr1-1 ! print "(60I2)", (Int(Hevent(iyr,:,ireg))) !ENddo !pause !Enddo Call Calc_HW_Char Fname="Heatwave_Chars.nc" Call Write_Char_to_NC (Fname) CONTAINS !--------------------------------------------------------------------- Subroutine Write_Char_to_NC (Fname) Integer :: iyr, ireg Real, Allocatable :: Years(:), Regs(:) Character (Len=*) :: Fname Character (Len=50) :: var_ln Allocate (Years(Nyr1), Regs(Nreg)) Do iyr=1, Nyr1; Years(iyr)=iyr+Yr1-1; Enddo Do ireg=1, Nreg; Regs(ireg)=ireg; Enddo var_ln="Number of HW events" Call wr_ncdf_2d (Fname,Years,Regs,Nyr1,Nreg, HWN, "nhw", var_ln, "-", "year", "reg", Miss) var_ln="Mean Duration of HW events" Call wr_ncdf_2d_apd (Fname, HWDm, "dmean", Nyr1,Nreg, var_ln, "day", "year", "reg", Miss) var_ln="Maximum Duration of HW events" Call wr_ncdf_2d_apd (Fname, HWDx, "dmx", Nyr1,Nreg, var_ln, "day", "year", "reg", Miss) var_ln="Number of HW days" Call wr_ncdf_2d_apd (Fname, HDays, "nday", Nyr1,Nreg, var_ln, "days/year", "year", "reg", Miss) var_ln="HW magnitude (average temperature across all for events)" Call wr_ncdf_2d_apd (Fname, HWM, "hwmag", Nyr1,Nreg, var_ln, "C-deg", "year", "reg", Miss) var_ln="HW amplitude (average temperature of hottest days across all events)" Call wr_ncdf_2d_apd (Fname, HWA, "hwamp", Nyr1,Nreg, var_ln, "C-deg", "year", "reg", Miss) End Subroutine Write_Char_to_NC !--------------------------------------------------------------------- Subroutine Calc_HW_Char Integer :: ireg, ist, iyr, it, ispl Real, Allocatable :: Spells(:), X(:), SDate(:), EDate(:) Real :: N_Spell, Imax, Nhw, Dmx, Dmean, Hdy, Txx, Tx0, Txm, NTx ! N_Spell: So dot trong toan chuoi ! Spells(N_Spell): Cac phan tu la Do dai tung dot Allocate (HWN(Nyr1,Nreg),HWDx(Nyr1,Nreg),HWDm(Nyr1,Nreg),HDays(Nyr1,Nreg),HWM(Nyr1,Nreg),HWA(Nyr1,Nreg)) Allocate (Spells(Ntime), X(Ntime), SDate(Ntime), EDate(Ntime)) !Real, Allocatable :: HWN(Nyr1,Nreg) ! Freq (Nyear, Nreg) - Number of HW events !Real, Allocatable :: HWDx(:,:) ! Max Duration: the length of the longest event (day) !Real, Allocatable :: HWDm(:,:) ! Mean Duration (day) !Real, Allocatable :: HDays(:,:) ! (Nyear, Nreg) - Number of hot days (day) !the number of heatwave days in the sum of the duration of all events !Real, Allocatable :: HWM(:,:) ! (Nyear, Nreg) - Magnitude (C) !the average heatwave magnitude is the average temperature across all for events (HWM) !Real, Allocatable :: HWA(:,:) ! (Nyear, Nreg) - Amplitude (C) !the heatwave amplitude is the hottest day of the event with the hottest average (HWA) ! Hevent(iyr,it,ireg) ! Hot event occurs over a region Do ireg=1, Nreg Do iyr=1, Nyr1 X(1:Ntime)=Hevent(iyr,1:Ntime,ireg) Call ESpells (X, Ntime, Imax, Spells, N_Spell, SDate, EDate, Miss) !print*,"reg,year=",ireg,iyr+Yr1-1 Nhw=0; Dmx=0; Dmean=0; Hdy=0; Txx=0; Txm=0; Ntx=0 Do it=1, Int(N_Spell) ! For each spell if (Spells(it)>=HW_len_min) then ! if spell is HW (Spells(it) is num of days of it-th event) Nhw=Nhw+1 ! Number of HW events Hdy=Hdy+Spells(it) ! Number of HW days Dmean=Dmean+Spells(it) ! Number of HW days (duration) if (Spells(it)>Dmx) Dmx=Spells(it) ! Maximum duration Tx0=0 Do ispl=Int(SDate(it)), Int(Edate(it)) ! Cac ngay trong 1 HW Do ist=1, Nstn(ireg) ! For each station in the region if (Tx(ist,iyr,ispl,ireg)>=q(ist,ireg)) then ! Hot day Txm=Txm+Tx(ist,iyr,ispl,ireg) ! Sum of Temp of hot days Ntx=Ntx+1 ! Num of hot day if (Tx(ist,iyr,ispl,ireg)>Tx0) Tx0=Tx(ist,iyr,ispl,ireg) ! Max Tx of the region in each HW event !print*," Tx0: ", Tx0 endif Enddo Enddo !print*," Nhw, Tx0: ",Nhw, Tx0 Txx=Txx+Tx0 ! Sum of Max Tx of all HW events endif ! print*,it, Spells(it) Enddo if (Nhw>0) Dmean=Dmean/Nhw ! Do keo dai trung binh cua cac HWs if (Ntx>0) Txm=Txm/Ntx ! Tx trung binh cua cac tram trong tat ca cac HW events if (Nhw>0) Txx=Txx/Nhw ! Mean of Max Tx over all HW events print "(2I5,A,10F5.1)",ireg,iyr+Yr1-1," Nhw,Dmean,Dmx=",Nhw,Dmean,Dmx,Hdy, Txm, Txx HWN(iyr,ireg) = Nhw ! So HW/year HWDm(iyr,ireg) = Dmean ! HW mean duration HWDx(iyr,ireg) = Dmx ! HW maximum duration HDays(iyr,ireg)= Hdy ! Total number of HW days HWM(iyr,ireg) = Txm ! HW magnitude (trung binh cac Tx>q cua tat ca cac tram) HWA(iyr,ireg) = Txx ! HW amplitude (trung binh Txx cac HW) Enddo pause Enddo End Subroutine Calc_HW_Char !--------------------------------------------------------------------- Subroutine Calc_Hot_Event Integer :: ist, it, N, ireg, NHSt Real :: X(Nst), Imax1 Hevent(:,:,:)=0 Do ireg=1, Nreg NHSt=Int(Nstn(ireg)*H_Area) ! So tram trong vung can thoa man dieu kien cua HW-day print*,"Num Stn HD,Ntime: ",NHSt,Ntime Do iyr=1, Nyr1 Do it=1, Ntime N=0 Do ist=1, Nstn(ireg) if (HD(ist,iyr,it,ireg)==1) N=N+1 Enddo if (N>=NHSt) Hevent(iyr,it,ireg)=1 ! Hot event occurs over a region Enddo Enddo Enddo End Subroutine Calc_Hot_Event !--------------------------------------------------------------------- Subroutine Calc_percentile_and_HD (p) Integer :: ireg, ist, iyr, it, N !--------------------- HD(:,:,:,:) = 0 Do ireg=1, Nreg Do ist=1, Nstn(ireg) it=0 Do iyr=1, Nyr1 Do iday=1, DayofYear if (Tx(ist,iyr,iday,ireg)/=Miss) then it=it+1; X(it)=Tx(ist,iyr,iday,ireg) endif Enddo Enddo !print*,X(1:it) N=it if (.not.Empty_Ar1(X,N,Miss).AND.N>Ntime/2) then q(ist,ireg) = percentile(X,N,p, Miss) ! q=X[F(x)=p] endif print*,"Stn,q=",Stname(ist,ireg),q(ist,ireg) Do iyr=1, Nyr1 Do iday=1, DayofYear if (Tx(ist,iyr,iday,ireg)>=q(ist,ireg)) HD(ist,iyr,iday,ireg)=1 Enddo Enddo Enddo Enddo End Subroutine Calc_percentile_and_HD !---------------------------------------------------------------------- Subroutine Extract_Tx Integer :: ireg, ist, ist1, imon, iday, iyr, iyr1 Do ireg=1, Nreg ist1=0 Do ist=1, Nst if (reg(ist)==ireg) then iyr1=0; ist1=ist1+1; Stname(ist1,ireg)=Stn(ist) Do iyr=1, Nyrs if (Int(Year(iyr))>=Yr1.AND.Int(Year(iyr))<=Yr2) then iyr1=iyr1+1; it=0 Do imon=Mon1, Mon2 ! Apr - Sep Nday=DayOfMon(imon, Int(Year(iyr))) Do iday=1, Nday it=it+1 Tx(ist1,iyr1,it,ireg)=dat(ist,iyr,imon,iday) Enddo Enddo !Ntime=it !print*,Stn(ist), Year(iyr) !print "(20F6.1)", Tx(nst1,iyr1,:) endif !pause Enddo endif Enddo Nstn(ireg)=ist1 Enddo End Subroutine Extract_Tx !---------------------------------------------------------------------- Subroutine Read_Tx (Fname) Character (Len=*) :: Fname Call rd_ncdf_4d (Fname,Station,Year,Mon,Day,Nst, Nyrs, Nmon, Nday,Dat,"Tx") Deallocate(Station) Call rd_ncdf_1d (Fname,Station, Nst, Lon, "Lon") Deallocate(Station) Call rd_ncdf_1d (Fname,Station, Nst, Lat, "Lat") Deallocate(Station) Call rd_ncdf_1d (Fname,Station, Nst, H, "h") Deallocate(Station) Call rd_ncdf_1d (Fname,Station, Nst, reg, "reg") Deallocate(Station) Call rd_ncdf_1d_Char (Fname,Station, Nst, Stn, "Stname") End Subroutine Read_Tx !----------------------------------------------------------------------- Logical Function Empty_Ar1 (X, N, amssng) ! Tra ve .True. neu X(N) khuyet toan bo SL Integer :: N, i, k Real :: X(N), amssng k = 0 ! Dem so phan tu khuyet cua X Do i=1, N if (X(i) == amssng) k = k + 1 Enddo if (k > N/10) then ; Empty_Ar1 = .True. ; else ; Empty_Ar1 = .False. ; endif End Function Empty_Ar1 !----------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/DayOfMon.f90" Include "/home/tanpv/Mylibs/Sub_Spells.f90" Include "/home/tanpv/Mylibs/Sub_Percentile.f90" !----------------------------------------------------------------------- END