Program Create_data Integer, Parameter :: NVung=7, NSta=14, NYr1=1961, NYr2=2007, NDay=31, NYto=8, M=12 Integer, Parameter :: NtrMx = 200 ! So tram cuc dai co SL !Integer, Parameter :: Nmax = (NYr2-NYr1+1)*NDay ! Do dai chuoi Daily data cho mot thang Character *3, Parameter :: Var(NYto) = (/"T2m", "Tx", "Tm", "R", "BH", "Um", "U13", "Vx"/) Character *11 :: StaName (NVung, NSta) Character *3 :: YT(NYto) Character (Len = 30) :: FInp, FOut, FOut1, Fmt Character (Len = 16) :: TenVung(NVung) Character (Len = 11) :: TenTram(NtrMx), VarName(NtrMx), Sample(Nmax) Integer :: NTram Real :: X(Nmax, NtrMx), W1(Nmax) Real*8 :: XX(Nmax, Mmax), X1(Nmax, Mmax), XX_R(Nmax,Mmax) Integer :: YYDD(Nmax), YY(Nmax) Real :: Con !---------------------------------------------------------------------------- Real, Dimension (NVung, NSta, NYto, NYr1:NYr2, M, NDay) :: Dat ! Data array Integer (Kind=1), Dimension (NVung, NSta, NYr1:NYr2, M, NDay) :: & H35, H37, HDry, HDry37, C15, C13, R25, R50 ! H35: Hot day T>35 ! H37: Hot day T>37 ! HDry: Hot and Dry day T>35 and RH<55 ! HDry37: Hot and Dry day T>37 and RH<55 ! C15: Cold day T<15 ! C13: Cold day T<13 ! R25: Heavy Rain R>25 ! R50: Heavy Rain R>50 Integer, Dimension (NVung, NSta, NYr1:NYr2, M+1) :: NH35, NH37, NHDry, NHDry37, NC15, NC13, NR25, NR50, NWork Real :: R_per(NVung*NSta, M), Per ! Luu Phan vi Per cua mua Real :: Rain_Dat(NVung*NSta, M, (NYr2-NYr1+1)*NDay) ! So lieu mua ngay tung thang cua cac tram Real :: W((NYr2-NYr1+1)*NDay) ! Mang lam viec Character *11 DSTram(NVung*NSta) Integer :: N_All_Sta !---------------------------------------------------------------------------- Dat(:,:,:,:,:,:) = -99.0 Open (4,file="SL_ALL.dat", form="Unformatted") print*," Reading..." Do iyt=1,NYto read(4) YT(iyt) Do iv=1,Nvung read(4) TenVung(iv) Do itr=1,NSta read(4) StaName(iv,itr) Do iyr= NYr1,NYr2 read (4) iyear Do iday=1,Nday read(4) idate, (Dat(iv,itr,iyt,iyr,imon,iday),imon=1,M) Enddo Enddo Enddo Enddo Enddo Close(4) ! Tinh percentile cua Rain Per = 0.90 Call Cal_Percentile(4, Per, "Percentile_R_Obs.txt") ! Tinh percentile cua Tx Per = 0.90 Call Cal_Percentile(2, Per, "Percentile_Tx_Obs.txt") ! Tinh percentile cua Tm Per = 0.10 Call Cal_Percentile(3, Per, "Percentile_Tm_Obs.txt") CONTAINS !----------------------------------------------------------------------- Subroutine Cal_Percentile(iyt, Per, FName) Character *(*) :: FName Integer :: iyt Call Extract_dat(iyt) Do ist = 1, N_All_Sta Do imon = 1, M nt = 0 Do it = 1, (NYr2-NYr1+1)*NDay if (Rain_Dat(ist,imon,it) /= -99.) then nt = nt + 1 W(nt) = Rain_Dat(ist,imon,it) endif Enddo if (.NOT.EmptyDat(W, nt)) then R_Per(ist,imon) = percentile_1(W,nt,per) else R_Per(ist,imon) = -99. endif Enddo Enddo ! Write out Percentiles Open (3,file=FName) Write (3,'(A20,12(A5,I2.2))')"Station_name", ("Mon_",j,j=1,M) Do ist = 1, N_All_Sta if (.NOT.EmptyDat(R_Per(ist,:),M)) & write(3,'(A20,12F7.1)') DSTram(ist), (R_Per(ist,j),j=1,M) Enddo Close(3) End Subroutine Cal_Percentile !----------------------------------------------------------------------- Subroutine Extract_dat(iyt) ! Mang Rain_Dat(NVung*NSta, M, (NYr2-NYr1+1)*NDay) chua SL mua ngay tung thang cua cac tram Integer :: iv,itr,iyt,iyr,imon,iday, ist, nt Rain_Dat(:, :, :) = -99. ist = 0 Do iv = 1, NVung Do itr = 1, NSta if (trim(StaName(iv,itr)) /= "") then ist = ist + 1 DSTram(ist) = trim(StaName(iv,itr)) Do imon = 1, M nt = 0 Do iyr = Nyr1, Nyr2 Do iday = 1, NDay nt = nt +1 Rain_Dat(ist,imon,nt) = Dat(iv,itr,iyt,iyr,imon,iday) Enddo Enddo Enddo endif Enddo Enddo N_All_Sta = ist ! Tong so tram co SL End Subroutine Extract_dat !----------------------------------------------------------------------- Function percentile_1(X,N,p) ! Tham khao: http://www.childrens-mercy.org/stats/definitions/percentiles.htm Integer :: N, i, j Real :: X(n), tmp, p, per ! (p value mustbe in [0-1]) Do i = 1, n-1 Do j = i+1, n if (x(j) < x(i)) then tmp = x(j) x(j) = x(i) x(i) = tmp endif Enddo Enddo tmp = p*(n+1) if (tmp < 1) then per = x(1) elseif (tmp > n) then per = x(n) else i = int(tmp) j = i+1 per = x(i) + (x(j)-x(i))*(tmp-real(i)) endif percentile_1 = per End Function Percentile_1 !----------------------------------------------------------------------- Subroutine WriteEventsToFile (Title1, Title2, X, N1, N2, M) Character *(*) Title1, Title2 Integer :: N1, N2, M, X(N1:N2,M), i, j Character*20 :: Fmt write(Fmt,'(A,I3.3,A)') '(I4,',M,'I6)' write(11,'(A20,A40)') Title1, Title2 do i=N1,N2 write(11,Fmt) i, (X(i, j), j=1,M) enddo END Subroutine WriteEventsToFile !----------------------------------------------------------------------- Subroutine NumOfHotDay Real, Parameter :: T_hot35 = 35., T_hot37 = 37., T_cold15 = 15., T_cold13 = 13., R_heavy25 = 25., & R_heavy50 = 50., RH55 = 55. Integer :: iv, ista, iyr, imon, iyto, it Real :: W1(Nday), W2(Nday) H35(:, :, :, :, :) = -99; H37(:, :, :, :, :) = -99; HDry(:, :, :, :, :) = -99; C15(:, :, :, :, :) = -99; C13(:, :, :, :, :) = -99; HDry37(:, :, :, :, :) = -99; R25(:, :, :, :, :) = -99; R50(:, :, :, :, :) = -99 NH35(:, :, :, :) = -99; NH37(:, :, :, :) = -99; NHDry(:, :, :, :) = -99; NC15(:, :, :, :) = -99; NC13(:, :, :, :) = -99; NHDry37(:, :, :, :) = -99; NR25(:, :, :, :) = -99; NR50(:, :, :, :) = -99 Do iv = 1, Nvung Do ista = 1, NSta if (trim(StaName(iv,itr)) /= "") then Do iyr = NYr1, NYr2 Do imon = 1,M NDayOfMon = DayOfMonth(imon,iyr) ! ---------------------------------------------- W1(1:NDayOfMon) = Dat(iv, ista, 1, iyr, imon, 1:NDayOfMon) ! T2m ! Cold day if (.NOT.EmptyDat(W1, NDayOfMon)) then NC15(iv, ista, iyr, imon) = 0; NC13(iv, ista, iyr, imon) = 0 Do it = 1, NDayOfMon if (W1(it) /= -99) then if (W1(it) <= T_cold15) then C15(iv, ista, iyr, imon, it) = 1 NC15(iv, ista, iyr, imon) = NC15(iv, ista, iyr, imon) + 1 endif if (W1(it) <= T_cold13) then C13(iv, ista, iyr, imon, it) = 1 NC13(iv, ista, iyr, imon) = NC13(iv, ista, iyr, imon) + 1 endif else C15(iv, ista, iyr, imon, it) = 0; C13(iv, ista, iyr, imon, it) = 0 endif Enddo endif ! ---------------------------------------------- W1(1:NDayOfMon) = Dat(iv, ista, 2, iyr, imon, 1:NDayOfMon) ! Tx ! Hot day if (.NOT.EmptyDat(W1, NDayOfMon)) then NH35(iv, ista, iyr, imon) = 0; NH37(iv, ista, iyr, imon) = 0 Do it = 1, NDayOfMon if (W1(it) /= -99) then if (W1(it) >= T_hot35) then H35(iv, ista, iyr, imon, it) = 1 NH35(iv, ista, iyr, imon) = NH35(iv, ista, iyr, imon) + 1 endif if (W1(it) >= T_hot37) then H37(iv, ista, iyr, imon, it) = 1 NH37(iv, ista, iyr, imon) = NH37(iv, ista, iyr, imon) + 1 endif else H35(iv, ista, iyr, imon, it) = 0; H37(iv, ista, iyr, imon, it) = 0 endif Enddo endif ! ---------------------------------------------- W1(1:NDayOfMon) = Dat(iv, ista, 2, iyr, imon, 1:NDayOfMon) ! Tx W2(1:NDayOfMon) = Dat(iv, ista, 7, iyr, imon, 1:NDayOfMon) ! RHm ! Hot Dry if ((.NOT.EmptyDat(W1, NDayOfMon)).AND.(.NOT.EmptyDat(W2, NDayOfMon))) then NHDry(iv, ista, iyr, imon) = 0; NHDry37(iv, ista, iyr, imon) = 0; Do it = 1, NDayOfMon if ((W1(it) /= -99).AND.(W2(it) /= -99)) then if ((W1(it) >= T_hot35).AND.(W2(it) <= RH55)) then HDry(iv, ista, iyr, imon, it) = 1 NHDry(iv, ista, iyr, imon) = NHDry(iv, ista, iyr, imon) + 1 endif if ((W1(it) >= T_hot37).AND.(W2(it) <= RH55)) then HDry37(iv, ista, iyr, imon, it) = 1 NHDry37(iv, ista, iyr, imon) = NHDry37(iv, ista, iyr, imon) + 1 endif else HDry(iv, ista, iyr, imon, it) = 0; HDry37(iv, ista, iyr, imon, it) = 0 endif Enddo endif ! ---------------------------------------------- W1(1:NDayOfMon) = Dat(iv, ista, 4, iyr, imon, 1:NDayOfMon) ! R ! Heavy Rainfall if (.NOT.EmptyDat(W1, NDayOfMon)) then NR25(iv, ista, iyr, imon) = 0; NR50(iv, ista, iyr, imon) = 0; Do it = 1, NDayOfMon if (W1(it) /= -99) then if (W1(it) >= R_heavy25) then R25(iv, ista, iyr, imon, it) = 1 NR25(iv, ista, iyr, imon) = NR25(iv, ista, iyr, imon) + 1 endif if (W1(it) >= R_heavy50) then R50(iv, ista, iyr, imon, it) = 1 NR50(iv, ista, iyr, imon) = NR50(iv, ista, iyr, imon) + 1 endif else R25(iv, ista, iyr, imon, it) = 0; R50(iv, ista, iyr, imon, it) = 0 endif Enddo endif ! ---------------------------------------------- Enddo ! Tinh tong cac su kien cho ca nam ! NH35(iv, ista, iyr, M+1) = 0; NH37(iv, ista, iyr, M+1) = 0; NHDry(iv, ista, iyr, M+1) = 0; ! NC15(iv, ista, iyr, M+1) = 0; NC13(iv, ista, iyr, M+1) = 0; NHDry37(iv, ista, iyr, M+1) = 0; ! NR25(iv, ista, iyr, M+1) = 0; NR50(iv, ista, iyr, M+1) = 0 ! Do imon = 1, M ! NH35(iv, ista, iyr, M+1) = NH35(iv, ista, iyr, M+1) + NH35(iv, ista, iyr, imon) ! NH37(iv, ista, iyr, M+1) = NH37(iv, ista, iyr, M+1) + NH37(iv, ista, iyr, imon) ! NHDry(iv, ista,iyr, M+1) = NHDry(iv, ista,iyr, M+1) + NHDry(iv, ista, iyr, imon) ! NHDry37(iv,ista,iyr,M+1) = NHDry37(iv,ista,iyr,M+1) + NHDry37(iv, ista, iyr, imon) ! NC15(iv, ista, iyr, M+1) = NC15(iv, ista, iyr, M+1) + NC15(iv, ista, iyr, imon) ! NC13(iv, ista, iyr, M+1) = NC13(iv, ista, iyr, M+1) + NC13(iv, ista, iyr, imon) ! NR25(iv, ista, iyr, M+1) = NR25(iv, ista, iyr, M+1) + NR25(iv, ista, iyr, imon) ! NR50(iv, ista, iyr, M+1) = NR50(iv, ista, iyr, M+1) + NR50(iv, ista, iyr, imon) ! Enddo Call Annual (NH35(iv, ista, iyr, 1:M), M, NH35(iv, ista, iyr, M+1)) Call Annual (NH37(iv, ista, iyr, 1:M), M, NH37(iv, ista, iyr, M+1)) Call Annual (NC15(iv, ista, iyr, 1:M), M, NC15(iv, ista, iyr, M+1)) Call Annual (NC13(iv, ista, iyr, 1:M), M, NC13(iv, ista, iyr, M+1)) Call Annual (NHDry(iv, ista, iyr, 1:M), M, NHDry(iv, ista, iyr, M+1)) Call Annual (NHDry37(iv, ista, iyr, 1:M), M, NHDry37(iv, ista, iyr, M+1)) Call Annual (NR25(iv, ista, iyr, 1:M), M, NR25(iv, ista, iyr, M+1)) Call Annual (NR50(iv, ista, iyr, 1:M), M, NR50(iv, ista, iyr, M+1)) Enddo endif Enddo Enddo END Subroutine NumOfHotDay !----------------------------------------------------------------------- Subroutine Annual (X, N, Annl) Integer :: N, Annl, X(N), i Annl = 0 Do i=1, N if (X(i) /= -99) then Annl = Annl + X(i) else Annl = -99 exit endif Enddo End Subroutine Annual !----------------------------------------------------------------------- Logical Function EmptyDat(X, N) ! Mang X khong rong neu co it nhat 1 so lieu quan trac Real :: X(N) Integer :: i Logical :: OK OK = .True. Do i=1, N if (X(i) /= -99) then OK = .False. Exit endif Enddo EmptyDat = OK End Function EmptyDat !----------------------------------------------------------------------- Logical Function EnoughDat(X, N) ! Mang X co day du SL neu khong bi khuyet phan tu nao Real :: X(N) Integer :: i Logical :: OK OK = .False. Do i=1,N if ((X(i) /= 0).AND.(X(i) /= -99)) then OK = .True. exit endif Enddo if (.NOT.OK) then EnoughDat = .False. Return endif OK = .True. Do i=1, N if (X(i) == -99) then OK = .False. Exit endif Enddo EnoughDat = OK End Function EnoughDat !----------------------------------------------------------------------- Integer Function LengthSeries(X, N) ! So phan tu cua mang X khac gia tri gan -99 Real :: X(N) Integer :: i, L L = 0 Do i=1, N if ((X(i) /= -99).AND.(X(i) /= 0)) then L = L+1 endif Enddo LengthSeries = L !/30 End Function LengthSeries !----------------------------------------------------------------------- Real Function MaxX(X, N) Real :: X(N), tmp Integer :: i Tmp = X(1) Do i=2, N if (X(i) > Tmp) Tmp = X(i) Enddo MaxX = Tmp End Function MaxX !----------------------------------------------------------------------- Integer Function DayOfMonth ( Mon, Year) Integer :: Mon, Year Select Case (Mon) Case (1,3,5,7,8,10,12) DayOfMonth = 31 Case (4,6,9,11) DayOfMonth = 30 Case (2) if ((MOD(Year, 4) == 0).AND.(MOD(Year, 100) /= 0)) then DayOfMonth = 29 else DayOfMonth = 28 endif if (MOD(Year, 400) == 0) DayOfMonth = 29 End Select End Function DayOfMonth !----------------------------------------------------------------------- END