Program Preprocessing Integer, Parameter :: NVar = 2, Nst_Max = 1000 ! So Grid Cells hoac so tram Integer, Parameter :: NMonth = 12, NWeek = 52 Integer, Parameter :: YBas1 = 1981, YBas2 = 2014, NYear = YBas2-YBas1+1, maxyrs = NYear Character (Len=4) :: YT(NVar) = (/"T2m","R"/) Integer, Parameter :: NPDI = 3 ! 1=PDSI, 2= Character (Len=21) :: Station(Nst_Max) Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NWeek) :: Dat_Mon, Dat_Week ! ivar=1: T2m ivar=2: Rain Real, Dimension (NVar, Nst_Max, NWeek) :: DatNor ! Normal Climate Condition Real :: AWC(Nst_Max), xlon(Nst_Max), xlat(Nst_Max), HH(Nst_Max) ! Outputs: Real, Dimension (Nst_Max,NYear*NWeek) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind Real, Dimension (Nst_Max,NYear*NWeek) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind Real, Dimension (Nst_Max,NYear*NWeek) :: wee_1_pdsi, wee_1_wplm, wee_1_zind Real, Dimension (Nst_Max,NYear*NWeek) :: wee_2_pdsi, wee_2_wplm, wee_2_zind Real, Dimension (Nst_Max,NYear*NWeek) :: Rain_Mon, Rain_Wee ! Integer, Parameter :: Nidx=14, Ntq = 7 ! Kich thuoc Ma tran Tuong quan (1..3: Mon_Ori; 4..6: Mon_Sel; 7: Rain Real, Dimension (Nst_Max,NYear*NWeek,Ntq) :: PMon, PWee ! Luu tat ca cac chuoi PDSI de tinh HSTQ Real, Dimension (Nst_Max,Ntq,Ntq) :: TQ_M, TQ_W ! Ma tran HSTQ Mon & Week Real, Dimension (Nst_Max,NYear*NWeek,Nidx) :: All_IDX ! Tat ca cac Chi so Real, Dimension (Nidx,Nst_Max,Nst_Max) :: TQ_Stn ! Tuong quan giua cac tram cua cac Chi so Real, Parameter :: NThre = 5 Real, Parameter :: Thre(NThre) = (/-0.5, -1.0, -2.0, -3.0, -4.0/) Integer, Dimension (Nst_Max,NYear*NWeek,Nidx,NThre) :: Leng_Spel ! So dot han ung voi do dai cac dot Integer, Dimension (Nst_Max,YBas1:YBas2,Nidx,NThre) :: Drought ! So thang/tuan han trong tung nam Real, Dimension (Nst_Max,NWeek,Nidx,NThre) :: Dr_MW,Dr_MW0 ! So nam han tung thang/tuan Integer, Dimension (Nst_Max,Nidx,NThre) :: NSpell ! Do dai cac dot han ung voi cac cap han !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, Dat_Dir, Res_Dir Character (Len=30) :: FMT Integer :: Nst, IO, Ntime Integer :: ist, iyr, imw, ivar, isys, system !------------------------------------------------------------------------------ Dat_Mon(:,:,:,:) = -99.; Dat_Week(:,:,:,:) = -99.; DatNor(:,:,:) = 0.; Rain_Mon(:,:)=-99; Rain_Wee(:,:)=-99 All_IDX(:,:,:) = -99; Leng_Spel(:,:,:,:) = -99; Drought(:,:,:,:) = -99 !BL_Dir = "/home/tanpv/PDSI/PDSI_2015-Oct/Data/" BL_Dir = "/work/projects/G09/OBS_VN/Monthly/TXT_R_T2m_4_Drought/" Res_Dir = "/home/tanpv/PDSI/PDSI/my_result/" Print*, "Reading Meteo data and station info ...." Call Read_data_Baseline("_Mon_from_Day_4_Drought_1981-2014.txt",Dat_Mon,NMonth) Call Read_data_Baseline("_Week_from_Day_4_Drought_1981-2014.txt",Dat_Week,NWeek) Do ist=1, Nst Ntime=0 Do iyr=YBas1, YBas2 Do imw=1, NWeek Ntime=Ntime+1 Rain_Mon(ist,Ntime) = Dat_Mon(2,ist,iyr,imw); Rain_Wee(ist,Ntime) = Dat_Week(2,ist,iyr,imw) Enddo Enddo Enddo Call putAWC(AWC, xlon, xlat, Nst) isys = system ("rm -Rf /home/tanpv/PDSI/PDSI/my_result") mon_ori_pdsi(:,:)=-99; mon_ori_wplm(:,:)=-99; mon_ori_zind(:,:)=-99; mon_sel_pdsi(:,:)=-99; mon_sel_wplm(:,:)=-99; mon_sel_zind(:,:)=-99; wee_1_pdsi(:,:)=-99; wee_1_wplm(:,:)=-99; wee_1_zind(:,:)=-99; wee_2_pdsi(:,:)=-99; wee_2_wplm(:,:)=-99; wee_2_zind(:,:)=-99; Do ist=1, Nst Print*, "Calculation for station: "//trim(station(ist)) Call Write_to_PDSI ("monthly", Dat_Mon, NMonth, ist) Call Write_to_PDSI ("weekly" , Dat_Week, NWeek, ist) isys = system ("./pdsi") ! PostProcessing for station Call Read_Table (trim(Res_Dir)//"monthly/original/"//"PDSI.tbl", NMonth, mon_ori_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/original/"//"WPLM.tbl", NMonth, mon_ori_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/original/"//"ZIND.tbl", NMonth, mon_ori_zind(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"PDSI.tbl", NMonth, mon_sel_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"WPLM.tbl", NMonth, mon_sel_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"ZIND.tbl", NMonth, mon_sel_zind(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/1/"//"PDSI.tbl", NWeek, wee_1_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/1/"//"WPLM.tbl", NWeek, wee_1_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/1/"//"ZIND.tbl", NWeek, wee_1_zind(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/2/"//"PDSI.tbl", NWeek/2, wee_2_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/2/"//"WPLM.tbl", NWeek/2, wee_2_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"weekly/2/"//"ZIND.tbl", NWeek/2, wee_2_zind(ist,:)) PMon(ist,:,1)=mon_ori_pdsi(ist,:); PMon(ist,:,2)=mon_ori_wplm(ist,:); PMon(ist,:,3)=mon_ori_zind(ist,:) PMon(ist,:,4)=mon_sel_pdsi(ist,:); PMon(ist,:,5)=mon_sel_wplm(ist,:); PMon(ist,:,6)=mon_sel_zind(ist,:) PMon(ist,:,7)=Rain_Mon(ist,:) PWee(ist,:,1)=wee_1_pdsi(ist,:); PWee(ist,:,2)=wee_1_wplm(ist,:); PWee(ist,:,3)=wee_1_zind(ist,:) PWee(ist,:,4)=wee_2_pdsi(ist,:); PWee(ist,:,5)=wee_2_wplm(ist,:); PWee(ist,:,6)=wee_2_zind(ist,:) PWee(ist,:,7)=Rain_Wee(ist,:) Call HSTQ (PMon, TQ_M, ist); Call HSTQ (PWee, TQ_W, ist) All_IDX(ist,:,1)=mon_ori_pdsi(ist,:); All_IDX(ist,:,2)=mon_ori_wplm(ist,:); All_IDX(ist,:,3)=mon_ori_zind(ist,:) All_IDX(ist,:,4)=mon_sel_pdsi(ist,:); All_IDX(ist,:,5)=mon_sel_wplm(ist,:); All_IDX(ist,:,6)=mon_sel_zind(ist,:) All_IDX(ist,:,7)=wee_1_pdsi(ist,:); All_IDX(ist,:,8)=wee_1_wplm(ist,:); All_IDX(ist,:,9)=wee_1_zind(ist,:) All_IDX(ist,:,10)=wee_2_pdsi(ist,:); All_IDX(ist,:,11)=wee_2_wplm(ist,:); All_IDX(ist,:,12)=wee_2_zind(ist,:) All_IDX(ist,:,13)=Rain_Mon(ist,:); All_IDX(ist,:,14)=Rain_Wee(ist,:) Enddo ! Print Out all Results: Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Ori_PDSI.txt", mon_ori_pdsi, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Ori_WPLM.txt", mon_ori_wplm, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Ori_ZIND.txt", mon_ori_zind, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Sel_PDSI.txt", mon_ori_pdsi, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Sel_WPLM.txt", mon_ori_wplm, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Mon_Sel_ZIND.txt", mon_ori_zind, Nst, NMonth) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_1_PDSI.txt", wee_1_pdsi, Nst, NWeek) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_1_WPLM.txt", wee_1_wplm, Nst, NWeek) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_1_ZIND.txt", wee_1_zind, Nst, NWeek) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_2_PDSI.txt", wee_2_pdsi, Nst, NWeek/2) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_2_WPLM.txt", wee_2_wplm, Nst, NWeek/2) Call Write_Table ("/home/tanpv/PDSI/PDSI/Results/Wee_2_ZIND.txt", wee_2_zind, Nst, NWeek/2) Call Write_HSTQ ("HSTQ_Mon.txt", TQ_M); Call Write_HSTQ ("HSTQ_Wee.txt", TQ_W) Call Calc_Write_HSTQ1 Call Calc_Spell Call Hovmoller_BL ("All_IDX_M", "../4Grads/", "Month") Call Hovmoller_BL ("All_IDX_W1", "../4Grads/", "Week1") Call Hovmoller_BL ("All_IDX_W2", "../4Grads/", "Week2") CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_Spell Integer :: N1, ist, idx, iyr, imw, it, ith, Nspl Character (Len=10) :: CapHan(Nthre) = (/"Chom Han","Han nhe","Han vua","Han nang","Rat nang"/) Character (Len=3) :: IDXs(Nidx)=(/"mop", "mow", "moz", "msp", "msw", "msz", "w1p", & "w1w", "w1z", "w2p", "w2w", "w2z", "mpr", "wpr"/) Integer :: w1(NYear*NWeek), w2(NYear*NWeek), NSpell_Max, Nmw, ic1, ic2 Logical :: LDrou(NThre) !------------------------------- Real :: Trend_Dr(Nst_Max,Nidx,NThre), Prob ! He so goc Sen cua cac cap han cua tung tram theo cac chi so Real :: wrk(NYear) !------------------------------- !Integer, Dimension (Nst_Max,NYear*NWeek,Nidx,NThre) :: Leng_Spel, Drought ! =1: Yes, =0: No ! Nst_Max,NYear,Nidx,NThre !Integer, Dimension (Nst_Max,Nidx,NThre) :: NSpell !Real, Dimension (Nst_Max,NYear*NWeek,Nidx) :: All_IDX ! Tinh trung binh cac chi so cho tat ca cac tram, gan gia tri vao tram thu (Nst+1) Do idx=1, Nidx-2 if (idx <= (Nidx-2)/2) then Nmw = Nmonth else if (idx <= (Nidx-2)/4*3) then Nmw = NWeek else Nmw = NWeek/2 endif Do it=1, NYear*Nmw All_IDX(Nst+1,it,idx) = 0 Do ist=1, Nst All_IDX(Nst+1,it,idx) = All_IDX(Nst+1,it,idx) + All_IDX(ist,it,idx) Enddo All_IDX(Nst+1,it,idx) = All_IDX(Nst+1,it,idx)/Nst Enddo Enddo Dr_MW(:,:,:,:) = 0; Dr_MW0(:,:,:,:) = 0 Do ist=1, Nst + 1 Do ith=1, NThre Do idx=1, Nidx-2 ! Tru hai chi so cuoi la Luong mua if (idx <= (Nidx-2)/2) then Nmw = Nmonth else if (idx <= (Nidx-2)/4*3) then Nmw = NWeek else Nmw = NWeek/2 endif it = 0; w1(:)=-99; w2(:)=-99 Do iyr=YBas1, YBas2 Drought(ist,iyr,idx,ith) = 0 Do imw=1, Nmw ! NWeek it = it+1; LDrou(:) = .False. ! it = 1,..., NYear * Nmw if (All_IDX(ist,it,idx) /= -99) then if (ith < NThre) then if (All_IDX(ist,it,idx)<=Thre(ith) .AND. All_IDX(ist,it,idx)>Thre(ith+1)) LDrou(ith)=.True. else if (All_IDX(ist,it,idx)<=Thre(ith)) LDrou(ith)=.True. endif if (LDrou(ith)) then ! (All_IDX(ist,it,idx) <= Thre(ith)) then Drought(ist,iyr,idx,ith) = Drought(ist,iyr,idx,ith)+1 w1(it) = 1 else w1(it) = 0 endif else w1(it) = -99 endif Enddo Enddo Call NSpell_and_Leng (w1, NYear*Nmw, w2, Nspl) ! w1(NYear*Nmw) ! Leng_Spel(ist,1:NYear*NWeek,idx,ith) = w2(1:NYear*NWeek); NSpell(ist,idx,ith) = Nspl Leng_Spel(ist,1:NYear*Nmw,idx,ith) = w2(1:NYear*Nmw); NSpell(ist,idx,ith) = Nspl ! Do it=1, NYear*Nmw !NYear*NWeek ! imw=1, Nmw LDrou(:) = .False. imw = mod(it, Nmw); if (imw == 0) imw = Nmw if (All_IDX(ist,it,idx) /= -99) then if (ith < NThre) then if (All_IDX(ist,it,idx)<=Thre(ith) .AND. All_IDX(ist,it,idx)>Thre(ith+1)) LDrou(ith)=.True. else if (All_IDX(ist,it,idx)<=Thre(ith)) LDrou(ith)=.True. endif if (LDrou(ith)) then ! (All_IDX(ist,it,idx) <= Thre(ith)) then ! if (All_IDX(ist,it,idx) <= Thre(ith)) then Dr_MW(ist,imw,idx,ith) = Dr_MW(ist,imw,idx,ith) + 1./NYear ! Tan suat han tung thang/tuan ! else ! Dr_MW0(ist,imw,idx,ith) = Dr_MW0(ist,imw,idx,ith) + 1!./NYear endif endif Enddo Enddo Enddo Enddo Station(Nst+1) = "Ave_Stn" ! Trung binh cac tram Open (3, file="Drought_Spells.txt", Status="Unknown") Open (4, file="Drought_Duration.txt", Status="Unknown") Open (2, file="Drought_Mon_Frac.txt", Status="Unknown") write (3, "(A)") " Cac dot han ung voi cac nguong" Do idx=1, Nidx-2 !NThre Do ith=1, NThre !Nidx-2 write (3, "(10x,4A)") "Cap Han: ",trim(CapHan(ith)), " - Tinh theo chi so: ", IDXs(idx) NSpell_Max = -9999 Do ist=1, Nst + 1 if (NSpell(ist,idx,ith) > NSpell_Max) NSpell_Max = NSpell(ist,idx,ith) Enddo write (3, "(A10,500I3)")" Length ",(ist, ist=1,NSpell_Max) Do ist=1, Nst + 1 if (NSpell(ist,idx,ith) >= 1) then write (3, "(A10,500I3)") trim(Station(ist)),(Leng_Spel(ist,it,idx,ith), it=1,NSpell_Max)! NSpell(ist,idx,ith)) endif Enddo Enddo write (3, "(A)") "----------------------------------------------------" Enddo !---------------------------------------------------- write (4,"(A)") " So thang/tuan han cua cac tram theo tung nam" write (2,"(A)") " Xac suat khi hau xay ra han trong tung thang/tuan" Do idx=1, Nidx-2 if (idx <= (Nidx-2)/2) then; Nmw = Nmonth; else; Nmw = NWeek; endif Do ith=1, NThre write (4, "(10x,4A)") "Cap Han: ",trim(CapHan(ith)), " - Tinh theo chi so: ", IDXs(idx) write (4, "(A5,100(1x,A7))") "Yrs", (trim(Station(ist)), ist=1, Nst+1) Do iyr=YBas1, YBas2 write (4, "(I4.4,1x,100(2x,I6))") iyr, (Drought(ist,iyr,idx,ith), ist=1, Nst+1) Enddo write (2, "(10x,4A)") "Cap Han: ",trim(CapHan(ith)), " - Tinh theo chi so: ", IDXs(idx) write (2, "(A5,100(1x,A7))") "M/W", (trim(Station(ist)), ist=1, Nst+1) Do imw=1, Nmw write (2, "(2x,I2.2,1x,100(1x,F7.1))") imw, (Dr_MW(ist,imw,idx,ith)*100,ist=1, Nst+1) ! write (2, "(2x,I2.2,1x,100(1x,F7.1))") imw, (Dr_MW0(ist,imw,idx,ith),ist=1, Nst+1) Enddo Enddo Enddo Close (3); Close(4); Close (2) ! Tinh xu the Do ist=1, Nst+1 Do idx=1, Nidx-2 Do ith=1, NThre wrk(1:NYear) = Drought(ist,YBas1:YBas2,idx,ith) !print "(10F8.1)", wrk(1:NYear) !pause Call MK_Sen(wrk,NYear,-99.,SenSlope,prob,N1) Trend_Dr(ist,idx,ith) = SenSlope*10. ! Per decade Enddo Enddo Enddo Open (3, file="Drought_Trend.txt", Status="Unknown") Do idx=1, Nidx-2 write (3, "(2A)") " Tinh theo chi so: ", IDXs(idx) write (3, "(A15,10A8)") "Stn\Dr.Lvs: ",(trim(CapHan(ith)), ith=1, NThre) Do ist=1, Nst+1 write (3,"(A15,10F8.3)") Station(ist), (Trend_Dr(ist,idx,ith), ith=1, NThre) Enddo Enddo Close (3) End Subroutine Calc_Spell !------------------------------------------------------------------------------ Subroutine NSpell_and_Leng (X, N, M, Imax) ! Input: Array X(N) with values of 0, 1, or -99 ! Output: Array M(Imax) with values of number of spell Integer :: N, i, k !Integer :: X(NYear*NWeek), M(NYear*NWeek), Imax INTEGER :: X(N), M(N), IMAX N = NYear*NWeek M(:) = 0; IMAX=0; I = 0 DO I = I + 1 ! Loop until end of series data IF (X(I)==1) THEN ! Event occur at first time of this spell K = 0 ! Number of days event occur continuosly DO IF (X(I)==0.OR.X(I)==-99.OR.I==N) THEN ! End of this spell IF (I==N.AND.X(I)==1) K=K+1 ! For the last spell IF (K/=0) M(K)=M(K)+1 ! Increase for this spell IMAX=IMAX+1 ! Increase number of spells ! PRINT*, K, M(K), IMAX EXIT END IF I = I + 1 K = K + 1 ! How many days of this spell ENDDO ELSE CYCLE END IF IF (I>=N) EXIT ENDDO I=IMAX DO IF (M(I)/=0) EXIT ! Determine exactly how many spells I=I-1 ENDDO IMAX = I !print*,"imax=",imax,"M=",M(1:imax) End Subroutine NSpell_and_Leng !------------------------------------------------------------------------------ Subroutine Write_HSTQ (Fname, X) Real, Dimension (Nst_Max,Ntq,Ntq) :: X Integer :: ist, j,k Character (Len=*) :: Fname Character (Len=2) :: Title(Ntq) = (/"op","ow","oz","sp","sw","sz","Pr"/) Open (3, file=trim(Fname), Status="Unknown") write (3, "(A)") " "//Fname Do ist=1, Nst write (3,"(A40)") trim(Station(ist)) write (3,"(8A8)") " ",(Title(j), j=1,Ntq) Do j=1,Ntq write (3,"(A4,4x,7F8.4)") Title(j),(X(ist,j,k),k=1,Ntq) Enddo Enddo Close (3) End Subroutine Write_HSTQ !------------------------------------------------------------------------------ Subroutine Calc_Write_HSTQ1 Integer :: idx, N, it, ic, j, k Real :: TBj, TBk, Sj, Sk, Rjk, tmp Character (Len=3) :: IDXs(Nidx)=(/"mop", "mow", "moz", "msp", "msw", "msz", "w1p", & "w1w", "w1z", "w2p", "w2w", "w2z", "mpr", "wpr"/) N = NYear*NWeek Do idx=1, Nidx Do j=1, Nst Do k=1, Nst TBj=0; TBk=0; Sj=0; Sk=0; Rjk=0; ic=0 Do it=1, N if (All_IDX(j,it,idx) /= -99. .AND. All_IDX(k,it,idx) /= -99.) then ic=ic+1 TBj=TBj+All_IDX(j,it,idx); TBk=TBk+All_IDX(k,it,idx); Sj=Sj+All_IDX(j,it,idx)*All_IDX(j,it,idx) Sk=Sk+All_IDX(k,it,idx)*All_IDX(k,it,idx); Rjk=Rjk + All_IDX(j,it,idx)*All_IDX(k,it,idx) endif Enddo tmp = ic; TBj=TBj/tmp; TBk=TBk/tmp; Sj=Sj/tmp-TBj*TBj; Sk=Sk/tmp-TBk*TBk; Rjk=Rjk/tmp TQ_Stn(idx,j,k) = (Rjk-TBj*TBk)/SQRT(Sj*Sk) Enddo Enddo Enddo Open (3, file="HSTQ_Stn.txt", Status="Unknown") write (3, "(A)") " He so Tuong quan cua cac Chi so giua cac tram" Do idx=1, Nidx write (3,"(A40)") "Index="//IDXs(idx) write (3,"(100(1x,A7)") " ",(trim(Station(j)),j=1, Nst) Do j=1,Nst write (3,"(A7,1x,100F8.4)") trim(Station(j)),(TQ_Stn(idx,j,k),k=1,Nst) Enddo Enddo Close (3) End Subroutine Calc_Write_HSTQ1 !------------------------------------------------------------------------------ Subroutine HSTQ (X, Y, ist) Real, Dimension (Nst_Max,NYear*NWeek,Ntq) :: X ! Luu tat ca cac chuoi PDSI de tinh HSTQ Real, Dimension (Nst_Max,Ntq,Ntq) :: Y Integer :: ist, N, it, ic, j, k Real :: TBj, TBk, Sj, Sk, Rjk, tmp N = NYear*NWeek Do j=1, Ntq Do k=1, Ntq TBj=0; TBk=0; Sj=0; Sk=0; Rjk=0; ic=0 Do it=1, N if (X(ist,it,j) /= -99. .AND. X(ist,it,k) /= -99.) then ic=ic+1 TBj=TBj+X(ist,it,j); TBk=TBk+X(ist,it,k); Sj=Sj+X(ist,it,j)*X(ist,it,j) Sk=Sk+X(ist,it,k)*X(ist,it,k); Rjk=Rjk + X(ist,it,j)*X(ist,it,k) endif Enddo tmp = ic TBj=TBj/tmp; TBk=TBk/tmp; Sj=Sj/tmp-TBj*TBj; Sk=Sk/tmp-TBk*TBk; Rjk=Rjk/tmp Y(ist,j,k) = (Rjk-TBj*TBk)/SQRT(Sj*Sk) Enddo Enddo End Subroutine HSTQ !------------------------------------------------------------------------------ Subroutine Write_Table (Fname, X, Nst, Nmon_week) Integer :: Nst, Nmon_week, Ntime, iyr, imw Character (Len=*) :: Fname Real :: X(Nst_Max,NYear*NWeek), Tmp Open (3, File=trim(Fname), Status="Unknown") write (Fmt, "(A,I3.3,A)") "3x,A3,1x,", Nst+1,"(1x,A7)" write (3, Fmt) "M/W", (Station(ist), ist=1, Nst), "Aver" write (Fmt,"(A,I3.3,A)") "I4.4,I2.2,", Nst+1,"F8.2" Ntime=0 Do iyr=YBas1, YBas2 Do imw=1, Nmon_week Ntime=Ntime+1; Tmp = 0. Do ist=1, Nst Tmp = Tmp + X(ist,Ntime) Enddo Tmp = Tmp/Real(Nst) write (3, Fmt) iyr,imw,(X(ist,Ntime), ist=1,Nst), Tmp Enddo Enddo Close (3) End Subroutine Write_Table !------------------------------------------------------------------------------ Subroutine Read_Table (Fname, Nmon_week, Outp) Integer :: Nmon_week, Ntime, iyr, imw, year Real :: wk(Nmon_week), Outp(Nmon_week*NYear) Character (Len=*) :: Fname Open (1, File=trim(Fname), status="Old") Ntime = 0 Do iyr=YBas1, YBas2 read (1, *) year, (wk(imw), imw=1, Nmon_week) Do imw=1, Nmon_week Ntime = Ntime+1; Outp(Ntime) = wk(imw) Enddo Enddo End Subroutine Read_Table !------------------------------------------------------------------------------ Subroutine Write_to_PDSI (Dat_Dir, DatBL, Nmon_week, ist) Integer :: ivar, ist, iyr, imw, isys, system Real :: Tmp Character (Len=*) :: Dat_Dir Character (Len=255) :: Fname Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NWeek) :: DatBL !Do ist=1, Nst Fname = "/home/tanpv/PDSI/PDSI/my_dat/" !//trim(Dat_Dir)//"/"//trim(Station(ist)) isys = system ('mkdir -p '//trim(Fname)) Open (1, file=trim(Fname)//"/"//"parameter", status="Unknown") write (1,*) AWC(ist), xlat(ist) Do ivar=1, NVar if (ivar == 1) then Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_T", status="Unknown") Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_T_normal",status="Unknown") else Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_P", status="Unknown") Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_P_normal",status="Unknown") endif Do imw=1, Nmon_week ! Tinh trung binh khi hau Tmp = 0 Do iyr=YBas1, YBas2 Tmp = Tmp + DatBL(ivar,ist,iyr,imw) Enddo DatNor(ivar,ist,imw) = Tmp/Real(NYear) Enddo write (4,"(52F8.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week) Do iyr=YBas1, YBas2 write (3,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week) Enddo Close (3); Close (4) Enddo Close (1) !Enddo End Subroutine Write_to_PDSI !------------------------------------------------------------------------------ Subroutine Read_data_Baseline(Fname,DatBL,Nmon_week) ! Data from Tan Integer :: IO, iyt, ist, iyr, year, imw, N Character (Len=*) :: FName Character (Len=21) :: Stname Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NWeek) :: DatBL !SSMax = 1. ; AWC=6.7 ; SUMax = AWC-SSMax ! Giang Prog. Do iyt=1, NVar ! =1: T2m; =2: R Open (1,file=trim(BL_Dir)//trim(YT(iyt))//trim(Fname)) ist = 1; IO = 0 Do !ist=1, Nst read (1, IOStat=IO, "(A)") Stname if (IO /= 0) exit ! print*,"--->",Stname Station(ist) = trim(Stname) read (1,IOStat=IO,"(2F8.3,F8.1)") xlon(ist), xlat(ist), HH(ist) if (IO /= 0) exit Do iyr=YBas1, YBas2 read (1,*) year, (DatBL(iyt,ist,iyr,imw), imw=1,Nmon_week) Do imw=1, Nmon_week ! Conver to Inch OR Fahrenheit if (iyt == 1) then DatBL(iyt,ist,iyr,imw) = DatBL(iyt,ist,iyr,imw)*1.8 + 32. else DatBL(iyt,ist,iyr,imw) = DatBL(iyt,ist,iyr,imw)/25.4 endif Enddo Enddo ! print*,ist," ",Station(ist),xlon(ist), xlat(ist), HH(ist) ist = ist + 1 Enddo Nst = ist - 1 Enddo End Subroutine Read_data_Baseline !------------------------------------------------------------------------------ SUBROUTINE putAWC(dat,lon,lat,nst) integer :: nst, ist real :: dat(nst), lon(nst), lat(nst) integer, parameter :: ncols = 1099, nrows = 549 integer :: icol, irow, itmp character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc" real, parameter :: xcor = -180.0, ycor = -89.9056530001 real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99 real :: indat(ncols,nrows), x(ncols), y(nrows) ! <-~~~=|==-----------------------------------------------------------