! 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, Parameter :: Ybl1=1981, Ybl2=2010 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) :: DirOut, Fname,Fname1,Fname2,Fname3 !--------------------------------------------------------------------- ! Heat wave: Real, Allocatable :: Tx0(:,:,:), Tx(:,:,:), HD(:,:,:) ! Tx and Hotdays (Nst, Nyear, Ntime, Nreg); Ntime=1/4-30/9 Real, Allocatable :: Tx01(:,:,:) 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) Real, Allocatable :: HWS(:,:) ! (Nyear, Nst) - Mean HW Severity - Sum of Tx-q Real, Allocatable :: HWSx(:,:) ! (Nyear, Nst) - Max HW Severity Real, Allocatable, Dimension (:,:) :: HWN_Reg, HWDx_Reg, HWDm_Reg, HDays_Reg,HWF_Reg, HWM_Reg, HWA_Reg,HWS_Reg,HWSx_Reg ! (Nyear,Nreg) !--------------------------------------------------------------------- ! Trens Analusis: Real, Allocatable, Dimension (:,:) :: HWN_Sen,HWDx_Sen,HWDm_Sen,HDays_Sen,HWF_Sen,HWM_Sen,HWA_Sen,HWS_Sen,HWSx_Sen ! (Nst,2) Real, Allocatable, Dimension (:,:,:) :: HW_Reg_Sen ! (Nreg=7,NChar=9,2=Sen+Prob) !********************************************************************* Real, Allocatable :: q(:) ! Percentile (ist,ireg) Real, Allocatable :: qw(:,:) ! Percentile for daily using window data (Nst,Ntime), Ntime=1/4-30/9 Real, Allocatable :: q25(:), q75(:) ! (Nst) Annual Tx percentiles for BL period at each station Integer, Parameter :: Mwd = 31, Mwd2 = 15 !Integer, Parameter :: Mwd = 15, Mwd2 = 7 Real :: NSt_reg(Nreg) ! So tram trong tung vung Real, Allocatable :: X(:) ! Working array Integer :: Ntime0,Ntime, Nst1, Nyr1, it, DayofYear, iyr1 Integer :: Nstn(Nreg) Real, Parameter :: p=0.90 ! 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 Integer, Parameter :: NChar=9 ! So characteristics !--------------------------------------------------------------------- Fname="/home/tanpv/Data/OBS/OBS_Stn/Daily/NC/Tx_Daily.nc" Call Read_Tx(Fname) Yr1=1980; Yr2=2018; Nyr1=(Yr2-Yr1+1); DayofYear=0 NyBL=(Ybl2-Ybl1+1) Do imon=Mon1, Mon2 ! Tong so ngay tu 1/4-30/9, khong co thang 2 nen khong co nam nhuan DayofYear=DayofYear+DayOfMon(imon,2019) Enddo Ntime0=0 Do imon=Mon1-1, Mon2+1 ! Tong so ngay tu 1/3-31/10, khong co thang 2 nen khong co nam nhuan Ntime0=Ntime0+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,Ntime),HD(Nst,Nyr1,Ntime)) Allocate (Tx0(Nst,NyBL,Ntime0),Tx01(Nst,NyBL,Ntime)) Call Extract_Tx Allocate (qw(Nst,Ntime),X(NyBL*Mwd)) !,Hevent(Nyr1,Ntime)) Allocate (q25(Nst),q75(Nst)) ! IQR print*,"Nst,Nyr1,NyBL,Ntime,Ntime0=",Nst,Nyr1,NyBL,Ntime,Ntime0 Call Calc_percentile_and_HD (p) !print*,"Nst,Nyr1,Ntime=",Nst,Nyr1,Ntime; pause DirOut="/home/tanpv/Progs/HW/RES/" Call Calc_HW_Char Call Calc_Trend Call Write_GrADS_Stn ("Sensl",1,"sensl") ! (Grads_File,ivar,varn) Call Write_GrADS_Stn ("prob", 2,"prob") Fname= trim(DirOut)//"Heatwave_Chars_Stn.nc"; Fname1=trim(DirOut)//"Percentiles.nc"; Fname2=trim(DirOut)//"Heatwave_Chars_Reg.nc" Fname3=trim(DirOut)//"Trend_Char_Reg.nc" Call Write_Char_to_NC (Fname,Fname1,Fname2,Fname3) print"(7I5)", (Int(Nst_reg(ireg)),ireg=1,Nreg) CONTAINS !--------------------------------------------------------------------- Subroutine Write_Char_to_NC (Fname,F_qw,F_reg,F_reg_trend) Integer :: iyr, it, ist, ireg Real, Allocatable :: Years(:), StnID(:),X(:),Regn(:),XX(:,:),xchar(:) Character (Len=*) :: Fname,F_qw,F_reg,F_reg_trend Character (Len=250) :: var_ln Allocate (Years(Nyr1), StnID(Nst), X(Nst), Regn(Nreg), xchar(Nchar)) Do iyr=1, Nyr1; Years(iyr)=iyr+Yr1-1; Enddo Do ist=1, Nst; StnID(ist)=ist; Enddo Do ireg=1,Nreg; Regn(ireg)=ireg; Enddo Do it=1,Nchar; xchar(it)=it; Enddo ! Station's Chars: var_ln="Number of HW events" Call wr_ncdf_2d (Fname,StnID,Years,Nst,Nyr1, HWN, "nhw", var_ln, "-", "stns","year", Miss) var_ln="Mean Duration of HW events" Call wr_ncdf_2d_apd (Fname, HWDm, "dmean", Nst,Nyr1, var_ln, "day", "stns","year", Miss) var_ln="Maximum Duration of HW events" Call wr_ncdf_2d_apd (Fname, HWDx, "dmx", Nst,Nyr1, var_ln, "day", "stns","year", Miss) var_ln="Number of Hot days" Call wr_ncdf_2d_apd (Fname, HDays, "hotday", Nst,Nyr1, var_ln, "days/year", "stns","year", Miss) var_ln="Number of HW days" Call wr_ncdf_2d_apd (Fname, HWF, "hwday", Nst,Nyr1, var_ln, "days/year", "stns","year", Miss) var_ln="HW magnitude (average temperature across all for events)" Call wr_ncdf_2d_apd (Fname, HWM, "hwmag", Nst,Nyr1, var_ln, "C-deg", "stns","year", Miss) var_ln="HW amplitude (average temperature of hottest days across all events)" Call wr_ncdf_2d_apd (Fname, HWA, "hwamp", Nst,Nyr1, var_ln, "C-deg", "stns","year", Miss) var_ln="Mean HW severity - defined as sum of differences between Tx and Percentile used as HW threshold" Call wr_ncdf_2d_apd (Fname, HWS, "hws", Nst,Nyr1, var_ln, "C-deg", "stns","year", Miss) var_ln="Max HW severity - defined as sum of differences between Tx and Percentile used as HW threshold" Call wr_ncdf_2d_apd (Fname, HWSx, "hwsx", Nst,Nyr1, var_ln, "C-deg", "stns","year", Miss) Call wr_ncdf_1d_apd (Fname,Lon,"lon", Nst, "Station's Longtitude", "DegE", "stns", Miss) Call wr_ncdf_1d_apd (Fname,Lat,"lat", Nst, "Station's latitude", "DegN", "stns", Miss) ! Station's Trends: X(1:Nst)=HWN_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwn_sl", Nst, "Sen slopes of HWN","hw/yr", "stns", Miss) X(1:Nst)=HWN_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwn_pb", Nst, "P-values of HWN", "-", "stns", Miss) X(1:Nst)=HWDx_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwdx_sl", Nst, "Max duration of HWs","day", "stns", Miss) X(1:Nst)=HWDx_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwdx_pb", Nst, "P-values of HWDx", "-", "stns", Miss) X(1:Nst)=HWDm_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwdm_sl", Nst, "Mean duration of HWs","day", "stns", Miss) X(1:Nst)=HWDm_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwdm_pb", Nst, "P-values of HWDm", "-", "stns", Miss) X(1:Nst)=HDays_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hdays_sl", Nst, "Number of Hot days","day", "stns", Miss) X(1:Nst)=HDays_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hdays_pb", Nst, "P-values of HDays", "-", "stns", Miss) X(1:Nst)=HWF_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwf_sl", Nst, "Number of HW days (freq)","day", "stns", Miss) X(1:Nst)=HWF_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwf_pb", Nst, "P-values of HWF", "-", "stns", Miss) X(1:Nst)=HWM_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwm_sl", Nst, "HW Magnitude","C", "stns", Miss) X(1:Nst)=HWM_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwm_pb", Nst, "P-values of HWM", "-", "stns", Miss) X(1:Nst)=HWA_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwa_sl", Nst, "HW Amplitude","C", "stns", Miss) X(1:Nst)=HWA_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwa_pb", Nst, "P-values of HWA", "-", "stns", Miss) X(1:Nst)=HWS_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hws_sl", Nst, "Mean HW severity","C", "stns", Miss) X(1:Nst)=HWS_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hws_pb", Nst, "P-values of HWS", "-", "stns", Miss) X(1:Nst)=HWSx_Sen(1:Nst,1) Call wr_ncdf_1d_apd (Fname,X,"hwsx_sl", Nst, "Max HW severity","C", "stns", Miss) X(1:Nst)=HWSx_Sen(1:Nst,2) Call wr_ncdf_1d_apd (Fname,X,"hwsx_pb", Nst, "P-values of HWSx", "-", "stns", Miss) Call wr_ncdf_1d_apd_char (Fname, Stn,"stname", Nst, "Station's Name", "-", "stns", Miss) ! Sub-regions' Trends: Allocate (XX(Nreg,Nchar)) var_ln="Trends of change in HW characteristics in the subregions: 1-HWN, 2-HWDx, 3-HWDm, 4-HWF, 5-HWM, 6-HWA, 7-HDays, 8-HWS, 9-HWSx" Do ireg=1,Nreg; XX(ireg,1:Nchar)=HW_Reg_Sen(ireg,1:Nchar,1); Enddo Call wr_ncdf_2d (F_reg_trend,Regn,xchar,Nreg,Nchar,XX,"sensl",var_ln,"1/year","reg","char",Miss) Do ireg=1,Nreg; XX(ireg,1:Nchar)=HW_Reg_Sen(ireg,1:Nchar,2); Enddo var_ln="Probability - P values" Call wr_ncdf_2d_apd (F_reg_trend,XX,"prob",Nreg,Nchar,var_ln,"-","reg","char",Miss) ! Region's Chars: var_ln="Number of HW events" Call wr_ncdf_2d (F_reg,Years,Regn,Nyr1,Nreg, HWN_Reg, "nhw", var_ln, "-", "year","reg", Miss) var_ln="Mean Duration of HW events" Call wr_ncdf_2d_apd (F_reg, HWDm_Reg, "dmean", Nyr1,Nreg, var_ln, "day", "year","reg", Miss) var_ln="Maximum Duration of HW events" Call wr_ncdf_2d_apd (F_reg, HWDx_Reg, "dmx", Nyr1,Nreg, var_ln, "day", "year","reg", Miss) var_ln="Number of HW days" Call wr_ncdf_2d_apd (F_reg, HWF_Reg, "hwdays", Nyr1,Nreg, var_ln, "day", "year","reg", Miss) var_ln="Number of Hot days" Call wr_ncdf_2d_apd (F_reg, HDays_Reg,"hdays", Nyr1,Nreg, var_ln, "day", "year","reg", Miss) var_ln="HW magnitude (average temperature across all for events)" Call wr_ncdf_2d_apd (F_reg, HWM_Reg, "hwmag", Nyr1,Nreg, var_ln, "C", "year","reg", Miss) var_ln="HW amplitude (average temperature of hottest days across all events)" Call wr_ncdf_2d_apd (F_reg, HWA_Reg, "hwamp", Nyr1,Nreg, var_ln, "C", "year","reg", Miss) var_ln="Mean HW severity - defined as sum of differences between Tx and Percentile used as HW threshold" Call wr_ncdf_2d_apd (F_reg, HWS_Reg, "hws", Nyr1,Nreg, var_ln, "C", "year","reg", Miss) var_ln="Maximum HW severity - defined as sum of differences between Tx and Percentile used as HW threshold" Call wr_ncdf_2d_apd (F_reg, HWSx_Reg, "hwsx", Nyr1,Nreg, var_ln, "C", "year","reg", Miss) ! Station's Percentiles: Deallocate (Years,XX); Allocate (Years(Ntime),XX(Ntime,Nst)) Do it=1,Ntime; Years(it)=it; Enddo ! Days since 1/4-30/9 Do ist=1,Nst; Do it=1,Ntime; XX(it,ist)=qw(ist,it); Enddo; Enddo var_ln="Percentiles of daily Tx for each calendar day using a centred 31-day-long time window of Baseline (1981-2010)" Call wr_ncdf_2d (F_qw,Years,StnID,Ntime,Nst, XX, "q", var_ln, "C", "day","stn", Miss) End Subroutine Write_Char_to_NC !--------------------------------------------------------------------- Subroutine Calc_HW_Char Integer :: ireg, ist, iyr, it, ispl, n1,n2,n3,n4,n5,n6,n7,n8,n9 Real, Allocatable :: Spells(:), X(:), SDate(:), EDate(:) Real :: N_Spell, Imax, Nhw, Dmx, Dmean, Hdy, Hdy1, Txx, Tx0, Txm, NTx, Hse,Hsex,Hse0 Allocate (HWN(Nst,Nyr1),HWDx(Nst,Nyr1),HWDm(Nst,Nyr1),HDays(Nst,Nyr1)) Allocate (HWM(Nst,Nyr1),HWA(Nst,Nyr1),HWF(Nst,Nyr1),HWS(Nst,Nyr1),HWSx(Nst,Nyr1)) Allocate (Spells(Ntime), X(Ntime), SDate(Ntime), EDate(Ntime)) HWN(:,:)=Miss; HWDm(:,:)=Miss; HWDx(:,:)=Miss; HWS(:,:)=Miss HDays(:,:)=Miss; HWM(:,:)=Miss; HWA(:,:)=Miss; HDays(:,:)=Miss Do ist=1, Nst Do iyr=1, Nyr1 X(1:Ntime)=HD(ist,iyr,1:Ntime) !Hevent(iyr,1:Ntime,ireg) Call ESpells (X, Ntime, Imax, Spells, N_Spell, SDate, EDate, Miss) Nhw=0; Dmx=0; Dmean=0; Hdy=0; Txx=0; Txm=0; Ntx=0; Hse=0; Hsex=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; Hse0=0 Do ispl=Int(SDate(it)), Int(Edate(it)) ! Cac ngay trong 1 HW if (Tx(ist,iyr,ispl)/=Miss.AND. & qw(ist,ispl)/=Miss.AND. & Tx(ist,iyr,ispl)>=qw(ist,ispl)) then ! HW day Hse=Hse+(Tx(ist,iyr,ispl)-qw(ist,ispl)) ! HW severity = Sum of Tx-q Hse0=Hse0+(Tx(ist,iyr,ispl)-qw(ist,ispl)) ! Find the Max of HW severity Txm=Txm+Tx(ist,iyr,ispl) ! Sum of Temp of HW days Ntx=Ntx+1 ! Num of HW day if (Tx(ist,iyr,ispl)>Tx0) Tx0=Tx(ist,iyr,ispl) ! Max Tx of the region in each HW event endif Enddo Txx=Txx+Tx0 ! Sum of Max Tx of all HW events if (Hsex0) 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 if (Nhw>0) Hse=Hse/Nhw ! Mean Severity HWN(ist,iyr) = Nhw ! So HW/year HWDm(ist,iyr) = Dmean ! HW mean duration HWDx(ist,iyr) = Dmx ! HW maximum duration HWF(ist,iyr) = Hdy ! Total number of HW days HWM(ist,iyr) = Txm ! HW magnitude (trung binh cac Tx>q cua tat ca cac tram) HWA(ist,iyr) = Txx ! HW amplitude (trung binh Txx cac HW) HWS(ist,iyr) = Hse ! HW Mean severity HWSx(ist,iyr) = Hsex ! Max HW severity Hdy1=0 ! So ngay nong Do it=1,Ntime if (X(it)==1) Hdy1=Hdy1+1 ! Num of Hot days Enddo HDays(ist,iyr)=Hdy1 Enddo Enddo ! Reginal mean: Allocate (HWN_Reg(Nyr1,Nreg), HWDx_Reg(Nyr1,Nreg), HWDm_Reg(Nyr1,Nreg),HWS_Reg(Nyr1,Nreg)) Allocate (HDays_Reg(Nyr1,Nreg),HWF_Reg(Nyr1,Nreg), HWM_Reg(Nyr1,Nreg), HWA_Reg(Nyr1,Nreg)) ! (Nyear,Nreg) Allocate (HWSx_Reg(Nyr1,Nreg)) HWN_Reg=Miss; HWDx_Reg=Miss; HWDm_Reg=Miss; HDays_Reg=Miss; HWF_Reg=Miss; HWM_Reg=Miss; HWA_Reg=Miss; HWS_Reg=Miss; HWSx_Reg=Miss Do iyr=1,Nyr1 Do ireg=1,Nreg Nhw=0; Dmean=0; Dmx=0; Hdy=0;Txm=0; Txx=0; Hdy1=0; Hse=0; Hsex=0 n1=0; n2=0; n3=0; n4=0; n5=0; n6=0; n7=0; n8=0; n9=0 Do ist=1,Nst if (Reg(ist)==ireg) then if (HWN(ist,iyr)>0) then n1=n1+1; Nhw=Nhw+HWN(ist,iyr) endif if (HWDx(ist,iyr)>0) then n2=n2+1; Dmx=Dmx+HWDx(ist,iyr) endif if (HWDm(ist,iyr)>0) then n3=n3+1; Dmean=Dmean+HWDm(ist,iyr) endif if (HWF(ist,iyr)>0) then n4=n4+1; Hdy=Hdy+HWF(ist,iyr) endif if (HWM(ist,iyr)>0) then n5=n5+1; Txm=Txm+HWM(ist,iyr) endif if (HWA(ist,iyr)>0) then n6=n6+1; Txx=Txx+HWA(ist,iyr) endif if (HDays(ist,iyr)>0) then n7=n7+1; Hdy1=Hdy1+HDays(ist,iyr) endif if (HWS(ist,iyr)>0) then n8=n8+1; Hse=Hse+HWS(ist,iyr) endif if (HWSx(ist,iyr)>0) then n9=n9+1; Hsex=Hsex+HWSx(ist,iyr) endif endif Enddo if (n1>0) HWN_Reg (iyr,ireg)=Nhw/n1 if (n2>0) HWDx_Reg (iyr,ireg)=Dmx/n2 if (n3>0) HWDm_Reg (iyr,ireg)=Dmean/n3 if (n4>0) HWF_Reg (iyr,ireg)=Hdy/n4 if (n5>0) HWM_Reg (iyr,ireg)=Txm/n5 if (n6>0) HWA_Reg (iyr,ireg)=Txx/n6 if (n7>0) HDays_Reg(iyr,ireg)=Hdy1/n7 if (n8>0) HWS_Reg (iyr,ireg)=Hse/n8 if (n9>0) HWSx_Reg (iyr,ireg)=Hsex/n9 Enddo Enddo print*,"Calc_HW_Char...Done" End Subroutine Calc_HW_Char !--------------------------------------------------------------------- Subroutine Calc_Trend Integer :: ist,iyr,it, N0, N1 Real :: X(Nyr1), SenSlp, Prob Allocate (HWN_Sen(Nst,2),HWDx_Sen(Nst,2),HWDm_Sen(Nst,2),HDays_Sen(Nst,2)) Allocate (HWF_Sen(Nst,2),HWM_Sen(Nst,2),HWA_Sen(Nst,2),HWS_Sen(Nst,2)) ! (Nst,2) Allocate (HWSx_Sen(Nst,2)) N0=Nyr1/2 ! Dieu kien de tinh SenSLop Do ist=1,Nst X(1:Nyr1)=HWN(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWN_Sen(ist,1)=SenSlp; HWN_Sen(ist,2)=prob ! 1 X(1:Nyr1)=HWDx(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWDx_Sen(ist,1)=SenSlp; HWDx_Sen(ist,2)=prob ! 2 X(1:Nyr1)=HWDm(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWDm_Sen(ist,1)=SenSlp; HWDm_Sen(ist,2)=prob ! 3 X(1:Nyr1)=HDays(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HDays_Sen(ist,1)=SenSlp; HDays_Sen(ist,2)=prob ! 4 X(1:Nyr1)=HWF(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWF_Sen(ist,1)=SenSlp; HWF_Sen(ist,2)=prob ! 5 X(1:Nyr1)=HWM(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWM_Sen(ist,1)=SenSlp; HWM_Sen(ist,2)=prob ! 6 X(1:Nyr1)=HWA(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWA_Sen(ist,1)=SenSlp; HWA_Sen(ist,2)=prob ! 7 X(1:Nyr1)=HWS(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWS_Sen(ist,1)=SenSlp; HWS_Sen(ist,2)=prob ! 8 X(1:Nyr1)=HWSx(ist,1:Nyr1) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HWSx_Sen(ist,1)=SenSlp; HWSx_Sen(ist,2)=prob ! 9 Enddo ! For sub-regions: Allocate (HW_Reg_Sen(Nreg,Nchar,2)) Do ireg=1,Nreg X(1:Nyr1)=HWN_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,1,1)=SenSlp; HW_Reg_Sen(ireg,1,2)=prob X(1:Nyr1)=HWDx_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,2,1)=SenSlp; HW_Reg_Sen(ireg,2,2)=prob X(1:Nyr1)=HWDm_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,3,1)=SenSlp; HW_Reg_Sen(ireg,3,2)=prob X(1:Nyr1)=HWF_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,4,1)=SenSlp; HW_Reg_Sen(ireg,4,2)=prob X(1:Nyr1)=HWM_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,5,1)=SenSlp; HW_Reg_Sen(ireg,5,2)=prob X(1:Nyr1)=HWA_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,6,1)=SenSlp; HW_Reg_Sen(ireg,6,2)=prob X(1:Nyr1)=HDays_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,7,1)=SenSlp; HW_Reg_Sen(ireg,7,2)=prob X(1:Nyr1)=HWS_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,8,1)=SenSlp; HW_Reg_Sen(ireg,8,2)=prob X(1:Nyr1)=HWSx_Reg(1:Nyr1,ireg) Call MK_Sen(X,Nyr1,N0,Miss,SenSlp,prob,N1); HW_Reg_Sen(ireg,9,1)=SenSlp; HW_Reg_Sen(ireg,9,2)=prob Enddo !Open (3,file=trim(DirOut)//"Trend.txt",status="Unknown") !write (3,"(1x,A3,A15,8(4x,A12,4x))") "No","Stn","HWN","HWDx","HWDm","HDays","HWF","HWM","HWA","HWS" !Do ist=1,Nst !if (q(ist)/=Miss) & ! write (3,"(I3,1x,A15,10(F10.6,F10.3))") ist,Stn(ist),(HWN_Sen(ist,i), i=1,2),(HWDx_Sen(ist,i), i=1,2),(HWDm_Sen(ist,i), i=1,2), & ! (HDays_Sen(ist,i), i=1,2),(HWF_Sen(ist,i), i=1,2),(HWM_Sen(ist,i), i=1,2),(HWA_Sen(ist,i), i=1,2),(HWS_Sen(ist,i), i=1,2) !Enddo !Close(3) End Subroutine Calc_Trend !--------------------------------------------------------------------- Subroutine Write_GrADS_Stn (Grads_File,ivar,varn) !(X,Lon, Lat, Nst, Nidx, Grads_File, Out_Dir, Field) Integer :: ivar, ist, ic, Nidx, i Character (Len=8) :: STID Character *(*) :: Grads_File, varn !, Out_Dir, Field !Real :: X(Nst, Nidx), Lon(Nst), Lat(Nst) Integer, Parameter :: Nvar=9 Real, Allocatable :: X(:,:) ! (Nst,2,Nvar) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG !---------------------------------------------------------------------------------- Allocate (X(Nst,Nvar)) print*,"In GrADS: Nst=",Nst Do ist=1,Nst X(ist,1)=HWN_Sen (ist,ivar); X(ist,2)=HWDx_Sen(ist,ivar); X(ist,3)=HWDm_Sen(ist,ivar) X(ist,4)=HDays_Sen(ist,ivar); X(ist,5)=HWF_Sen (ist,ivar); X(ist,6)=HWM_Sen (ist,ivar) X(ist,7)=HWA_Sen (ist,ivar); X(ist,8)=HWS_Sen (ist,ivar); X(ist,9)=HWSx_Sen(ist,ivar) Enddo Open (10,file=trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') TIM = 0.0 ! Co the bien doi tu -0.5 .. +0.5; thuong dat bang 0. NFLAG = 1 ! Khong doi ! STID: character(len=8) ==> Station name ! SenSL, Prob: Do i = 1,Nvar ! Nvar chi so, moi chi so duoc xem nhu 1 timeslice NLEV = 1 ! =0 o ban ghi cuoi cung cua 1 timeslice; =1 cho tat ca cac timeslice khac ic = 0 Do ist=1, Nst ic = ic+1 STID = Stn(ist)(1:8) ! Chi lay 8 ky tu dau tien RLAT = Lat(ist); RLON = Lon(ist) WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Var = X(ist,i) WRITE (10) Var ! Gia tri cua bien Enddo ! before new time: NLEV = 0 WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Enddo ! ivar Close(10) !---------------------------------------------------------------------! ! Write control file Station data ! !---------------------------------------------------------------------! open (2001,file=trim(Grads_File)//'.ctl') write (2001,'(2A)') 'dset ','^'//trim(Grads_File)//'.dat' write (2001,'(A)') 'dtype station' write (2001,'(2A)') 'stnmap ','^'//trim(Grads_File)//'.dat'//'.map' write (2001,'(A,F8.1)') 'undef ',Miss write (2001,'(A)') 'title Factor Analysis' write (2001,'(A,I3,1x,A)') 'tdef ',Nvar,' linear jan1980 1mo' ! write (2001,'(A,I3,1x,A)') 'tdef ',1,'linear jan1980 1mo' write (2001,'(A)') 'vars 1' write (2001,'(A)') trim(varn)//' 0 99 Sen slops' write (2001,'(A)') 'endvars' close (2001) !Enddo ! ivar Deallocate (X) End Subroutine Write_GrADS_Stn !--------------------------------------------------------------------- Subroutine Calc_percentile_and_HD (p) Integer :: ireg, ist, iyr, it, N, Nst1, ist1, it1, it2, it0 !Integer :: Mwd = 31, Mwd2 = 15 Real, Allocatable :: Tx1(:,:,:), HD1(:,:,:), qw1(:,:), lon1(:), lat1(:), Reg1(:) Real, Allocatable :: XX(:), X1D(:), q251(:), q751(:) Integer :: N30y ! = NyBL*Ntime = 30*183 !--------------------- HD(:,:,:) = 0; Nst1=0 Do ist=1, Nst ! Calculating qw Do iday=1,Ntime ! Apr-Sep it1=Mwd2+1+iday; it2=it1+Mwd-1; it0=0 Do iyr=1, NyBL ! Cho thoi ky BL Do it=it1,it2 if (Tx0(ist,iyr,it)/=Miss) then it0=it0+1; X(it0)=Tx0(ist,iyr,it) ! X - SL cua NyBL nam endif !print*,it1,it2,it0,X(it0) Enddo Enddo N=it0 ! Dung luong mang X !pause if (.not.Empty_Ar1(X,N,Miss).AND.N>Mwd*NyBL/2) then qw(ist,iday) = percentile(x, N, p, Miss) else qw(ist,iday) = Miss endif Enddo Enddo ! Calculate qw(Nst,Ntime) ! Calculating annual q25, q75 for BL period: Allocate (XX(NyBL),X1D(Ntime)) Do ist=1, Nst Do iyr=1,NyBL X1D(1:Ntime) = Tx01(ist,iyr,1:Ntime) XX(iyr) = Maxval(X1D) Enddo if (.not.Empty_Ar1(XX,NyBL,Miss)) then q25(ist) = percentile(XX, NyBL, 0.25, Miss) q75(ist) = percentile(XX, NyBL, 0.75, Miss) else q25(ist) = Miss; q75(ist) = Miss endif Enddo ! Deallocate (XX,X1D); Allocate (XX(Ntime)) Do ist=1,Nst ! Calculating HD Do iday=1,Ntime if (qw(ist,iday)/=Miss) then Do iyr=1, Nyr1 if (Tx(ist,iyr,iday)>=qw(ist,iday)) HD(ist,iyr,iday)=1 if (Tx(ist,iyr,iday)==Miss) HD(ist,iyr,iday)=Miss Enddo else HD(ist,:,iday)=Miss ! Ngay iday cua tat ca cac nam endif Enddo XX(1:Ntime)=qw(ist,1:Ntime) if (.not.Empty_Ar1(XX,Ntime,Miss)) Nst1=Nst1+1 Enddo ! Xac dinh lai so tram, danh sach tram, Tx va HD: Allocate (Tx1(Nst1,Nyr1,DayofYear),HD1(Nst1,Nyr1,DayofYear),Stname(Nst1), qw1(Nst1,Ntime),lon1(Nst1),lat1(Nst1),Reg1(Nst1)) Allocate (q251(Nst1),q751(Nst1)) ist1=0 Do ist=1, Nst XX(:) = qw(ist,:) if (.not.Empty_Ar1 (XX, Ntime, Miss)) then ist1=ist1+1 qw1(ist1,:) = qw(ist,:) Stname(ist1) = Stn(ist) Reg1(ist1) = Reg(ist) q251(ist1) = q25(ist) q751(ist1) = q75(ist) !print*,Stname(ist1)," h=",h(ist)," q=", q1(ist1) lon1(ist1) = lon(ist); lat1(ist1) = lat(ist) Do iyr=1,Nyr1 Tx1(ist1,iyr,1:Ntime) = Tx(ist,iyr,1:Ntime) HD1(ist1,iyr,1:Ntime) = HD(ist,iyr,1:Ntime) Enddo endif Enddo !print*,"Nst1,ist1=",Nst1,ist1;pause Deallocate (qw,Tx,HD,Stn, lon,lat,Reg,q25,q75); Nst=Nst1 ! New Nst value Allocate (qw(Nst,Ntime),Tx(Nst,Nyr1,Ntime),HD(Nst,Nyr1,Ntime),Stn(Nst),lon(Nst),lat(Nst),Reg(Nst)) Allocate (q25(Nst),q75(Nst)) qw = qw1; Tx = Tx1; HD = HD1; Stn = Stname; Lon = Lon1; Lat = Lat1; Reg = Reg1 ! New q25 = q251; q75 = q751 Deallocate (qw1,Tx1,HD1,Stname,Lon1,Lat1,Reg1, XX, q251, q751) Do ireg=1,Nreg NSt_reg(ireg) = 0 Do ist=1,Nst if (Reg(ist)==ireg) NSt_reg(ireg)=NSt_reg(ireg)+1 Enddo Enddo End Subroutine Calc_percentile_and_HD !---------------------------------------------------------------------- Subroutine Extract_Tx Integer :: ireg, ist, ist1, imon, iday, iyr, iyr1, it, it1 Do ist=1, Nst ! Mang Tx(Nst,Nyr1,NumOfDay_Apr_Sep) iyr1=0; Do iyr=1, Nyrs if (Int(Year(iyr))>=Yr1.AND.Int(Year(iyr))<=Yr2) then ! Chi tinh Nyr1 nam (Yr1-Yr2) iyr1=iyr1+1; it=0 Do imon=Mon1, Mon2 ! Apr - Sep Nday=DayOfMon(imon, Int(Year(iyr))) Do iday=1, Nday it=it+1 ! Tong so ngay tu Apr-Sep Tx(ist,iyr1,it)=dat(ist,iyr,imon,iday) Enddo Enddo endif Enddo ! iyr=1, Nyrs ! Tx0: For calculating qw; Tx01: For q25, q75 iyr1=0 Do iyr=1,Nyrs if (Int(Year(iyr))>=Ybl1.AND.Int(Year(iyr))<=Ybl2) then ! Chi tinh NyBL nam (Ybl1-Ybl2) iyr1=iyr1+1; it=0 Do imon=Mon1-1, Mon2+1 ! Mar - Oct Nday=DayOfMon(imon, Int(Year(iyr))) Do iday=1, Nday it=it+1 ! Tong so ngay tu Mar-Oct Tx0(ist,iyr1,it)=dat(ist,iyr,imon,iday) Enddo Enddo it=0 ! Tx01 Do imon=Mon1, Mon2 ! Apr - Sep Nday=DayOfMon(imon, Int(Year(iyr))) Do iday=1, Nday it=it+1 ! Tong so ngay tu Mar-Oct Tx01(ist,iyr1,it)=dat(ist,iyr,imon,iday) Enddo Enddo endif Enddo ! iyr=1,NyBL 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") ! Avoid all wrong data: Where (Dat>=44) Dat=Miss End Subroutine Read_Tx !----------------------------------------------------------------------- Logical Function Empty_Ar1 (X, N, Miss) ! Tra ve .True. neu X(N) khuyet toan bo SL Integer :: N, i, k Real :: X(N), Miss k = 0 ! Dem so phan tu khuyet cua X Do i=1, N if (X(i) == Miss) k = k + 1 Enddo if (k > N/2) 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" Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90" !----------------------------------------------------------------------- END