Program FA_Monthly Integer, Parameter :: NYr1=1980, NYr2=2014, Mmon=12, Nmon = 12, Nsea = 4 ! 4 mua (12-2,3-5,6-8,9-11) Integer, Parameter :: NtrMx = 16 ! So tram cuc dai co SL (So chi so IDX) Character (Len=50), Parameter :: Fname_idx(NtrMx) = & (/"DMI.txt","BEST.txt","MEI.txt","Nino3.4.txt","ONI.txt","SOI.txt","WP.txt","WPC850.txt", & "QBO.txt","RINDO_SLPA.txt","CPAC850.txt","PDO.txt","TNI.txt","PNA.txt","REQSOI.txt","WHWP.txt"/) Include 'FTNT.INC' Character *11 :: StaName (NtrMx) Character (Len = 80) :: FInp, FOut, FOut1, Grads_File, Out_Dir, Out_Var Character (Len = 30) :: Fmt, StartTime Character (Len = 11) :: TenTram(NtrMx), VarName(NtrMx), Sample(Nmax) Real :: Lon(NtrMx), Lat(NtrMx), H(NtrMx) Integer :: NTram, IO, yr, Ntr, iyr, ist, imon Real*8 :: X(Nmax, NtrMx), W1(Nmax) Real*8 :: XX(Nmax, Mmax), X1(Nmax, Mmax, Mmon), XX_R(Nmax,Mmax), Xsea(Nmax, Mmax, Nsea) Real*8 :: XMon(Mmon-1,NtrMx) ! Trung binh thang nhieu nam cua tat ca cac tram Integer :: YYDD(Nmax), YY(Nmax) Real :: Con Logical :: NT1_TT2 Character (Len=4) :: Var ! Ten bien R=11, T2m=14 Common /StationInfo/ TenTram, VarName, Sample Real :: Gr_ENS(Nmax, Mmax) ! Mang chua SL cho Grads (ploting ensemble) ! Nmax: Ens ; Mmax: Time !---------------------------------------------------------------------------- TenTram(1)="DMI"; TenTram(2)="BEST"; TenTram(3)="MEI"; TenTram(4)="Nino3.4"; TenTram(5)="ONI"; TenTram(6)="SOI"; TenTram(7)="WP"; TenTram(8)="WPC850" TenTram(9)="QBO"; TenTram(10)="RINDO_SLPA"; TenTram(11)="CPAC850"; TenTram(12)="PDO"; TenTram(13)="TNI"; TenTram(14)="PNA"; TenTram(15)="REQSOI"; TenTram(16)="WHWP" !---------------------------------------------------------------------------- Nyear = Nyr2-Nyr1+1; Ntr = NtrMx; Out_Dir="Out/" Out_Var="_" Call Read_Clim_IDX print*, "Call Read_Clim_IDX ... Done!" !Call Monthly_Variation !print*, "Call Monthly_Variation... Done!" Call All_Timeseries print*, "All_Timeseries... Done!" !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------------------! Subroutine Read_Clim_IDX Integer :: iidx, iyr, imon, year, IO, it Character (Len=50) :: Fname Real :: X(Nmon) Logical :: Exists !Allocate(ClimIDX(Nidx, Nyear, Nmon), Dat_IDX(Ntime, Nidx)) ! X1(Nmax, Mmax, Mmon) Do iidx=1, Ntr Fname="../../Nafosted/Data/ClimIDX/"//trim(Fname_idx(iidx)) Inquire(file=trim(Fname), Exist=Exists) if (.Not.Exists) then print*,"File ",trim(Fname)," does not exist" cycle endif Open(1, file=trim(Fname), Status="Old") read(1,*) ! Bo dong dau Do read(1,*,IOStat=IO) year, (X(imon),imon=1,Nmon) if (IO /= 0) exit Do iyr=NYr1, NYr2 if (year==iyr) then; YY(iyr-NYr1+1) = iyr; X1(iyr-NYr1+1,iidx,1:Nmon)=X(1:Nmon); endif Enddo Enddo Close (1) Enddo End Subroutine Read_Clim_IDX !----------------------------------------------------------------------- Subroutine Monthly_Variation Do k=1,Nmon ! Mmon=1,2,...,12 correspond to Jan, Feb,..., Dec Print*,"Tinh cho thang: ",k write (Fmt,'(I2.2)') k ! Fout = trim(Out_Dir)//trim(Out_Var)//"PTNT_"//trim(Fmt)//".txt" ; Open (3,file=Fout,status="Unknown") write (Grads_File,'(A,I2.2)') trim(Out_Var)//'PTNT_',k CON = 1. !----------------- Phuong an lay mang luoi tram la cac bien ------------ VarName(1:Ntr) = TenTram(1:Ntr) Do i=1,Nyear write(Sample(i),'(I4)') YY(i) XX(i,1:Ntr) = X1(i,1:Ntr,k) Enddo NT1_TT2 = .true.; write (StartTime,"(A,I4,A)") "jan", NYr1, " 1yrs" CALL PTNT(XX, Nyear, Ntr, CON, VarName, Sample, Grads_File, Out_Dir, NT1_TT2, StartTime) ! Var la cac tram Close(3) Enddo ! k=1,Nmon End Subroutine Monthly_Variation !----------------------------------------------------------------------- Subroutine All_Timeseries Integer :: it ! Do k=1,Nmon ! Mmon=1,2,...,12 correspond to Jan, Feb,..., Dec Print*,"Tinh cho toan chuoi " write (Fmt,'(A)') "All" ! Fout = trim(Out_Dir)//trim(Out_Var)//"PTNT_"//trim(Fmt)//".txt" ; Open (3,file=Fout,status="Unknown") write (Grads_File,'(A,I2.2)') trim(Out_Var)//'PTNT_All' CON = 1. !----------------- Phuong an lay mang luoi tram la cac bien ------------ VarName(1:Ntr) = TenTram(1:Ntr); it=0 Do i=1,Nyear Do imon=1, Nmon it=it+1; write(Sample(it),'(I4.4,I2.2)') YY(i), imon XX(it,1:Ntr) = X1(i,1:Ntr,imon) Enddo Enddo NT1_TT2 = .true.; write (StartTime,"(A,I4,A)") "jan", NYr1, " 1mo" CALL PTNT(XX, Nyear*Nmon, Ntr, CON, VarName, Sample, Grads_File, Out_Dir, NT1_TT2, StartTime) ! Var la cac tram NT1_TT2 = .false.; write (Grads_File,'(A)') trim(Out_Var)//'PTNT_All_PCA' CALL PTNT(XX, Nyear*Nmon, Ntr, CON, VarName, Sample, Grads_File, Out_Dir, NT1_TT2, StartTime) ! Var la cac tram Close(3) ! Enddo ! k=1,Nmon End Subroutine All_Timeseries !----------------------------------------------------------------------- END ! Main Prog !----------------------------------------------------------------------- Subroutine Write_GrADS (Ntr, Nyear, K, Grads_File, Out_Dir, NT1_TT2, StartTime) Integer, Parameter :: NYr1=1980, NYr2=2014, Mmon=12 Integer, Parameter :: NtrMx = 20 ! So tram cuc dai co SL INCLUDE 'FTNT.INC' real*8 :: TaiTrong(MMax,MMax) ! MOI COT LA TAI TRONG CUA MOT NHAN TO REAL*8 :: NhanTo(NMax,MMax) ! MOI COT LA MOT NHAN TO REAL*8 :: TT_DAQUAY(MMAX,MMAX), NT_DAQUAY(NMAX,MMAX) Real, Allocatable :: Timeseries1(:,:), Timeseries2(:,:), Load1(:,:), Load2(:,:) Character (Len = 11) :: TenTram(NtrMx), VarName(NtrMx), Sample(Nmax) Real :: Lon(NtrMx), Lat(NtrMx), H(NtrMx), Var Integer :: ift, ist, iyr ! Bien chay theo Factor, Station, Year Integer :: Ntime, Nsta Character (Len=8) :: STID Character *(*) :: Grads_File, Out_Dir, StartTime Logical :: NT1_TT2 Common /StationInfo/ TenTram, VarName, Sample Common /Factors/ TaiTrong, NhanTo, TT_DAQUAY, NT_DAQUAY Real :: RLAT, RLON, TIM, tmp Integer :: NLEV,NFLAG, irec print*,"Ntr,Nyear=",Ntr,Nyear !---------------------------------------------------------------------------------- if (NT1_TT2) then ! Factors are Timeseries Allocate (Timeseries1(NMax,MMax), Timeseries2(NMax,MMax)) Allocate (Load1(MMax,MMax), Load2(MMax,MMax)) Timeseries1(:,:) = NhanTo(:,:); Timeseries2(:,:) = NT_DAQUAY(:,:) Load1(:,:) = TaiTrong(:,:); Load2(:,:) = TT_DAQUAY(:,:) Ntime = Nyear; NSta = Ntr else ! Loadings are Timeseries Allocate (Timeseries1(MMax,MMax), Timeseries2(MMax,MMax)) Allocate (Load1(NMax,MMax), Load2(NMax,MMax)) Timeseries1(:,:) = TaiTrong(:,:); Timeseries2(:,:) = TT_DAQUAY(:,:) Load1(:,:) = NhanTo(:,:); Load2(:,:) = NT_DAQUAY(:,:) Ntime = Ntr; NSta = Nyear endif ! Plotting time series Open (33,file=trim(Out_Dir)//trim(Grads_File)//'_timeseries.dat',FORM='Binary',access="direct",recl=4) irec=0 Do iyr=1,Ntime !Nyear Do ift=1,K irec=irec+1 ; tmp = Timeseries1(iyr,ift) ; write (33,rec=irec) tmp Enddo Do ift=1,K irec=irec+1 ; tmp = Timeseries2(iyr,ift) ; write (33,rec=irec) tmp Enddo Enddo Close (33) Open (4, file=trim(Out_Dir)//trim(Grads_File)//"_timeseries.ctl") write (4,"(A)") "dset ^"//trim(Grads_File)//"_timeseries.dat" write (4,"(A)") "title Chuoi thoi gian cac Nhan to" write (4,"(A)") "undef -99." write (4,"(A)") "xdef 1 linear 21 0.5" write (4,"(A)") "ydef 1 linear 105 0.5" write (4,"(A)") "zdef 1 linear 1000 1" write (4,"(A,I4,A)") "tdef ",Ntime," linear "//trim(StartTime) !00z01Jan1970 1yrs" write (4,"(A,I3)") "vars ",2*K Do ift=1,K write (4,"(A,I2.2,A,A)") "f",ift," 0 99 "," Nhan to chua quay" Enddo Do ift=1,K write (4,"(A,I2.2,A,A)") "fq",ift," 0 99 "," Nhan to da quay" Enddo write (4,"(A)") "endvars" Close(4) print*,"Calling Write_GrADS Done !!!!" !---------------------------------------------------------------------------------- End Subroutine Write_GrADS !----------------------------------------------------------------------- ! CHUONG TRINH GOI TINH TAI TRONG NHAN TO VA CAC NHAN TO SUBROUTINE PTNT(XX,N,M,CON,VarName, Sample,Grads_File, Out_Dir, NT1_TT2, StartTime) ! INCLUDE 'FTNT.INC' real*8 :: B(MMax),D(MMax),S(MMax),XBAR(MMax) real*8 :: T(MMax),V(MMax*MMax),R(MMax*MMax),TV(MMax*MMax) real*8 :: xx(NMax,MMax), x(NMax*MMax), Y(NMax,MMax) real*8 :: TriRieng(MMax),VTORieng(MMax,MMax) ! MOI COT LA MOT VECTO real*8 :: TaiTrong(MMax,MMax) ! MOI COT LA TAI TRONG CUA MOT NHAN TO REAL*8 :: NhanTo(NMax,MMax) ! MOI COT LA MOT NHAN TO ! REAL*8 :: TT_DAQUAY(MMAX,MMAX), NT_DAQUAY(NMAX,MMAX) REAL*8 :: TMP,tmp1 Character*20 :: Fmt real :: Con, Threshold Character (Len = 11) :: VarName(MMax), Sample(MMax) Integer :: M, N, I, J, K, IO, Nthre Character *(*) :: Grads_File, Out_Dir, StartTime Common /Factors/ TaiTrong, NhanTo, TT_DAQUAY, NT_DAQUAY Threshold = 0.995 ! Giu lai 95% luong thong tin WRITE(*,*)' BAT DAU C.T FACTOR' !------------------------------------------------------------------------ k=0 do j=1,m do i=1,n k=k+1 x(k)=xx(i,j) enddo enddo write (3, "(A)") "PHAN TICH NHAN TO" WRITE(3,'(A17,I3)')'Dung luong mau = ',N WRITE(3,'(A17,I3)')'So bien = ',M IO=1 CALL CORRE(N,M,IO,x,xbar,s,V,R,d,b,T) ! Tinh Ma tran Do lech (Y) cua X Tmp=SQRT(REAL(N)) do i=1,N do j=1,M Y(i,j) = (XX(i,j)-XBAR(j))/S(J)/TMP enddo enddo write (Fmt, '(A,I3.3,A)') '(',M,'(1x,A7))' write (3, Fmt) (VarName(j),j=1,M) write (3, '(A)') 'Trung binh ' write (Fmt,'(A1,I3.3,A5)') "(", M, "F8.1)" WRITE(3,Fmt)(XBAR(J),J=1,M) write (3, '(A)') 'Do lech chuan ' WRITE (3,Fmt)(S(J),J=1,M) WRITE(3,'(A)') 'Ma tran tuong quan ' write (Fmt,'(A6,I3.3,A5)') "(I3,", M, "F8.4)" Do I=1,M DO J=1,M if ((I-J) < 0) then L=I+(J*J-J)/2 else L=J+(I*I-I)/2 endif D(J)=R(L) ENDDO WRITE(3,Fmt)I,(D(J),J=1,M) ENDDO !--------------------------------------------- MV=0 CALL EIGEN(R,V,M,MV) ! Luu cac tri rieng va Vecto rieng tmp = 0. ! Tong cac tri rieng do I=1,M L=I+(I*I-I)/2 TriRieng(I) = R(L) tmp = tmp+TriRieng(I) ! Tong cac tri rieng enddo tmp1 = 0; Nthre = 0 ! De so voi Tong cac tri rieng Do i=1,M tmp1 = tmp1+TriRieng(I); Nthre=Nthre+1 if (tmp1/tmp > Threshold .OR. Nthre > 3) exit ! Vuot nguong nhung khong qua 4 nhan to Enddo CON = TriRieng(Nthre) ! Xac dinh lai nguong gia tri cua Tri rieng duoc giu lai L=0 do J=1,M do I=1,M L=L+1 VTORieng(i,j) = V(L) enddo enddo WRITE(*,*)' TRI RIENG 1,2,3,4,5=',R(1),R(3),R(6),R(10),R(15) write(*,'(A,F10.3)') "Gia tri nguong CON: ", CON write(*,'(A,F10.3)') "So phan tram luong thong tin duoc giu lai: ",tmp1/tmp*100. ! WRITE(*,'(A\)')' VAO LAI THAM SO (CON)=' ! READ(*,*)CON CALL TRACE(M,R,CON,K,D) DO I=1,K L=I+(I*I-I)/2 S(I)=R(L) ENDDO WRITE(3,'(A)') 'Cac gia tri rieng ' ! Mang R(..) write (Fmt,'(A,I3.3,A5)') "(", K, "F8.2)" WRITE(3,Fmt)(S(J),J=1,K) ! Mang S ~ R WRITE(3,'(A)') 'Ty so tich luy cua cac tri rieng ' WRITE(3,Fmt)(D(J),J=1,K) WRITE(3,'(A)') 'Cac vec to rieng' write (Fmt,'(A,I3.3,A5)') "(A,I3,", M, "F8.3)" L=0 DO J=1,K DO I=1,M L=L+1 D(I)=V(L) ENDDO WRITE(3,Fmt) "Vector ",j, (D(I),I=1,M) ENDDO CALL LOAD(M,K,R,V) ! Vecto Tai trong = Vecto rieng * SQRT(Gia tri rieng) WRITE(3,'(A,I2,A)') 'Ma tran cac Tai trong nhan to (',K,' Nhan to )' write (Fmt,'(A,I2.2,A5)') "(I4,1x,A11,1x", K, "F8.4)" DO I=1,M DO J=1,K L=M*(J-1)+I D(J)=V(L) ENDDO ! Luu Ma tran Tai trong Nhan to (K nhan to) do j=1,K TaiTrong(I,J) = D(j) enddo WRITE(3,Fmt) i, VarName(i), (D(J),J=1,K) ENDDO If (K > 1) Then CALL VARMX(M,K,V,NC,TV,B,T,D) NV=NC+1 ! In Ket qua kiem tra (Unit=4) ! WRITE(4,'(2A10)') 'Vong lap ','Phuong sai' ! DO I=1,NV ! NC=I-1 ! WRITE(4,'(I4,6X,F20.10)') NC, TV(I) ! ENDDO WRITE(3,'(A,I2,A)') 'Ma tran cac Tai trong nhan to da quay (',K,' Nhan to)' write (Fmt,'(A,I2.2,A5)') "(I4,1x,A11,1x", K, "F8.4)" DO I=1,M DO J=1,K L=M*(J-1)+I TT_DAQUAY(I,J)=V(L) S(J)=V(L) ENDDO WRITE(3,Fmt) i, VarName(i), (S(J),J=1,K) ENDDO ! Thu so sanh ! WRITE(3,'(A)') 'Ma tran cac nhan to da quay (de so sanh)' ! Do i=1,M ! WRITE(3,Fmt) i, VarName(i), (TT_DAQUAY(I,J),j=1,K) ! Enddo ! Tinh Ma tran Nhan to (F) & NTO DA QUAY do i=1,N do j=1,K NhanTo(i,j)=0.0 NT_DAQUAY(I,J)=0.0 do L=1,M NhanTo(i,j)=NhanTo(i,j)+Y(i,L)*TaiTrong(L,j) NT_DAQUAY(I,J)=NT_DAQUAY(I,J)+Y(I,L)*TT_DAQUAY(L,J) enddo NhanTo(i,j)=NhanTo(i,j)/(TriRieng(j)) NT_DAQUAY(i,j)=NT_DAQUAY(i,j)/(TriRieng(j)) enddo enddo ! In Ket qua ! write(3,*)' Tai trong Nhan to' ! write (Fmt,'(A,I2.2,A5)') "(I2,1x,A11,1x", K, "F8.4)" ! do i=1,M ! write(3,Fmt) i, VarName(i), (TaiTrong(i,j),j=1,K) ! enddo write(3,*)' Ma tran Nhan to' do i=1,N ! write(3,Fmt) i, Sample(i), (NhanTo(i,j),j=1,K) write(3,"(I4,1x,A11,1x,10F8.4)") i, Sample(i), (NhanTo(i,j),j=1,K) enddo write(3,*)' Ma tran Nhan to da quay' do i=1,N ! write(3,Fmt) i, Sample(i), (NT_DaQuay(i,j),j=1,K) write(3,"(I4,1x,A11,1x,10F8.4)") i, Sample(i), (NT_DaQuay(i,j),j=1,K) enddo ! Write out for GrADS print*,"Calling Write_GrADS ..." Call Write_GrADS (M, N, K, Grads_File, Out_Dir, NT1_TT2, StartTime) !-------------------- Else WRITE(3,'(A,I3,A)') 'Chi co',K,' nhan to duoc giu lai. Khong thuc hien phep quay' EndIF END SUBROUTINE PTNT !============================================================== ! CHUONG TRINH TIM K TRI RIENG SUBROUTINE TRACE(M,R,CON,K,D) INCLUDE 'FTNT.INC' REAL*8 R(MMAX*MMAX), D(MMAX) Integer :: M, K, I, L ! WRITE(*,*)' BAT DAU C.T TRACE' FM=M L=0 DO I=1,M L=L+I D(I)=R(L) ENDDO K=0 DO I=1,M IF (D(I) < CON) Exit K=K+1 D(I)=D(I)/FM ENDDO DO I=2,K D(I)=D(I)+D(I-1) ENDDO END SUBROUTINE TRACE ! !----------------------------------------------------------------------- ! CHUONG TRINH TINH TAI TRONG NHAN TO SUBROUTINE LOAD(M,K,R,V) INCLUDE 'FTNT.INC' REAL*8 R(MMAX*MMAX),V(MMAX*MMAX),SQ Integer :: M, K, I, J, L ! WRITE(*,*)' BAT DAU C.T LOAD' L=0 JJ=0 DO J=1,K JJ=JJ+J SQ=SQRT(R(JJ)) DO I=1,M L=L+1 V(L)=SQ*V(L) ENDDO ENDDO END SUBROUTINE LOAD !----------------------------------------------------------------------- ! CHUONG TRINH TINH TRI RIENG VA VEC TO RIENG ! SUBROUTINE EIGEN(A,R,N,MV) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),R(MMAX*MMAX) REAL*8 ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2,SINCS,RANGE Integer :: N, I, J, M ! WRITE(*,*)' BAT DAU EIGEN' 5 RANGE=1.E-6 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0. IF(I-J) 20,15,20 15 R(IJ)=1. 20 CONTINUE 25 ANORM=0. DO 35 I=1,N DO 35 J=1,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR)130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GOTO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GOTO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-N) 135,140,135 135 M=M+1 GOTO 60 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GOTO 55 150 IF(IND-1) 160,155,160 155 IND=0 GOTO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE WRITE(*,*)' KET THUC EIGEN' RETURN END SUBROUTINE EIGEN ! !----------------------------------------------------------------------- ! CHUONG TRINH QUAY NHAN TO SUBROUTINE VARMX(M,K,A,NC,TV,H,F,D) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),TV(NMAX),H(MMAX),F(MMAX),D(MMAX) Integer :: M, K, I, J WRITE(*,*)' BAT DAU C.T VARMX' EPS=0.0000116 TVLT=0. LL=K-1 NV=1 NC=0 FN=M FFN=FN*FN CONS=.7071066 DO 110 I=1,M H(I)=0 DO 110 J=1,K L=M*(J-1)+I 110 H(I)=H(I)+A(L)*A(L) DO 120 I=1,M 115 H(I)=SQRT(H(I)) DO 120 J=1,K L=M*(J-1)+I 120 A(L)=A(L)/H(I) GOTO 132 130 NV=NV+1 TVLT=TV(NV-1) 132 TV(NV)=0. DO 150 J=1,K AA=0. BB=0. LB=M*(J-1) DO 140 I=1,M L=LB+I CC=A(L)*A(L) AA=AA+CC 140 BB=BB+CC*CC 150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN IF(NV-51) 160 ,430,430 160 IF((TV(NV)-TVLT)-(1.E-7)) 170,170,190 170 NC=NC+1 IF(NC-3)190,190 ,430 190 DO 420 J=1,LL L1=M*(J-1) II=J+1 DO 420 K1=II,K L2=M*(K1-1) AA=0 BB=0 CC=0 DD=0 DO 230 I=1,M L3=L1+I L4=L2+I U=(A(L3)+A(L4))*(A(L3)-A(L4)) T=A(L3)*A(L4) T=T+T CC=CC+(U+T)*(U-T) DD=DD+2*U*T AA=AA+U 230 BB=BB+T T=DD-2*AA*BB/FN B=CC-(AA*AA-BB*BB)/FN IF(T-B) 280,240,320 240 IF((T+B)-EPS)420,250,250 250 COS4T=CONS SIN4T=CONS GOTO 350 280 TAN4T=ABS(T)/ABS(B) IF(TAN4T-EPS)300,290,290 290 COS4T=1/SQRT(1+TAN4T*TAN4T) SIN4T=TAN4T*COS4T GOTO 350 300 IF(B) 310,420,420 310 SINP=CONS COSP=CONS GOTO 400 320 CTN4T=ABS(T/B) IF(CTN4T-EPS) 340,330,330 330 SIN4T=1/SQRT(1+CTN4T*CTN4T) COS4T=CTN4T*SIN4T GOTO 350 340 COS4T=0. SIN4T=1. 350 COS2T=SQRT((1+COS4T)/2) SIN2T=SIN4T/(2*COS2T) 355 COST=SQRT((1+COS2T)/2) SINT=SIN2T/(2*COST) IF(B) 370,370,360 360 COSP=COST SINP=SINT GOTO 380 370 COSP=CONS*COST+CONS*SINT 375 SINP=ABS(CONS*COST-CONS*SINT) 380 IF(T) 390,390,400 390 SINP=-SINP 400 DO 410 I=1,M L3=L1+I L4=L2+I AA=A(L3)*COSP+A(L4)*SINP A(L4)=-A(L3)*SINP+A(L4)*COSP 410 A(L3)=AA 420 CONTINUE GOTO 130 430 DO 440 I=1,M DO 440 J=1,K L=M*(J-1)+I 440 A(L)=A(L)*H(I) NC=NV-1 DO 450 I=1,M 450 H(I)=H(I)*H(I) DO 470 I=1,M F(I)=0. DO 460 J=1,K L=M*(J-1)+I 460 F(I)=F(I)+A(L)*A(L) 470 D(I)=H(I)-F(I) RETURN END SUBROUTINE VARMX ! !----------------------------------------------------------------------- ! SUBROUTINE CORRE(N,M,I0,X,XBAR,STD,RX,R,B,D,T) INCLUDE 'FTNT.INC' REAL*8 X(NMAX*MMAX),XBAR(MMAX),STD(MMAX),RX(MMAX*MMAX) REAL*8 R(MMAX*MMAX),B(MMAX),D(MMAX),T(MMAX) Integer :: N, M, L, J, K, IO, I IO = 1 DO 100 J=1,M B(J)=0. 100 T(J)=0. K=(M*M+M)/2 DO 102 I=1,K 102 R(I)=0. FN=N L=0 IF(I0)105,127,105 105 DO 108 J=1,M DO 107 I=1,N L=L+1 107 T(J)=T(J)+X(L) XBAR(J)=T(J) 108 T(J)=T(J)/FN DO 115 I=1,N JK=0 L=I-N DO 110 J=1,M L=L+N D(J)=X(L)-T(J) 110 B(J)=B(J)+D(J) DO 115 J=1,M DO 115 K=1,J JK=JK+1 115 R(JK)=R(JK)+D(J)*D(K) GOTO 205 127 IF(N-M)130,130,135 130 KK=N GOTO 137 135 KK=M 137 DO 140 I=1,KK ! CALL DATA(M,D) DO 140 J=1,M T(J)=T(J)+D(J) L=L+1 140 RX(L)=D(J) FRR=KK DO 150 J=1,M XBAR(J)=T(J) 150 T(J)=T(J)/FRR L=0 DO 180 I=1,KK JK=0 DO 170 J=1,M L=L+1 170 D(J)=RX(L)-T(J) DO 180 J=1,M B(J)=B(J)+D(J) DO 180 K=1,J JK=JK+1 180 R(JK)=R(JK)+D(J)*D(K) IF(N-KK) 205,205,185 185 KK=N-KK DO 200 I=1,KK JK=0 ! CALL DATA(M,D) DO 190 J=1,M XBAR(J)=XBAR(J)+D(J) D(J)=D(J)-T(J) 190 B(J)=B(J)+D(J) DO 200 J=1,M DO 200 K=1,M JK=JK+1 200 R(JK)=R(JK)+D(J)*D(K) 205 JK=0 DO 210 J=1,M XBAR(J)=XBAR(J)/FN DO 210 K=1,J JK=JK+1 210 R(JK)=R(JK)-B(J)*B(K)/FN JK=0 DO 220 J=1,M JK=JK+J 220 STD(J)=SQRT(ABS(R(JK))) DO 230 J=1,M DO 230 K=J,M JK=J+(K*K-K)/2 L=M*(J-1)+K RX(L)=R(JK) L=M*(K-1)+J RX(L)=R(JK) IF(STD(J)-STD(K)) 225,222,225 222 R(JK)=1.0 GOTO 230 225 STD2=STD(J)*STD(K) ! PRINT*,STD2 R(JK)=R(JK)/STD2 ! 225 R(JK)=R(JK)/(STD(J)*STD(K)) ! WRITE(3,500) R(JK) ! 500 FORMAT(5X,F10.4) 230 CONTINUE ! WRITE(*,*)' Dung luong mau=',fn FN=SQRT(FN-1.) DO 240 J=1,M 240 STD(J)=STD(J)/FN L=-M DO 250 I=1,M L=L+M+1 250 B(I)=RX(L) ! WRITE(*,*)' KET THUC CORRE' RETURN END SUBROUTINE CORRE !----------------------------------------------------------------------- ! SUBROUTINE DATA(M,D) RETURN END SUBROUTINE DATA !----------------------------------------------------------------------- ! SUBROUTINE FACTOR (X,F,A,N,M,K) ! !...TINH F(N,K) TU PHUOG TRINH X(N,M) = F(N,K)*AT(K,M) !...VOI AT LA CHUYEN VI CUA A(M,K) ! INCLUDE 'FTNT.INC' INTEGER N,M,K, J, I REAL*8 X(NMAX,MMAX), F(NMAX,MMAX), A(MMAX,MMAX), XA(NMAX,MMAX) REAL*8 AT(MMAX,MMAX), A1(MMAX,MMAX), AA(MMAX*MMAX),A2(MMAX,MMAX) !...CHUYEN VI DO I = 1,M DO J = 1,K AT(J,I) = A(I,J) ENDDO ENDDO !...TICH AT(K,M)*A(M,K) = A1(K,K) DO I = 1,K DO J = 1,K A1(I,J) = 0.0 DO L = 1,M A1(I,J) = A1(I,J) + AT(I,L)*A(L,J) ENDDO ENDDO ENDDO !...CHUYEN THANH MANG 1 CHIEU DE TINH NGHICH DAO L = 0 DO J = 1,K DO I = 1,K L = L + 1 AA(L) = A1(I,J) ENDDO ENDDO !...NGHICH DAO MA TRAN AA CALL MINV(AA,K,D) !...LUU AA MOT CHIEU VAO A1 HAI CHIEU L = 0 DO J = 1,K DO I = 1,K L = L + 1 A1(I,J) = AA(L) ENDDO ENDDO !...TICH X(N,M)*A(M,K)*A1(K,K) = F(N,K) = XA(N,K)*A1(K,K) ! DO I = 1,N ! DO J = 1,K ! XA(I,J) = 0.0 ! DO L = 1,M ! XA(I,J) = XA(I,J) + X(I,L)*A(L,J) ! ENDDO ! ENDDO ! ENDDO ! ! DO I = 1,N ! DO J = 1,K ! F(I,J) = 0.0 ! DO L = 1,M ! F(I,J) = F(I,J) + XA(I,L)*A1(L,J) ! ENDDO ! ENDDO ! ENDDO !! !...TICH X(N,M)*A(M,K)*A1(K,K) = F(N,K) = X(N,K)*A2(K,K) DO I = 1,M DO J = 1,K A2(I,J) = 0.0 DO L = 1,K A2(I,J) = A2(I,J) + A(I,L)*A1(L,J) ENDDO ENDDO ENDDO ! DO I = 1,N DO J = 1,K F(I,J) = 0.0 DO L = 1,M F(I,J) = F(I,J) + X(I,L)*A2(L,J) ENDDO ENDDO ENDDO RETURN END SUBROUTINE FACTOR !----------------------------------------------------------------------- SUBROUTINE MINV(A,N,D) ! ,L,M) INCLUDE 'FTNT.INC' REAL*8 A(MMAX*MMAX),L(MMAX*MMAX),M(MMAX*MMAX) Integer :: N, k, I D=1. NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE J=L(K) IF(J-K)35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI)=HOLD 35 I=M(K) IF(I-K)45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) 40 A(JI)=HOLD 45 IF(BIGA)48,46,48 46 D=0. RETURN 48 DO 55 I=1,N IF(I-K)50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K)60,65,60 60 IF(J-K)62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K)70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE D=D*BIGA A(KK)=1./BIGA 80 CONTINUE K=N 100 K=(K-1) IF(K)150,150,105 105 I=L(K) IF(I-K)120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) 110 A(JI)=HOLD 120 J=M(K) IF(J-K)100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) 130 A(JI)=HOLD GOTO 100 150 RETURN END SUBROUTINE MINV !END