! All Subroutines !************************* Subroutine Write_2_Thanh Integer :: i Do i = 1, Nidx Open (77,file="Thanh_"//trim(ECE_Var(i))//".txt", Status="unknown") Do ist=1, Nst if (SenSlp(ist,i) /= RMiss) then if (NSample(ist,i) >= NumYear) then ! Lay nhung tram co 20 nam SL tro len if (Prob(ist,i) /= RMiss.AND.Prob(ist,i) <= 0.1) then; isgnt=1; else; isgnt=0; endif Write(77,"(A18,2F18.3,F8.4,2I8)"),Station(ist),Lon(ist),Lat(ist),SenSlp(ist,i),isgnt,NSample(ist,i) endif endif Enddo Close(77) Enddo End Subroutine Write_2_Thanh !---------------------------------------------------------------------- Subroutine Write_CSIRO ! Write to the file Open (66,File="N_ECE_Days_CSIRO.txt",Status="Unknown") ! Number of extreme days write (66,"(A)") "Number of extreme days in each year at each station" write (66,"(A)") "(HD90p: IPCC Hot days, CD10p: IPCC Cold days, RD95p: IPCC Heavy rain days)" write (66,"(A)") "(RH35c: Tx>35 Hot days, RC15c: T2m<15 Cold days, CD15c: Tm<15 Cold days, RR50m: VN Heavy rain days)" Do ist=1, Nst if (trim(Station(ist))/="") then write (66,"(A)")"Station: "//trim(Station(ist)) write (66,"(3(A6,F7.1))") "TX90p=",TX90p(ist), "TN10p=",TN10p(ist), "R95p=",R95p(ist) write (66,"(A4,7A7)") "Year", "HD90p", "RH35c", "CD10p", "CD15c", "RC15c", "RD95p", "RR50m" Do iyr=Yr1, Yr2 if (HD90p(ist,iyr)/=RMiss.OR.CD10p(ist,iyr)/=RMiss.OR.RD95p(ist,iyr)/=RMiss.OR.RH35c(ist,iyr)/=RMiss.OR. & RC15c(ist,iyr)/=RMiss.OR.RR50c(ist,iyr)/=RMiss.OR.CD15c(ist,iyr)/=RMiss) then write (66,"(I4,7F7.1)") iyr, HD90p(ist,iyr),RH35c(ist,iyr), CD10p(ist,iyr),CD15c(ist,iyr), & RC15c(ist,iyr), RD95p(ist,iyr),RR50c(ist,iyr) endif Enddo endif Enddo Close (66) Open (66,File="Trend_ECE_Sen_Slope.txt",Status="Unknown") ! Trend write (66,"(A)") "Trend of Change (Sen's slope values) in the Number of extreme days of each station" write (66,"(A)") "(HD90p: IPCC Hot days, CD10p: IPCC Cold days, RD95p: IPCC Heavy rain days)" write (66,"(A)") "(RH35c: VN Hot days, RC15c: VN Cold days, RR50m: VN Heavy rain days)" write (66,"(A)") " (Unit is days per decade)" write (66,*) write (66,"(A15,6(A10,A9,1x))") "Station", "HD90p","Prob", "RH35c","Prob", "CD10p","Prob", "RC15c","Prob", "RD95p","Prob", "RR50m","Prob" Do ist=1, Nst if (trim(Station(ist))/="") then Wk(1:6) = SenSlp(ist,9:14) if (.not. Empty_Arr1 (Wk,6)) then write (66,"(A15,6(F10.4,F9.3,1x))") Station(ist), (SenSlp(ist,i),Prob(ist,i), i=9,14) endif endif Enddo Close (66) Open (66,File="Trend_ECE_Linear_Slope.txt",Status="Unknown") ! Trend write (66,"(A)") "Trend of Change (Linear slope values) in the Number of extreme days of each station" write (66,"(A)") "(HD90p: IPCC Hot days, CD10p: IPCC Cold days, RD95p: IPCC Heavy rain days)" write (66,"(A)") "(RH35c: VN Hot days, RC15c: VN Cold days, RR50m: VN Heavy rain days)" write (66,"(A)") " (Unit is days per decade; Prob>0: Satisfy 10% significance)" write (66,*) write (66,"(A15,6(A10,A9,1x))") "Station", "HD90p","Prob", "RH35c","Prob", "CD10p","Prob", "RC15c","Prob", "RD95p","Prob", "RR50m","Prob" Do ist=1, Nst if (trim(Station(ist))/="") then Wk(1:6) = A1(ist,9:14) if (.not. Empty_Arr1 (Wk,6)) then write (66,"(A15,6(F10.4,F9.3,1x))") Station(ist), (A1(ist,i),F_F10(ist,i), i=9,14) endif endif Enddo Close (66) End Subroutine Write_CSIRO !---------------------------------------------------------------------- Subroutine Cal_Some_Extreme_Percentile ! Calculate ECEs ! 1:R 2:Tx 3:Tm 4:T2m 5:Um 6:U13 7:Vx 8:BH !------------------------------------------------ For CSIRO --------------------- TX90p(:)=RMiss; TN10p(:)=RMiss; R95p(:)=RMiss; HD90p(:,:)=RMiss; CD10p(:,:)=RMiss; RD95p(:,:)=RMiss RH35c(:,:)=RMiss; RC15c(:,:)=RMiss; RR50c(:,:)=RMiss; RH32(:,:)=RMiss; RR01_Ann(:,:)=RMiss RH35(:,:)=RMiss; RH37(:,:)=RMiss; RC15(:,:)=RMiss; RC13(:,:)=RMiss; RR50(:,:)=RMiss; RR100(:,:)=RMiss Do ist=1, Nst Call GetDataBaseline(ist, 2, Ndys) ! Tx TX90p(ist) = percentile(Wx, Ndys, P90) Call GetDataBaseline(ist, 3, Ndys) ! Tm TN10p(ist) = percentile(Wx, Ndys, P10) Call GetDataBaseline(ist, 1, Ndys) ! R R95p (ist) = percentile(Wx, Ndys, P95) Do iyr=Yr1, Yr2 HD90p(ist,iyr) = Extreme_Days(2, ist, iyr, TX90p(ist), .true.) CD10p(ist,iyr) = Extreme_Days(3, ist, iyr, TN10p(ist), .false.) RD95p(ist,iyr) = Extreme_Days(1, ist, iyr, R95p (ist), .true.) CD15c(ist,iyr) = Extreme_Days(3, ist, iyr, 15., .false.) RH35c(ist,iyr) = Extreme_Days(2, ist, iyr, 35., .true.) RC15c(ist,iyr) = Extreme_Days(4, ist, iyr, 15., .false.) RR50c(ist,iyr) = Extreme_Days(1, ist, iyr, 50., .true.) RH35(ist,iyr) = RH35c(ist,iyr) RC15(ist,iyr) = RC15c(ist,iyr) RH32(ist,iyr) = Extreme_Days(2, ist, iyr, 32., .true.) RH37(ist,iyr) = Extreme_Days(2, ist, iyr, 37., .true.) if (RH37(ist,iyr)>50) print"(A20,A,F7.1)", trim(Station(ist))," ", RH37(ist,iyr) RC13(ist,iyr) = Extreme_Days(4, ist, iyr, 13., .false.) RR01_Ann(ist,iyr) = Extreme_Days(1, ist, iyr, 0.1, .true.) RR50(ist,iyr) = RR50c(ist,iyr) RR100(ist,iyr) = Extreme_Days(1, ist, iyr, 100., .true.) Enddo Enddo End Subroutine Cal_Some_Extreme_Percentile !---------------------------------------------------------------------- Subroutine Cal_Extrem_Max_Min_Values (ivar, Xmon, Xmon_cli, Xyear, Xcli, Maxx) Integer :: ivar, ist, iyr, imon, iday, N, N1, N2 Real :: Xmon(Nst, Yr1:Yr2, Nmon), Xmon_cli(Nst, Nmon), Xyear(Nst, Yr1:Yr2), Xcli(Nst) Real :: X2((Yr2-Yr1+1)), X1(Nmon*Nday), X(Nday) Logical :: Maxx Xmon = RMiss; Xmon_cli = RMiss; Xyear = RMiss; Xcli = RMiss Do ist=1, Nst if (.not.Empty_Sta (ivar,ist)) then N2=0 Do iyr=Yr1, Yr2 N1=0 Do imon=1, Nmon N=0 Do iday=1, Nday if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ! if (DatAll(ivar,ist,iyr,imon,iday)>43.AND.ivar>1) then ! print"(A21,A,A6,i4.4,2(A,I2.2))",trim(Station(ist))," var=",Var(ivar),iyr,"/",imon,"/",iday ! endif N=N+1; X(N) = DatAll(ivar,ist,iyr,imon,iday) N1=N1+1; X1(N1) = DatAll(ivar,ist,iyr,imon,iday) endif Enddo if (N>=DayMon) then ! tren 25 ngay trong thang co SL if (Maxx) then; Xmon(ist,iyr,imon)=Xmax(X,N); else; Xmon(ist,iyr,imon)=Xmin(X,N); endif endif Enddo if (N1>=DayYear) then ! Tren 330 ngay trong nam co SL if (Maxx) then; Xyear(ist,iyr)=Xmax(X1,N1); else; Xyear(ist,iyr)=Xmin(X1,N1); endif endif if (Xyear(ist,iyr) /= RMiss) then N2=N2+1; X2(N2) = Xyear(ist,iyr) endif Enddo if (N2>=NumYear) then ! Tren 20 nam co SL if (Maxx) then; Xcli(ist)=Xmax(X2,N2); else; Xcli(ist)=Xmin(X2,N2); endif endif ! Tinh XMon_Cli Do imon=1, Nmon N2=0 Do iyr=Yr1, Yr2 if (Xmon(ist,iyr,imon) /= RMiss) then N2=N2+1; X2(N2)=Xmon(ist,iyr,imon) endif if (N2>=NumYear) then ! Tren 20 nam co SL if (Maxx) then; XMon_Cli(ist,imon)=Xmax(X2,N2); else; XMon_Cli(ist,imon)=Xmin(X2,N2); endif endif Enddo Enddo endif Enddo End Subroutine Cal_Extrem_Max_Min_Values !---------------------------------------------------------------------- Subroutine GetDataBaseline(ist, ivar, Ndys) ! Get data of ivar, ist and return into the Wx(Ndys) Integer :: ist, ivar, iyr, imon, iday, NdayOfMon, ic, Ndys ic = 0 Do iyr = Yr1Base, Yr2Base Do imon = 1,Nmon NdayOfMon = DayOfMonth(imon,iyr) Do iday = 1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ic = ic+1; Wx(ic) = DatAll(ivar,ist,iyr,imon,iday) endif Enddo Enddo Enddo Ndys = ic End Subroutine GetDataBaseline !---------------------------------------------------------------------- Subroutine GetDataOneYear(ist, ivar, iyr, Ndys) ! Get data of ivar, ist, iyr and return into the Wx(Ndys) if Ndys > 330 days Integer :: ist, ivar, iyr, imon, iday, NdayOfMon, ic, Ndys ic = 0 Do imon = 1,Nmon NdayOfMon = DayOfMonth(imon,iyr) Do iday = 1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ic = ic+1; Wx(ic) = DatAll(ivar,ist,iyr,imon,iday) endif Enddo Enddo if (ic>=DayYear) then ; Ndys = ic; else ; Ndys = 0; endif End Subroutine GetDataOneYear !---------------------------------------------------------------------- Subroutine ECE_Annual_Mean (X,Y) Integer :: ist, iyr, ic Real :: X(Nst, Yr1:Yr2), Y(Nst), Tmp Do ist=1, Nst Tmp = 0; ic = 0; Y(ist) = RMiss Do iyr=Yr1, Yr2 if (X(ist,iyr) /= RMiss) then ic=ic+1; Tmp = Tmp+X(ist,iyr) endif Enddo ! if (real(ic)/real(Yr2-Yr1+1)>=0.90) Y(ist) = Tmp/real(ic) ! 90% so nam co SL if (ic >= NumYear) then; Y(ist) = Tmp/real(ic); else; Y(ist) = RMiss; endif ! It nhat 20 nam co co SL ! if (Y(ist)/=RMiss.and.Y(ist)>50) print"(A20,A,F7.1)",trim(station(ist))," ", Y(ist) Enddo End Subroutine ECE_Annual_Mean !------------------------------------------------------------------------------------- Logical Function Empty_Sta (ivar,ist) Integer :: ist, iyr, imon, iday, ivar Logical :: Empty Empty = .true. Do iyr=Yr1, Yr2 Do imon=1, Nmon Do iday=1, Nday if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then Empty = .false.; Exit endif Enddo Enddo Enddo Empty_Sta = Empty End Function Empty_Sta !------------------------------------------------------------------------------------- Logical Function EmptyDat(X, N) ! Mang X khong rong neu co it nhat 1 so lieu quan trac Integer :: N, i Real :: X(N) Logical :: OK OK = .True. Do i=1, N if (X(i) /= RMiss) then OK = .False. Exit endif Enddo EmptyDat = OK End Function EmptyDat !------------------------------------------------------------------------------------- Logical Function Empty_Arr1 (W1,N) Integer :: N,i Real :: W1(N) Logical :: Empty Empty = .true. Do i=1,N if (W1(i) /= RMiss) then Empty = .false.; Exit endif Enddo Empty_Arr1 = Empty End Function Empty_Arr1 !------------------------------------------------------------------------------------- Subroutine Reinit_all NH35(:, :, :) = IMiss; NH37(:, :, :) = IMiss; NC15(:, :, :) = IMiss; NC13(:, :, :) = IMiss; NR50(:, :, :) = IMiss; NR100(:, :, :) = IMiss T2m(:,:) = RMiss; R(:,:) = RMiss; RH35(:,:) = RMiss; RH37(:,:) = RMiss; RC15(:,:) = RMiss; RC13(:,:) = RMiss Tmon(:,:,:) = RMiss; Rmon(:,:,:) = RMiss ! 2013-05-06 RR50(:,:) = RMiss; RR100(:,:) = RMiss; SenSlp(:,:) = RMiss; Prob(:,:) = RMiss ECE_Ave(:,:) = RMiss; STT_Tram = RMiss ! H32(:, :, :, :) = IMiss; R01(:, :, :, :) = IMiss ! Cafe NH32(:,:,:) = IMiss; NR01(:,:,:) = IMiss ! Cafe RR01_Ann(:,:) = RMiss; RR01_Dry(:,:) = RMiss; RR01_Wet(:,:) = RMiss; ! Cafe Eva(:,:) = RMiss; Eva_Dry(:,:) = RMiss; Eva_Wet(:,:) = RMiss ! Cafe Eva_mon(:,:,:) = RMiss T2m_Cli(:) = RMiss; Tx_Cli(:) = RMiss; Tm_Cli(:) = RMiss; Tx_mon_Cli(:,:) = RMiss; Tm_mon_Cli(:,:) = RMiss ! H35(:, :, :, :) = IMiss; H37(:, :, :, :) = IMiss; C15(:, :, :, :) = IMiss; C13(:, :, :, :) = IMiss; ! R50(:, :, :, :) = IMiss; R100(:, :, :, :) = IMiss End Subroutine Reinit_all !------------------------------------------------------------------------------------- Subroutine Reinit_ECE Integer :: ist, iyr Real :: W(Yr1:Yr2) Do ist = 1, Nst W(:) = RH35(ist,:); if (Sta_No_ECE(W)) RH35(ist,:) = RMiss W(:) = RH37(ist,:); if (Sta_No_ECE(W)) RH37(ist,:) = RMiss W(:) = RC15(ist,:); if (Sta_No_ECE(W)) RC15(ist,:) = RMiss W(:) = RC13(ist,:); if (Sta_No_ECE(W)) RC13(ist,:) = RMiss W(:) = RR50(ist,:); if (Sta_No_ECE(W)) RR50(ist,:) = RMiss W(:) = RR100(ist,:); if (Sta_No_ECE(W)) RR100(ist,:) = RMiss !Cafe W(:) = RH32(ist,:); if (Sta_No_ECE(W)) RH32(ist,:) = RMiss W(:) = RR01_Ann(ist,:); if (Sta_No_ECE(W)) RR01_Ann(ist,:) = RMiss Enddo End Subroutine Reinit_ECE !---------------------------------------------------------------------------------- Logical Function Sta_No_ECE (X) Integer :: ist, iyr Real :: X(Yr1:Yr2) Logical :: No No = .true. Do iyr=Yr1, Yr2 if (X(iyr) /= 0. .And. X(iyr) /= RMiss) then No = .false.; Exit Endif Enddo Sta_No_ECE = No End Function Sta_No_ECE !---------------------------------------------------------------------------------- Subroutine Cal_Linear_Trend Integer :: iyr, ist Do iyr=Yr1, Yr2 Wt(iyr-Yr1+1)=iyr ! Time array Enddo ! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ! T2m RH32 RH35 RH37 R Eva RR01_Ann RR50 R_Wet R_Dry RR01_Wet RR01_Dry Eva_Wet Eva_Dry ! 15 16 17 18 19 20 21 22 23 24 25 26 ! RC15 RC13 HD90p CD10p RD95p RR100 Tx TXx Tm TNn Rx1day T2X Do ist=1, Nst !print*,Station(ist) if (.not.Empty_Sta (4,ist)) then ! T2m Wk(1:NYear) = T2m(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,1)=HA1; F_F10(ist,1)=F Wk(1:NYear) = RC15(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,15)=HA1; F_F10(ist,15)=F Wk(1:NYear) = RC13(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,16)=HA1; F_F10(ist,16)=F Wk(1:NYear) = T2X(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,26)=HA1; F_F10(ist,26)=F endif if (.not.Empty_Sta (1,ist)) then ! R Wk(1:NYear) = R(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,5)=HA1; F_F10(ist,5)=F A1(ist,5)=Decade_Trend(A1(ist,5),R_Cli(ist)) Wk(1:NYear) = RR50(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,8)=HA1; F_F10(ist,8)=F Wk(1:NYear) = RR100(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,20)=HA1; F_F10(ist,20)=F Wk(1:NYear) = RD95p(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,19)=HA1; F_F10(ist,19)=F Wk(1:NYear) = RR01_Ann(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,7)=HA1; F_F10(ist,7)=F Wk(1:NYear) = R_Wet(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,9)=HA1; F_F10(ist,9)=F A1(ist,9)=Decade_Trend(A1(ist,9),R_Wet_Cli(ist)) Wk(1:NYear) = R_Dry(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,10)=HA1; F_F10(ist,10)=F A1(ist,10)=Decade_Trend(A1(ist,10),R_Dry_Cli(ist)) Wk(1:NYear) = RR01_Wet(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,11)=HA1; F_F10(ist,11)=F Wk(1:NYear) = RR01_Dry(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,12)=HA1; F_F10(ist,12)=F Wk(1:NYear) = Rx1day(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,25)=HA1; F_F10(ist,25)=F endif if (.not.Empty_Sta (2,ist)) then ! Tx Wk(1:NYear) = RH32(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,2)=HA1; F_F10(ist,2)=F Wk(1:NYear) = RH35(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,3)=HA1; F_F10(ist,3)=F Wk(1:NYear) = RH37(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,4)=HA1; F_F10(ist,4)=F Wk(1:NYear) = HD90p(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,17)=HA1; F_F10(ist,17)=F Wk(1:NYear) = Tx(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,21)=HA1; F_F10(ist,21)=F Wk(1:NYear) = TXx(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,22)=HA1; F_F10(ist,22)=F endif if (.not.Empty_Sta (3,ist)) then ! Tm Wk(1:NYear) = CD10p(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,18)=HA1; F_F10(ist,18)=F Wk(1:NYear) = Tm(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,23)=HA1; F_F10(ist,23)=F Wk(1:NYear) = TNn(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,24)=HA1; F_F10(ist,24)=F endif if (.not.Empty_Sta (8,ist)) then ! Evapotranspiration Wk(1:NYear) = Eva(ist,Yr1:Yr2); Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,6)=HA1; F_F10(ist,6)=F A1(ist,6)=Decade_Trend(A1(ist,6),Eva_Cli(ist)) Wk(1:NYear) = Eva_Dry(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,14)=HA1; F_F10(ist,14)=F A1(ist,14)=Decade_Trend(A1(ist,14),Eva_Dry_Cli(ist)) Wk(1:NYear) = Eva_Wet(ist,Yr1:Yr2);Call Linear_Trend(Wt,Wk,NYear,HA1,F,RMiss); A1(ist,13)=HA1; F_F10(ist,13)=F A1(ist,13)=Decade_Trend(A1(ist,13),Eva_Wet_Cli(ist)) endif Enddo ! Convert Trendsto the Unit per decade Do i=1, Nidx Do ist=1, Nst if (A1(ist,i) /= RMiss) A1(ist,i)=A1(ist,i)*10. Enddo Enddo End Subroutine Cal_Linear_Trend !---------------------------------------------------------------------- Real Function Decade_Trend(trend, Cli) Real :: trend, Cli, tmp if (trend /= RMiss .AND. Cli /= RMiss) then Decade_Trend = trend/Cli*100. else Decade_Trend = RMiss endif End Function Decade_Trend !---------------------------------------------------------------------- Subroutine Cal_Sen_Annual Integer :: ist ! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ! T2m RH32 RH35 RH37 R Eva RR01_Ann RR50 R_Wet R_Dry RR01_Wet RR01_Dry Eva_Wet Eva_Dry ! 15 16 17 18 19 20 21 22 23 24 25 26 ! RC15 RC13 HD90p CD10p RD95p RR100 Tx TXx Tm TNn Rx1day T2X Do ist=1, Nst if (.not.Empty_Sta (4,ist)) then ! T2m Wk(1:NYear) = T2m(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,1)=SenSlope; Prob(ist,1)=probly; NSample(ist,1)=N1 Wk(1:NYear) = RC15(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,15)=SenSlope; Prob(ist,15)=probly NSample(ist,15)=N1 Wk(1:NYear) = RC13(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,16)=SenSlope; Prob(ist,16)=probly NSample(ist,16)=N1 Wk(1:NYear) = T2X(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,26)=SenSlope; Prob(ist,26)=probly NSample(ist,26)=N1 endif if (.not.Empty_Sta (1,ist)) then ! R Wk(1:NYear) = R(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,5)=SenSlope; Prob(ist,5)=probly SenSlp(ist,5)=Decade_Trend(SenSlp(ist,5),R_Cli(ist)); NSample(ist,5)=N1 !>>> New Wk(1:NYear) = RR50(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,8)=SenSlope; Prob(ist,8)=probly NSample(ist,8)=N1 Wk(1:NYear) = RR100(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,20)=SenSlope; Prob(ist,20)=probly NSample(ist,20)=N1 Wk(1:NYear) = RR01_Ann(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,7)=SenSlope; Prob(ist,7)=probly NSample(ist,7)=N1 Wk(1:NYear) = R_Dry(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,10)=SenSlope; Prob(ist,10)=probly SenSlp(ist,10)=Decade_Trend(SenSlp(ist,10),R_Dry_Cli(ist)) !>>> New NSample(ist,10)=N1 Wk(1:NYear) = R_Wet(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,9)=SenSlope; Prob(ist,9)=probly SenSlp(ist,9) =Decade_Trend(SenSlp(ist, 9),R_Wet_Cli(ist)) !>>> New NSample(ist,9)=N1 Wk(1:NYear) = RR01_Dry(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,12)=SenSlope; Prob(ist,12)=probly NSample(ist,12)=N1 Wk(1:NYear) = RR01_Wet(ist,Yr1:Yr2);Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,11)=SenSlope; Prob(ist,11)=probly NSample(ist,11)=N1 Wk(1:NYear) = RD95p(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,19)=SenSlope; Prob(ist,19)=probly NSample(ist,19)=N1 Wk(1:NYear) = Rx1day(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,25)=SenSlope; Prob(ist,25)=probly SenSlp(ist,25)=Decade_Trend(SenSlp(ist,25),Rx1day_Cli(ist)) NSample(ist,25)=N1 endif if (.not.Empty_Sta (2,ist)) then ! Tx Wk(1:NYear) = RH35(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,3)=SenSlope; Prob(ist,3)=probly NSample(ist,3)=N1 Wk(1:NYear) = RH37(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,4)=SenSlope; Prob(ist,4)=probly NSample(ist,4)=N1 Wk(1:NYear) = RH32(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,2)=SenSlope; Prob(ist,2)=probly NSample(ist,2)=N1 Wk(1:NYear) = HD90p(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,17)=SenSlope; Prob(ist,17)=probly NSample(ist,17)=N1 Wk(1:NYear) = Tx(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,21)=SenSlope; Prob(ist,21)=probly NSample(ist,21)=N1 Wk(1:NYear) = TXx(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,22)=SenSlope; Prob(ist,22)=probly NSample(ist,22)=N1 endif if (.not.Empty_Sta (3,ist)) then ! Tm Wk(1:NYear) = CD10p(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,18)=SenSlope; Prob(ist,18)=probly NSample(ist,18)=N1 Wk(1:NYear) = Tm(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,23)=SenSlope; Prob(ist,23)=probly NSample(ist,23)=N1 Wk(1:NYear) = TNn(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,24)=SenSlope; Prob(ist,24)=probly NSample(ist,24)=N1 endif if (.not.Empty_Sta (8,ist)) then ! Evapotranspiration Wk(1:NYear) = Eva(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,6)=SenSlope; Prob(ist,6)=probly SenSlp(ist,6)=Decade_Trend(SenSlp(ist,6),Eva_Cli(ist)) !>>> New NSample(ist,6)=N1 Wk(1:NYear) = Eva_Dry(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,14)=SenSlope; Prob(ist,14)=probly SenSlp(ist,14)=Decade_Trend(SenSlp(ist,14),Eva_Dry_Cli(ist)) !>>> New NSample(ist,14)=N1 Wk(1:NYear) = Eva_Wet(ist,Yr1:Yr2); Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); SenSlp(ist,13)=SenSlope; Prob(ist,13)=probly SenSlp(ist,13)=Decade_Trend(SenSlp(ist,13),Eva_Wet_Cli(ist)) !>>> New NSample(ist,13)=N1 endif Enddo ! Convert Trendsto the Unit per decade Where (Prob > 0.1) Prob = RMiss Do i=1, Nidx Do ist=1, Nst if (SenSlp(ist,i) /= RMiss) SenSlp(ist,i)=SenSlp(ist,i)*10. Enddo Enddo End Subroutine Cal_Sen_Annual !---------------------------------------------------------------------- Subroutine Cal_Sen_Mon Integer :: ist, imon SenTmon(:,:)=RMiss; SenRmon(:,:)=RMiss; SenEva_mon(:,:)=RMiss; SenRR01_mon(:,:)=RMiss Do ist=1, Nst Do imon=1, Nmon Wk(1:NYear) = Tmon(ist,Yr1:Yr2,imon) if (.not.Empty_Arr1(Wk,NYear)) then Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); if (SenSlope/=RMiss) SenTmon(ist,imon)=SenSlope*10. else SenTmon(ist,imon)=RMiss endif Wk(1:NYear) = Rmon(ist,Yr1:Yr2,imon) if (.not.Empty_Arr1(Wk,NYear)) then Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); if (SenSlope/=RMiss) SenRmon(ist,imon)=Decade_Trend(SenSlope,Rmon_Cli(ist,imon))*10. !!>> New else SenRmon(ist,imon)=RMiss endif Wk(1:NYear) = Eva_mon(ist,Yr1:Yr2,imon) if (.not.Empty_Arr1(Wk,NYear)) then Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); if (SenSlope/=RMiss) SenEva_mon(ist,imon)=Decade_Trend(SenSlope,Eva_mon_Cli(ist,imon))*10 !!>> New else SenEva_mon(ist,imon)=RMiss endif Wk(1:NYear) = NR01(ist,Yr1:Yr2,imon) if (.not.Empty_Arr1(Wk,NYear)) then Call MK_Sen(Wk,NYear,RMiss,SenSlope,probly,N1); if (SenSlope/=RMiss) SenRR01_mon(ist,imon)=SenSlope*10 !!>> (days/decade, New else SenRR01_mon(ist,imon)=RMiss endif Enddo print*,ist," ",Station(ist),SenTmon(ist,:) Enddo End Subroutine Cal_Sen_Mon !------------------------------------------------------------------------------------- Subroutine Annual_Means (ivar, X, X_Cli, Nst, Yr1, Yr2, Cumulative) Integer :: ist, iyr, imon, ivar, iday, NdayOfMon, NdayOfYear, ic, Nst, Yr1, Yr2 Real :: Tmp Real :: X(NSt,Yr1:Yr2), X_Cli(Nst) Logical :: Cumulative Do ist=1, Nst Do iyr=Yr1, Yr2 ic=0; Tmp=0; NDayOfYear=0; X(ist, iyr) = RMiss Do imon=1, Nmon NdayOfMon = DayOfMonth(imon,iyr); NDayOfYear = NDayOfYear+NdayOfMon Do iday=1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ic = ic+1; Tmp = Tmp+DatAll(ivar,ist,iyr,imon,iday) endif Enddo Enddo if (ic >= DayYear) then if (Cumulative) then X(ist,iyr) = Tmp/real(ic)*NDayOfYear else X(ist,iyr) = Tmp/real(ic) endif endif Enddo ic = 0; Tmp = 0. Do iyr=Yr1, Yr2 if (X(ist,iyr) /= RMiss) then ic = ic+1; Tmp = Tmp + X(ist,iyr) endif Enddo if (ic>=NumYear) then; X_Cli(ist) = Tmp/real(ic); else; X_Cli(ist) = RMiss; endif ! TB nhieu nam cua >=20 nam Enddo End Subroutine Annual_Means !------------------------------------------------------------------------------------- Subroutine Monthly_Means (ivar, Xmon, Xmon_Cli, Nst, Yr1, Yr2, Cumulative) Integer :: ist, iyr, imon, ivar, iday, NdayOfMon, NdayOfYear, ic, Nst, Yr1, Yr2, NYear, NYr80 Real :: Tmp Real :: Xmon(NSt,Yr1:Yr2,Nmon), Xmon_Cli(Nst, Nmon) Logical :: Cumulative Xmon(:,:,:) = RMiss; Xmon_Cli(:,:) = RMiss; NYr80 = NumYear !Nint (0.8*NYear) Do ist=1, Nst ! R if (.not.Empty_Sta (ivar,ist)) then Do iyr=Yr1, Yr2 Do imon=1, Nmon ic = 0; Tmp = 0.; NdayOfMon = DayOfMonth(imon,iyr) Do iday=1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ic = ic+1; Tmp = Tmp+DatAll(ivar,ist,iyr,imon,iday) endif Enddo if (ic>=DayMon) then ! Trung binh thang phai co it nhat 25 ngay SL if (Cumulative) then Xmon(ist,iyr,imon) = Tmp/real(ic)*NdayOfMon else Xmon(ist,iyr,imon) = Tmp/real(ic) endif else Xmon(ist,iyr,imon) = RMiss endif Enddo Enddo endif Enddo ! Tinh Trung binh nhieu nam Do ist=1, Nst Do imon=1, Nmon ic = 0; Tmp = 0. Do iyr=Yr1, Yr2 if (Xmon(ist,iyr,imon) /= RMiss) then; ic = ic+1; Tmp = Tmp + Xmon(ist,iyr,imon); endif Enddo ! if (ic>0) then; Xmon_Cli(ist,imon) = Tmp/real(ic); else; Xmon_Cli(ist,imon) = RMiss; endif ! Dam bao co it nhat 20 so nam co SL if (ic >= NumYear) then; Xmon_Cli(ist,imon) = Tmp/real(ic); else; Xmon_Cli(ist,imon) = RMiss; endif Enddo Enddo End Subroutine Monthly_Means !------------------------------------------------------------------------------------- Subroutine Daily_Means (ivar, Xdy_Cli, Nst, Yr1, Yr2) Integer :: ist, iyr, NYear, NYr80, imon, ivar, iday, idy, NdayOfMon, ic, Nst, Yr1, Yr2 Real :: Tmp Real :: Xdy_Cli(Nst, NDayOfYear) Logical :: Cumulative Xdy_Cli(:,:) = RMiss; NYr80 = NumYear ! Nint (0.8*NYear) Do ist=1, Nst ! R if (.not.Empty_Sta (ivar,ist)) then idy = 0 ! Dem so ngay trong nam Do imon=1, Nmon NdayOfMon = DayOfMonth(imon,2011) ! Khong phai nam nhuan Do iday=1, NdayOfMon idy = idy+1; ic = 0; Tmp = 0. Do iyr=Yr1, Yr2 if (DatAll(ivar,ist,iyr,imon,iday) /= RMiss) then ic = ic+1; Tmp = Tmp+DatAll(ivar,ist,iyr,imon,iday) endif Enddo if (ic>=NumYear) then ! Phai co it nhat 20 nam SL Xdy_Cli(ist,idy) = Tmp/real(ic) else Xdy_Cli(ist,idy) = RMiss endif Enddo Enddo endif Enddo End Subroutine Daily_Means !------------------------------------------------------------------------------------- Real Function Extreme_Days(ivar, ist, iyr, Thres, Above) Integer :: ist, ivar, iyr, imon, iday, NdayOfMon, ic, NMiss Real :: Thres Logical :: Above ! =.true. if value exceeds Thres otherwise =.false. ic = 0; NMiss = 0 Do imon=1, Nmon NdayOfMon = DayOfMonth(imon,iyr) Do iday=1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) == RMiss) then NMiss = NMiss + 1 else if (Above) then if (DatAll(ivar,ist,iyr,imon,iday) >= Thres) ic=ic+1 else if (DatAll(ivar,ist,iyr,imon,iday) <= Thres) ic=ic+1 endif endif Enddo Enddo if (NMiss <= NDayOfYear-NumYear) then ; Extreme_Days = ic; else; Extreme_Days = RMiss; endif ! Nhieu nhat 35 ngay trong nam khuyet SL End Function Extreme_Days !------------------------------------------------------------------------- Real Function Extreme_Days_In_Month(ivar, ist, iyr, imon, Thres, Above) Integer :: ist, ivar, iyr, imon, iday, NdayOfMon, ic, NMiss Real :: Thres Logical :: Above ! =.true. if value exceeds Thres otherwise =.false. ic = 0; NMiss = 0 NdayOfMon = DayOfMonth(imon,iyr) Do iday=1, NdayOfMon if (DatAll(ivar,ist,iyr,imon,iday) == RMiss) then NMiss = NMiss + 1 else if (Above) then if (DatAll(ivar,ist,iyr,imon,iday) >= Thres) ic=ic+1 else if (DatAll(ivar,ist,iyr,imon,iday) <= Thres) ic=ic+1 endif endif Enddo if (NMiss <= NdayOfMon-DayMon) then ; Extreme_Days_In_Month = ic; else; Extreme_Days_In_Month = RMiss; endif ! Nhieu nhat 5 ngay/thang khuyet SL End Function Extreme_Days_In_Month !---------------------------------------------------------------------------------- Subroutine Num_Of_Extreme_Days Real, Parameter :: T_hot32=32.,T_hot35=35., T_hot37=37., T_cold15=15., T_cold13=13. Real, Parameter :: R_heavy50=50., R_heavy100=100., R_01=0.1 Integer :: ist, iyr, imon, ivar, iday, NdayOfMon, ic, ic1 Real :: W1(Nday) ! (/"R", "Tx", "Tm", "T2m", "Um", "U13", "Vx", "BH"/) ! 1 2 3 4 5 6 7 8 Do ist=1, Nst if (.not.Empty_Sta (4,ist)) then Do iyr=Yr1, Yr2 Do imon=1, Nmon ! So ngay ret dam, ret hai trong thang NC13(ist, iyr, imon) = Extreme_Days_In_Month(4, ist, iyr, imon, T_cold13, .false.) NC15(ist, iyr, imon) = Extreme_Days_In_Month(4, ist, iyr, imon, T_cold15, .false.) Enddo Enddo endif Enddo ! ---------------------------------------------- Do ist=1, Nst if (.not.Empty_Sta (2,ist)) then Do iyr=Yr1, Yr2 Do imon=1, Nmon ! So ngay nang nong trong thang NH32(ist, iyr, imon) = Extreme_Days_In_Month(2, ist, iyr, imon, T_hot32, .true.) NH35(ist, iyr, imon) = Extreme_Days_In_Month(2, ist, iyr, imon, T_hot35, .true.) NH37(ist, iyr, imon) = Extreme_Days_In_Month(2, ist, iyr, imon, T_hot37, .true.) Enddo Enddo endif Enddo ! ---------------------------------------------- Do ist=1, Nst if (.not.Empty_Sta (1,ist)) then Do iyr=Yr1, Yr2 Do imon=1, Nmon ! So ngay mua, mua lon va rat lon trong thang NR01(ist, iyr, imon) = Extreme_Days_In_Month(1, ist, iyr, imon, R_01, .true.) NR50(ist, iyr, imon) = Extreme_Days_In_Month(1, ist, iyr, imon, R_heavy50, .true.) NR100(ist, iyr, imon) = Extreme_Days_In_Month(1, ist, iyr, imon, R_heavy100, .true.) Enddo Enddo endif Enddo End Subroutine Num_Of_Extreme_Days !------------------------------------------------------------------------------------- Real Function Xmax (X, N) Integer :: N, i Real :: X(N), tmp tmp = -9999. Do i=1, N if (X(i)>tmp) tmp = X(i) Enddo Xmax = tmp End Function Xmax !------------------------------------------------------------------------------------- Real Function Xmin (X, N) Integer :: N, i Real :: X(N), tmp tmp = 99999. Do i=1, N if (X(i)>>> South to North ! if (IDX(ist) == 0) Cycle if (IDX(ist) == 1) then if (Leng_Series(X(ist,1:Nmon),Nmon) >= 11) then ! Phai co du 11 thang SL write (99,"(A)",Advance="No") trim(Station(ist))//"|" ! Ylabs ic = ic+1; Y(ic,1:Nmon)=X(ist,1:Nmon) ! print*,ic," ",trim(Station(ist)) endif endif Enddo write (99,"(A,I7)") '"', ic ! Ylabs Do ist=1,ic Do imon=1,Nmon irec=irec+1 ; tmp = Y(ist,imon) ; write (33,rec=irec) tmp Enddo Enddo Close (33) Open (4, file=trim(Out_Dir)//trim(Grads_File)//"_monthly.ctl") write (4,"(A)") "dset ^"//trim(Grads_File)//"_monthly.dat" write (4,"(A)") "title Monthly values" write (4,"(A,F10.3)") "undef ",RMiss ! write (4,"(A)") "options yrev" write (4,"(A,I3,A)") "xdef ",Nmon," linear 1 1" write (4,"(A,I3,A)") "ydef ",ic," linear 1 1" write (4,"(A)") "zdef 1 linear 1000 1" write (4,"(A,I4,A,I4,A)") "tdef ",1," linear 00z01Jan",1960," 1mon" write (4,"(A,I3)") "vars ", 1 ! ic write (4,"(A)") trim(Field)//" 0 99 Variable" write (4,"(A)") "endvars" Close(4) End Subroutine Monthly_to_Grads !------------------------------------------------------------- Subroutine Daily_to_Grads (X, Nst, Nmon, Grads_File, Out_Dir, Field) Integer :: Nst, Nmon, ist, imon, irec, ic, IDX(Nst) Character (Len=*) :: Grads_File, Out_Dir, Field !Real :: X(Nst, Nmon), Y(Nst, Nmon), Tmp Real :: X(Nst, NDayOfYear), Y(Nst, NDayOfYear), Tmp Logical :: Empty Do ist = 1, Nst if (.not.Empty_Arr1 (X(ist,1:Nmon),Nmon)) then IDX(ist)=1 ! print*,"IDX ",ist,"=",IDX(ist)," ", Station(ist) ! print"(12F7.1)",(X(ist,imon),imon=1,Nmon) else IDX(ist)=0 endif Enddo Open (33,file=trim(Out_Dir)//trim(Grads_File)//'_daily.dat', & FORM='Binary',access="direct",recl=4, status="unknown") ! write (99,"(A)") "Ylabs for Monthly to Grads" write (99,"(A)") trim(Field) ! Ten truong trong file Ylabs cho grads write (99,"(A)",Advance="No") '"' ! Dau " o dau dong ! print*,trim(Field) !! irec=0; ic = 0 Do ist = Nst, 1, -1 ! 1,Nst >>> South to North if (IDX(ist) == 1) then if (Leng_Series (X(ist,1:Nmon),Nmon) >= DayYear) then write (99,"(A)",Advance="No") trim(Station(ist))//"|" ! Ylabs ic = ic+1; Y(ic,1:Nmon)=X(ist,1:Nmon) ! print*,ic," ",trim(Station(ist)) endif endif Enddo write (99,"(A,I7)") '"',ic ! Ylabs Do ist=1,ic Do imon=1,Nmon irec=irec+1 ; tmp = Y(ist,imon) ; write (33,rec=irec) tmp Enddo Enddo Close (33) Open (4, file=trim(Out_Dir)//trim(Grads_File)//"_daily.ctl") write (4,"(A)") "dset ^"//trim(Grads_File)//"_daily.dat" write (4,"(A)") "title Daily values" write (4,"(A,F10.3)") "undef ",RMiss ! write (4,"(A)") "options yrev" write (4,"(A,I3,A)") "xdef ",Nmon," linear 1 1" write (4,"(A,I3,A)") "ydef ",ic," linear 1 1" write (4,"(A)") "zdef 1 linear 1000 1" write (4,"(A,I4,A,I4,A)") "tdef ",1," linear 00z01Jan",1960," 1mon" write (4,"(A,I3)") "vars ", 1 ! ic write (4,"(A)") trim(Field)//" 0 99 Variable" write (4,"(A)") "endvars" Close(4) End Subroutine Daily_to_Grads !---------------------------------------------------------------------------------- Subroutine Write_GrADS (X, Grads_File, Out_Dir, Field) Integer :: ivar, ist, ic Character (Len=8) :: STID Character *(*) :: Grads_File, Out_Dir, Field Real :: X(Nst, Nidx) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG ! Character (Len=2) :: ifile ! Character (Len=255) :: Grads_File !---------------------------------------------------------------------------------- ! OPEN (10,file=trim(Out_Dir)//trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') OPEN (10,file=trim(Out_Dir)//trim(Grads_File)//'.dat',FORM='UNFORMATTED') 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 ! ! print*,"Nst=",Nst Do ivar = 1,Nidx ! Nvar chi so, moi chi so duoc xem nhu 1 timeslice ! write (ifile,"(I2.2)") ivar; Grads_File=trim(Gr_File)//"_"//ifile ! OPEN (10,file=trim(Out_Dir)//trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') NLEV = 1 ! =0 o ban ghi cuoi cung cua 1 timeslice; =1 cho tat ca cac timeslice khac ic = 0 Do ist=1, Nst ! if (X(ist,ivar) /= RMiss) then ic = ic+1 STID = Station(ist)(1:8) ! Chi lay 8 ky tu dau tien RLAT = Lat(ist); RLON = Lon(ist) ! print*,STID WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Var = X(ist,ivar) !if (ivar==4.and.Var>50) print"(A20,A,F7.1)",trim(station(ist))," ",Var WRITE (10) Var ! Gia tri cua bien ! endif Enddo ! before new time: NLEV = 0 WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG ! print*,"Grasd=", ic Enddo ! ivar Close(10) !---------------------------------------------------------------------! ! Write control file Station data ! !---------------------------------------------------------------------! open (2001,file=trim(Out_Dir)//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 ',RMiss write (2001,'(A)') 'title Factor Analysis' write (2001,'(A,I3,1x,A)') 'tdef ',Nidx,'linear jan1980 1mo' ! write (2001,'(A,I3,1x,A)') 'tdef ',1,'linear jan1980 1mo' write (2001,'(A)') 'vars 1' write (2001,'(A)') trim(Field)//' 0 99 Factor' write (2001,'(A)') 'endvars' close (2001) !Enddo ! ivar End Subroutine Write_GrADS !------------------------------------------------------------------------------------- Integer Function Leng_Series (W,N) Integer :: N, i, ic Real :: W(N) ic=0 Do i=1,N if (W(i) /= RMiss) ic=ic+1 Enddo Leng_Series = ic End Function Leng_Series !------------------------------------------------------------------------------------- Subroutine Hovmoller (X, Grads_File, Out_Dir, Field) Integer :: ist1, ist, iyr, NYr, NYr1, irec, Nc, ic, IDX(Nst), NY Character (Len=*) :: Grads_File, Out_Dir, Field Real :: X(Nst, Yr1:Yr2), Tmp, TB(Nst) Real, Allocatable :: W(:), Y(:,:) Logical :: Empty ! Define Stations that data are available IDX(:) = 0; Nc = 0 ; NYr1 = Yr2-Yr1+1; TB(:) = RMiss Do ist = 1, Nst Empty = .true.; ic = 0 Do iyr = Yr1, Yr2 if (X(ist,iyr) /= RMiss) then ! Empty = .false.; exit ic = ic + 1 ! Dem so nam co SL endif Enddo ! if (real(ic)/real(NYr1) >= 0.7) Empty = .false. ! So nam co SL duoi 70% if (ic >=NumYear) Empty = .false. if (.not.Empty) then IDX(ist) = 1; Nc = Nc+1 endif Enddo if (Nc <=0 ) return ! Khong co SL de ve !if (trim(Grads_File)=="Eva") print*,IDX NYr = Yr2-Yr1+1 Allocate (W(1:NYr), Y(1:NYr,1:Nc)) !, LonC(1:Nc), LatC(1:Nc)) Open (33,file=trim(Out_Dir)//trim(Grads_File)//'_Hovmoller.dat',FORM='UNFORMATTED', & access="direct",recl=4*NYr*Nc,status="unknown") write (99,"(A)") "Ylabs for Hovmoller Diagram" write (99,"(A)") trim(Field) ! Ten truong trong file Ylabs cho grads write (99,"(A)",Advance="No") '"' ! Dau " o dau dong ic = 0 Loop_Nst: Do ist = Nst,1,-1 ! 1, Nst >>> South to North if (IDX(ist) == 0) Cycle Loop_Nst !-------------------------------------------------- ! write (99,"(A)",Advance="No") trim(Station(ist))//"|" ! Ylabs NY = 0; tmp = 0. Do iyr=1980, 1999 if (X(ist,iyr) /= RMiss) then NY = NY+1; tmp = tmp+X(ist,iyr) endif Enddo if (NY>=16) TB(ist)=tmp/Real(NY) ! print"(I4,A,I4,F10.1)",IDX(ist),Station(ist),NY, TB(ist) Do iyr=Yr1, Yr2 if (X(ist,iyr)/=RMiss.AND.TB(ist)/=RMiss) then; X(ist,iyr) = X(ist,iyr) - TB(ist); else; X(ist,iyr) = RMiss; endif Enddo !-------------------------------------------------- ! if (Leng_Series(X(ist,Yr1:Yr2),Yr2-Yr1+1) >= NumYear.AND.Station(ist)(1:2) /= "P-") then if (Leng_Series(X(ist,Yr1:Yr2),Yr2-Yr1+1) >= NumYear) then ic = ic+1; Y(1:NYr,ic) = X(ist,Yr1:Yr2) write (99,"(A)",Advance="No") trim(Station(ist))//"|" ! Ylabs endif Enddo Loop_Nst write (99,"(A,I7)") '"',ic ! Ylabs write (33,rec=1) Y Close (33) Open (4, file=trim(Out_Dir)//trim(Grads_File)//"_Hovmoller.ctl") write (4,"(A)") "dset ^"//trim(Grads_File)//"_Hovmoller.dat" write (4,"(A)") "title Hovmoller Diagram" write (4,"(A,F10.3)") "undef ",RMiss write (4,"(A,I4,A,I4,A)") "xdef ",Nyr," linear ",Yr1," 1" write (4,"(A,I4,A)") "ydef ",ic," linear 1 1" ! write (4,"(A)") "options yrev" write (4,"(A)") "zdef 1 linear 1000 1" write (4,"(A,I4,A,I4,A)") "tdef ",1," linear 00z01Jan",Yr1," 1yrs" write (4,"(A,I3)") "vars 1" write (4,"(A,A,A)") trim(Field)," 0 99 "," Variable " write (4,"(A)") "endvars" Close(4) End Subroutine Hovmoller !------------------------------------------------------------- Subroutine R_Season (MDry1, MDry2, MWet1, MWet2) ! RR01_Dry (Nst,Yr1:Yr2) : So ngay mua trong mua kho ! RR01_Wet (Nst,Yr1:Yr2) : So ngay mua trong mua mua ! R_Dry (Nst,Yr1:Yr2) : Tong luong mua trong mua kho ! R_Wet (Nst,Yr1:Yr2) : Tong luong mua trong mua mua ! Rmon(ist, iyr, imon) : Luong mua cac thang ! NR01(ist, iyr, imon) : So ngay mua cac thang ! MDry1, MDry2, MWet1, MWet2 : Thang bat dau, ket thuc mua kho, mua mua Integer :: ic1, ic2, ic3, ic4, imon, NDryMon, NWetMon Integer :: MDry1, MDry2, MWet1, MWet2 Real :: Tmp1, Tmp2, Tmp3, Tmp4 RR01_Dry(:,:) = RMiss; RR01_Wet(:,:) = RMiss; R_Dry(:,:) = RMiss; R_Wet(:,:) = RMiss Eva_Dry(:,:) = RMiss; Eva_Wet(:,:) = RMiss NDryMon = Nmon-MDry1+1+MDry2 ; NWetMon = MWet2-MWet1+1 Do ist=1, Nst Do iyr=Yr1, Yr2-1 ! Mua kho dau tien tinh cho nam dau, mua kho cuoi cung la nam Yr2-1 ic1 = 0; Tmp1 = 0; ic2 = 0; Tmp2 = 0; ic3 = 0; Tmp3 = 0 Do imon=1, Nmon if (imon >= MDry1) then if (NR01(ist,iyr,imon) /= RMiss) then ! So ngay mua ic1 = ic1+1; Tmp1 = Tmp1+NR01(ist,iyr,imon) endif if (Rmon(ist,iyr,imon) /= RMiss) then ! Luong mua ic2 = ic2+1; Tmp2 = Tmp2+Rmon(ist,iyr,imon) endif if (Eva_mon(ist,iyr,imon) /= RMiss) then ! Luong BH ic3 = ic3+1; Tmp3 = Tmp3+Eva_mon(ist,iyr,imon) endif else if (imon <= MDry2) then if (NR01(ist,iyr+1,imon) /= RMiss) then ic1 = ic1+1; Tmp1 = Tmp1+NR01(ist,iyr+1,imon) endif if (Rmon(ist,iyr+1,imon) /= RMiss) then ic2 = ic2+1; Tmp2 = Tmp2+Rmon(ist,iyr+1,imon) endif if (Eva_mon(ist,iyr+1,imon) /= RMiss) then ic3 = ic3+1; Tmp3 = Tmp3+Eva_mon(ist,iyr+1,imon) endif endif Enddo if (ic1>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp1 = Tmp1/Real(ic1); RR01_Dry(ist,iyr) = Tmp1*NDryMon ! Tong mua kho endif if (ic2>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp2 = Tmp2/Real(ic2); R_Dry(ist,iyr) = Tmp2*NDryMon ! Tong mua kho endif if (ic3>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp3 = Tmp3/Real(ic3); Eva_Dry(ist,iyr) = Tmp3*NDryMon ! Tong mua kho endif Enddo ! ------------- Wet Season ----------- Do iyr=Yr1, Yr2 ! Mua mua tinh cho tat ca cac nam ic1 = 0; Tmp1 = 0; ic2 = 0; Tmp2 = 0; ic3 = 0; Tmp3 = 0 Do imon=MWet1, MWet2 if (NR01(ist,iyr,imon) /= RMiss) then ic1 = ic1+1; Tmp1 = Tmp1+NR01(ist,iyr,imon) endif if (Rmon(ist,iyr,imon) /= RMiss) then ! Luong mua ic2 = ic2+1; Tmp2 = Tmp2+Rmon(ist,iyr,imon) endif if (Eva_mon(ist,iyr,imon) /= RMiss) then ! Luong BH ic3 = ic3+1; Tmp3 = Tmp3+Eva_mon(ist,iyr,imon) endif Enddo if (ic1>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp1 = Tmp1/Real(ic1); RR01_Wet(ist,iyr) = Tmp1*NWetMon ! Tong mua mua endif if (ic2>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp2 = Tmp2/Real(ic2); R_Wet(ist,iyr) = Tmp2*NWetMon ! Tong mua mua endif if (ic3>=NDryMon/3*2) then ! Du 2/3 so thang co SL Tmp3 = Tmp3/Real(ic3); Eva_Wet(ist,iyr) = Tmp3*NWetMon ! Tong mua mua endif Enddo Enddo ! Tinh Trung binh khi hau cua R va BH tron mua mua va mua kho Do ist=1, Nst ic1 = 0; Tmp1 = 0.; ic2 = 0; Tmp2 = 0.; ic3 = 0; Tmp3 = 0.; ic4 = 0; Tmp4 = 0. Do iyr=Yr1, Yr2 if (R_Dry(ist,iyr) /= RMiss) then ic1 = ic1+1; Tmp1 = Tmp1 + R_Dry(ist,iyr) endif if (Eva_Dry(ist,iyr) /= RMiss) then ic2 = ic2+1; Tmp2 = Tmp2 + Eva_Dry(ist,iyr) endif if (R_Wet(ist,iyr) /= RMiss) then ic3 = ic3+1; Tmp3 = Tmp3 + R_Wet(ist,iyr) endif if (Eva_Wet(ist,iyr) /= RMiss) then ic4 = ic4+1; Tmp4 = Tmp4 + Eva_Wet(ist,iyr) endif Enddo if (ic1>=NumYear) then; R_Dry_Cli(ist) = Tmp1/real(ic1); else; R_Dry_Cli(ist) = RMiss; endif if (ic2>=NumYear) then; Eva_Dry_Cli(ist) = Tmp2/real(ic2); else; Eva_Dry_Cli(ist) = RMiss; endif if (ic3>=NumYear) then; R_Wet_Cli(ist) = Tmp3/real(ic3); else; R_Wet_Cli(ist) = RMiss; endif if (ic4>=NumYear) then; Eva_Wet_Cli(ist) = Tmp4/real(ic4); else; Eva_Wet_Cli(ist) = RMiss; endif Enddo End Subroutine R_Season !----------------------------------------------------------------------- Subroutine Linear_Trend(X,Y,N,A1,F_F10,RMiss) Integer :: N, i, j, N1 Real :: X(N), Y(N), W1(N), W2(N), TBx, TBy, R, A0, A1, F10, F, U, Q, Yhq, F_F10 Real :: Ws1(N), Ws2(N) Real :: Sx, Sy, RR, RMiss !print*,Y TBx = 0; TBy = 0; RR = 0; N1 = 0; Do i=1,N if (Y(i) /= RMiss) then N1 = N1 + 1 TBx = TBx + X(i); TBy = TBy + Y(i) Sx = Sx + X(i)*X(i); Sy = Sy + Y(i)*Y(i) RR = RR + X(i)*Y(i) W1(N1) = X(i); W2(N1) = Y(i) endif Enddo if ((N1 < NumYear).OR.((TBx == 0).AND.(TBy == 0)).OR.(Sx==0).OR.(Sy==0)) then A0 = -99; A1 = -99; R = -99; F_F10 = -99; else TBx = TBx/Real(N1); TBy = TBy/Real(N1) Sx = SQRT(Sx/Real(N1)-TBx*TBx); Sy = SQRT(Sy/Real(N1)-TBy*TBy) R = (RR/Real(N1)-TBx*TBy)/(Sx*Sy) A1 = R*Sy/Sx A0 = TBy - A1*TBx F10 = FINV(0.10,1.0,Real(N1)) U = 0; Q = 0 Do i=1,N1 Yhq = A0 + A1*W1(i) U = U + (Yhq-TBy)*(Yhq-TBy) Q = Q + (W2(i)-Yhq)*(W2(i)-Yhq) Enddo if (Q == 0) then F = -99 else F = U/(Q/(N1-2)) endif if (F /=-99) then; F_F10 = F-F10; else; F_F10 = -99; endif endif End subroutine Linear_Trend Real FUNCTION FINV(P,N1,N2) Real, PARAMETER :: EPS=1.0E-6 REAL :: P,AP,N1,N2, A,B,C,P0,SS,SS1 IF (P.LT.0.0.OR.P.GT.1.0) THEN WRITE(*,*)' INVALID NUMERIC INPUT IN FINV FUNCTION' STOP ENDIF IF (P.EQ.0.0) THEN WRITE(*,*)' THE P VALUE TOO SMALL IN FINV FUNCTION' STOP ELSE IF (P.EQ.1) THEN FINV=0.0 RETURN ENDIF ! A=0.0 B=9999999.0 C=A AP=P 10 P0=FDIST(B,N1,N2) SS=ABS((P0-AP)/AP) SS1=ABS((B-C)/B) IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN IF (AP.GT.P0) THEN C=B B=(A+B)/2.0 GOTO 10 ELSE IF (AP.LT.P0) THEN B=(C+B)/2.0 GOTO 10 ENDIF ENDIF A=(C+B)/2.0 FINV=A END FUNCTION FINV Real FUNCTION FDIST(X0,N1,N2) REAL :: X0,N1,N2,X X=N2/(N2+N1*X0) FDIST=BETAI(0.5*N2,0.5*N1,X) END FUNCTION FDIST Real FUNCTION betai(a,b,x) REAL :: a,b,x REAL :: bt if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai' if(x.eq.0..or.x.eq.1.)then bt=0. else bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x)) endif if(x.lt.(a+1.)/(a+b+2.))then betai=bt*betacf(a,b,x)/a return else betai=1.-bt*betacf(b,a,1.-x)/b return endif END FUNCTION betai Real FUNCTION betacf(a,b,x) INTEGER, PARAMETER :: MAXIT=100 REAL :: a,b,x Real, PARAMETER :: EPS=3.e-7,FPMIN=1.e-30 INTEGER :: m,m2 REAL :: aa,c,d,del,h,qab,qam,qap qab=a+b qap=a+1. qam=a-1. c=1. d=1.-qab*x/qap if(abs(d).lt.FPMIN)d=FPMIN d=1./d h=d do 11 m=1,MAXIT m2=2*m aa=m*(b-m)*x/((qam+m2)*(a+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d h=h*d*c aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue pause 'a or b too big, or MAXIT too small in betacf' 1 betacf=h END FUNCTION betacf Real FUNCTION gammln(xx) REAL :: xx INTEGER :: j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & -.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) END FUNCTION gammln !----------------------------------------------------------------------- Subroutine MK_Sen(Y,N,Miss,SenSlope,prob,N1) Integer :: N,i,j,S, Ng, T(N), N1 !Real, Intent (In) :: X(N) Real :: X(N), Y(N) Real :: VarS, z, Miss, prob, SenSlope, Q(N*(N-1)/2) Logical :: Empty X(:) = Y(:) ! Dam bao mang dau vao khong bi thay doi Empty = .true.; N1=0 Do i=1, N if (X(i) /= Miss) then N1 = N1+1 endif Enddo if (N1 >= NumYear) Empty = .false. if (Empty) then SenSlope = Miss; prob = Miss ; Return else ! Calculate Sen's slope first Ng = 0 Do i=1, N-1 Do j=i+1, N if (X(j) /= Miss .And. X(i) /= Miss) then Ng = Ng+1 Q(Ng) = (X(j)-X(i))/Real(j-i) endif Enddo Enddo Call My_Sort (Q, Ng) if (mod(Ng,2) /= 0) then SenSlope = Q(Ng/2+1) else SenSlope = (Q(Ng/2) + Q(Ng/2+1))/2. endif S = 0 Do i=1,N-1 Do j=i+1,N if (X(j) /= Miss .And. X(i) /= Miss) then if (X(j) > X(i)) then S = S+1 else if (X(j) < X(i)) then S = S-1 endif endif Enddo Enddo Call My_Sort(X, N) Ng = 0; T(:) = 0; i=0 Do if (i>=N) exit i = i+1 if (i == 1) then Ng = Ng+1; T(Ng) = T(Ng)+1 else if (X(i)==X(i-1)) then T(Ng) = T(Ng)+1 else Ng = Ng+1; T(Ng) = T(Ng)+1 endif endif Enddo VarS = 0 if (X(1) == Miss) then ! Loai nhom dau tien neu co Miss Do i=2, Ng VarS = VarS + T(i)*(T(i)-1)*(2*T(i)+5) Enddo else Do i=1, Ng VarS = VarS + T(i)*(T(i)-1)*(2*T(i)+5) Enddo endif VarS = (N*(N-1)*(2*N+5)-VarS)/18. if (S > 0) then z = (S-1)/SQRT(VarS) else if (S < 0) then z = (S+1)/SQRT(VarS) else z = 0 endif prob = erfcc(abs(z)/1.4142136) if (prob>=1.) then !print "(2(A,F10.3))","S=",SenSlope," Pro=",prob SenSlope = Miss; prob = Miss endif endif End Subroutine MK_Sen Subroutine My_Sort(X,N) ! Sap xep X(n) thanh chuoi tang dan Integer :: N, i,j Real :: X(N), tmp 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 End Subroutine My_Sort Real FUNCTION erfcc (x) REAL :: x, tmp REAL :: t,z z = abs(x) t = 1./(1.+0.5*z) tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ & t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ & t*(1.48851587+t*(-.82215223+t*.17087277))))))))) if (x.lt.0.) tmp=2.-tmp erfcc=tmp END FUNCTION erfcc !----------------------------------------------------------------------- Real function percentile(x, length, per) integer length real x(length), per real xtos(length),bb,cc,oout integer nn nn=0 do i=1, length if(x(i) /= RMiss)then nn=nn+1 xtos(nn)=x(i) endif enddo if(nn.eq.0) then oout=RMiss else call sort(nn,xtos) bb=nn*per+per/3.+1/3. cc=real(int(bb)) if(int(cc).ge.nn) then oout=xtos(nn) else oout=xtos(int(cc))+(bb-cc)*(xtos(int(cc)+1)-xtos(int(cc))) endif endif percentile = oout end function percentile !----------------------------------------------------------------------- SUBROUTINE sort(n,arr) INTEGER n,M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,ir,j,jstack,k,l,istack(NSTACK) REAL a,temp jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 12 j=l+1,ir a=arr(j) do 11 i=j-1,1,-1 if(arr(i).le.a)goto 2 arr(i+1)=arr(i) 11 continue i=0 2 arr(i+1)=a 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp if(arr(l+1).gt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp endif if(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp endif if(arr(l+1).gt.arr(l))then temp=arr(l+1) arr(l+1)=arr(l) arr(l)=temp endif i=l+1 j=ir a=arr(l) 3 continue i=i+1 if(arr(i).lt.a)goto 3 4 continue j=j-1 if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp goto 3 5 arr(l)=arr(j) arr(j)=a jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in sort' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END SUBROUTINE sort