! 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) !--------------------------------------------------------------------- ! Trens Analusis: Real, Allocatable, Dimension (:,:) :: HWN_Sen,HWDx_Sen,HWDm_Sen,HDays_Sen,HWF_Sen,HWM_Sen,HWA_Sen ! (Nst,2) !********************************************************************* 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.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 !--------------------------------------------------------------------- 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 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 Ntime=DayofYear ! Tong so ngay cac thang tu Mon1-Mon2 cua cac nam tu Yr1 den Yr2 Allocate (Tx(Nst,Nyr1,DayofYear),HD(Nst,Nyr1,DayofYear)) Call Extract_Tx Allocate (q(Nst),X(Nyr1*Ntime)) !,Hevent(Nyr1,Ntime)) print*,"Nst,Nyr1,Ntime=",Nst,Nyr1,Ntime Call Calc_percentile_and_HD (p) print*,"Nst,Nyr1,Ntime=",Nst,Nyr1,Ntime; pause !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 Call Calc_Trend Call Write_GrADS_Stn ("Sensl",1,"sensl") ! (Grads_File,ivar,varn) Call Write_GrADS_Stn ("prob", 2,"prob") Fname="Heatwave_Chars_Stn.nc" Call Write_Char_to_NC (Fname) CONTAINS !--------------------------------------------------------------------- Subroutine Write_Char_to_NC (Fname) Integer :: iyr, ireg Real, Allocatable :: Years(:), StnID(:),X(:) Character (Len=*) :: Fname Character (Len=50) :: var_ln Allocate (Years(Nyr1), StnID(Nst), X(Nst)) Do iyr=1, Nyr1; Years(iyr)=iyr+Yr1-1; Enddo Do ist=1, Nst; StnID(ist)=ist; Enddo 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 HW days" Call wr_ncdf_2d_apd (Fname, HDays, "nday", 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) 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) ! Call wr_ncdf_1d_apd_char (Fname, Stn,"stname", Nst, "Station's Name", "-", "stn", Miss) ! Trends: ! 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)) ! (Nst,2) 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) Call wr_ncdf_1d_apd_char (Fname, Stn,"stname", Nst, "Station's Name", "-", "stns", 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, Hdy1, Txx, Tx0, Txm, NTx ! N_Spell: So dot trong toan chuoi ! Spells(N_Spell): Cac phan tu la Do dai tung dot Allocate (HWN(Nst,Nyr1),HWDx(Nst,Nyr1),HWDm(Nst,Nyr1),HDays(Nst,Nyr1)) Allocate (HWM(Nst,Nyr1),HWA(Nst,Nyr1),HWF(Nst,Nyr1)) Allocate (Spells(Ntime), X(Ntime), SDate(Ntime), EDate(Ntime)) Do ist=1, Nst if (q(ist) /= Miss) then 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 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 if (Tx(ist,iyr,ispl)>=q(ist)) then ! Hot day Txm=Txm+Tx(ist,iyr,ispl) ! Sum of Temp of hot days Ntx=Ntx+1 ! Num of hot 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 endif 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 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) Hdy1=0 ! So ngay nong Do it=1,Ntime if (X(it)==1) Hdy1=Hdy1+1 Enddo HDays(ist,iyr)=Hdy1 Enddo else HWN(ist,:)=Miss; HWDm(ist,:)=Miss; HWDx(ist,:)=Miss; HDays(ist,:)=Miss; HWM(ist,:)=Miss; HWA(ist,:)=Miss; HDays(ist,:)=Miss endif Enddo 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)) ! (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 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 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 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 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 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 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 Enddo Open (3,file="Trend.txt",status="Unknown") write (3,"(1x,A3,A15,7(4x,A12,4x))") "No","Stn","HWN","HWDx","HWDm","HDays","HWF","HWM","HWA" 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) 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=7 Real, Allocatable :: X(:,:) ! (Nst,2,Nvar) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG !---------------------------------------------------------------------------------- ! 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)) ! (Nst,2) 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) Enddo !OPEN (10,file=trim(Out_Dir)//trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') 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_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, Nst1, ist1 Real, Allocatable :: Tx1(:,:,:), HD1(:,:,:), q1(:), lon1(:), lat1(:) ! Allocate (Tx(Nst,Nyr1,DayofYear),HD(Nst,Nyr1,DayofYear)) !--------------------- HD(:,:,:) = 0; Nst1=0 Do ist=1, Nst !n(ireg) it=0 Do iyr=1, Nyr1 Do iday=1, DayofYear ! =Ntime if (Tx(ist,iyr,iday)/=Miss) then it=it+1; X(it)=Tx(ist,iyr,iday) ! X - SL cua Nyr1 nam endif Enddo Enddo !print*,X(1:it) N=it ! Dung luong mang X if (.not.Empty_Ar1(X,N,Miss).AND.N>Ntime*Nyr1/2) then ! q(ist) = percentile(X,N,p, Miss) ! q=X[F(x)=p] q(ist) = percentile(x, N, p, Miss) else q(ist) = Miss endif print*,"Stn,q=",Stn(ist),q(ist) if (q(ist) /= Miss) then Do iyr=1, Nyr1 Do iday=1, DayofYear if (Tx(ist,iyr,iday)>=q(ist)) HD(ist,iyr,iday)=1 if (Tx(ist,iyr,iday)==Miss) HD(ist,iyr,iday)=Miss Enddo Enddo Nst1=Nst1+1 ! So tram co percntile /= Miss (duoc xu ly) else HD(ist,:,:)=Miss !iyr,iday)=Miss endif Enddo ! Xac dinh lai so tram, danh sach tram, Tx va HD: Allocate (Tx1(Nst1,Nyr1,DayofYear),HD1(Nst1,Nyr1,DayofYear),Stname(Nst1), q1(Nst1),lon1(Nst1),lat1(Nst1)) ist1=0 Do ist=1, Nst if (q(ist) /= Miss) then ist1=ist1+1 q1(ist1) = q(ist) Stname(ist1) = Stn(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:DayofYear) = Tx(ist,iyr,1:DayofYear) HD1(ist1,iyr,1:DayofYear) = HD(ist,iyr,1:DayofYear) Enddo endif Enddo !print*,"Nst1,ist1=",Nst1,ist1;pause Deallocate (q,Tx,HD,Stn, lon,lat); Nst=Nst1 ! New Nst value Allocate (q(Nst),Tx(Nst,Nyr1,DayofYear),HD(Nst,Nyr1,DayofYear),Stn(Nst),lon(Nst),lat(Nst)) q = q1; Tx = Tx1; HD = HD1; Stn = Stname; Lon = Lon1; Lat = Lat1 ! New Deallocate (q1,Tx1,HD1,Stname,Lon1,Lat1) End Subroutine Calc_percentile_and_HD !---------------------------------------------------------------------- Subroutine Extract_Tx Integer :: ireg, ist, ist1, imon, iday, iyr, iyr1 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 !Ntime=it !print*,Stn(ist), Year(iyr) !print "(20F6.1)", Tx(nst1,iyr1,:) endif !pause Enddo 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>=43) 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