! Maximum Covariance Analysis of Drought Inde:wqx Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Nmon=12, Yr1=1980, Yr2=2020, Nidx=16 Integer, Parameter :: Nmode=6 !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40 Real, Parameter :: Slon1=70, Slon2=290, Slat1=-20, Slat2=40 Real, Allocatable :: Lon(:), Lat(:), xtime(:), Regs(:) Real, Allocatable :: Lon_sst(:), Lat_sst(:), Mask_sst(:,:) Real, Allocatable :: Mask(:,:), xlon(:), xlat(:), ylon(:), ylat(:) Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j Integer :: SeaPoint, Nsst, Nlon_sst, Nlat_sst Integer :: Yr_sta Real :: Miss_sst !---------------------------------------------------------------- Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ & "nino3", "nino34", "nino4", "oni", "censo", "rindo_slpa", & "soi", "reqsoi", "repac_slpa", & "cpac85", "whwp", "ammsst", "tna", "amonus", "pdo", "dmi"/) ! For MCA: Real, Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:) Real, Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:) Real, Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:) Real, Allocatable :: A(:,:), B(:,:) ! Timeseries of X & Y Real, Allocatable :: Uk(:,:), Vk(:,:), WWk(:,:), Wk(:) Real, Allocatable :: Ak(:,:), Bk(:,:), Corre(:), Probs(:) Real, Allocatable :: Accum(:), Percent(:) ! Accumulated percentage of Singular values Real, Allocatable :: XX(:), YY(:), Dr_IDX(:,:,:) ! (Nlon,Nlat,Ntime) Real, Allocatable :: U3D(:,:,:),V3D(:,:,:) ! Uk,Vk in 3D array Real, Allocatable :: Homo(:,:),Hete(:,:) ! Homo and Hetero maps (Nlon,Nlat,Nmode) Real, Allocatable :: p_Ho(:,:),p_He(:,:) ! Pval Integer :: Nt, Nx, Ny, isys, ilen Integer :: Mov_Ave_Step ! Buoc truot trong Moving average Real :: Sumx,Corr,prob Real, Allocatable :: xx1d(:), yy1d(:) Logical :: Calc_mtx, Normalize, SVD_MatLab, SDV_MatLab, IsAnomaly !-------------------------------------------------------------------- Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab, Dir_Out !Character (Len=15) :: Dat_src="Ens" ! "APH" ! "Ens" Character (Len=15) :: Dat_src="Stn" ! "Ens" Character (Len=15) :: varn="HWD" Character (Len=3) :: Normalize_Opt Character (Len=20), Allocatable :: Stname(:) !-------------------------------------------------------------------- ! LAPACK: Integer :: M, N, LDA, LDU, LDVT, LWMAX, INFO, LWORK DOUBLE PRECISION, Allocatable :: AA(:,:), UU(:,:), VT(:,:), SS(:), WORK(:) EXTERNAL :: DGESVD !-------------------------------------------------------------------- ! Mov_Ave_Step = 0 ! Smooth with moving average ! print*,"Set Options ..." ! Normalize = .false. !.true. ! SDV_MatLab = .true. ! Calc_mtx = .true. ! IsAnomaly = .true. ! Di thuong = Origin - Monthly mean Include "MCA_Options.par" if (Normalize) then Normalize_Opt="Cor"; else Normalize_Opt="Cov"; endif print*,"Normalize,Calc_mtx,SVD_MatLab,IsAnomaly: ",Normalize,Calc_mtx,SVD_MatLab,IsAnomaly !-------------------------------------------------------------------- Fname="Heatwave.txt" Call Read_HW_data (Fname) Fname="/work/projects/G09/SEACLID/ClimIDX/Clim_IDX_2022.nc" Call Read_ClimIDX (Fname) !,varn) Nt=Ntime; Nx=Landpoint; Ny=Nidx Nst=Landpoint Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(Ny), Sx(Nx), Sy(Ny)) !Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nt,Nx), B(Nt,Ny)) Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt)) Allocate (WW(Ny,Ny)) Allocate (Uk(Nx,Nmode), Vk(Ny,Nmode), WWk(Nmode,Nmode), Wk(Nmode)) !Allocate (Ak(Nt,Nmode), Bk(Nt,Nmode), Corre(Nmode),Probs(Nmode)) Allocate (Ak(Nmode,Nt), Bk(Nmode,Nt), Corre(Nmode),Probs(Nmode)) Allocate (Accum(Nmode),Percent(Nmode)) Allocate (U3D(Nlon,Nlat,Nmode),V3D(Nlon_sst,Nlat_sst,Nmode)) Allocate (XX(Nt),YY(Nt)) !-------------------------------------------------------------------- ! print*,"Centering data ..." ! Normalize = .false. !.true. ! SDV_MatLab = .true. ! Calc_mtx = .true. ! IsAnomaly = .true. ! Di thuong = Origin - Monthly mean !print*,Normalize,SDV_MatLab,Calc_mtx,IsAnomaly; pause ! if (Normalize) then ! Normalize_Opt="Cor"; ! else ! Normalize_Opt="Cov"; ! endif !-------------------------------------------------------------------- Call Centered_Data (Normalize) ! Su dung de tinh Ak, Bk Dir_Out="/home/tanpv/Progs/HW_Ha/NC/" if (Calc_mtx) then ! Chi tinh Ma tran Covariance ! Covariance matrix: print*,"Calculating Covariance matrix..." C = MatMul(Xc,Transpose(Yc))/Real(Nt-1) ! Write C to nc file: Allocate (xx1d(Nx), yy1d(Ny)) Do i=1,Nx; xx1d(i)=i; Enddo Do i=1,Ny; yy1d(i)=i; Enddo Call wr_ncdf_2d (trim(Dir_Out)//"C_Matrix.nc",xx1d,yy1d, & Nx,Ny,C,"cor","-", "-", "x", "y", Miss) Deallocate (xx1d,yy1d) ! Calculating SVD: if (SVD_MatLab) then isys = system ("./run_mlab.sh ") Fnu=trim(Dir_Out)//"U.nc" Fnv=trim(Dir_Out)//"V.nc" Fnab=trim(Dir_Out)//"S.nc" Call Read_UVW_MatLab (Fnu,Fnv,Fnab) else isys = system ("./run_r.sh") Fnu=trim(Dir_Out)//"u.nc" Fnv=trim(Dir_Out)//"v.nc" Fnab=trim(Dir_Out)//"w.nc" Call Read_UVW (Fnu,Fnv,Fnab) endif Call Calc_3D_AB !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ else ! Xu ky KQ sau khi thuc hien SVD bang R if (SVD_MatLab) then Fnu=trim(Dir_Out)//"U.nc" Fnv=trim(Dir_Out)//"V.nc" Fnab=trim(Dir_Out)//"S.nc" Call Read_UVW_MatLab (Fnu,Fnv,Fnab) else Fnu=trim(Dir_Out)//"u.nc" Fnv=trim(Dir_Out)//"v.nc" Fnab=trim(Dir_Out)//"w.nc" Call Read_UVW (Fnu,Fnv,Fnab) endif Call Calc_3D_AB ! print test whether U and V are orthogonal: print"(6F10.2)", MatMul(Transpose(Vk),Vk) endif !-------------------------------------------------------- ! Write to NC files: print*,"Writing to file ..." Fnu=trim(Dir_Out)//"U_HW_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnv=trim(Dir_Out)//"V_HW_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt)//".nc" Fnab=trim(Dir_Out)//"AB_HW_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt)//".nc" Call Write_to_NC (Fnu,Fnv,Fnab) ! Calculate and write Homogenerous and Heterogenerous correlation maps: Call Homo_And_Hetero_Map Fname = trim(Dir_Out)//"Ho_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fname1 = trim(Dir_Out)//"He_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnu = trim(Dir_Out)//"pHo_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Fnv = trim(Dir_Out)//"pHe_Cli_"//trim(Dat_src)//"_"//trim(varn)//"_"//trim(Normalize_Opt) Call Write_Homo_Hete_Maps_to_NC (Fname,Fname1,Fnu,Fnv) !-------------------------------------------------------- print*,"Wk", Wk print*,"Corr=", Corre print*,"Prob=",Probs print*,"Perc=",Percent print*,"Accum",Accum CONTAINS !----------------------------------------------------------------- Subroutine Read_UVW_MatLab (Fu,Fv,Fw) Character (Len=*) :: Fu,Fv,Fw Real, Allocatable :: UU(:,:) Real, Allocatable :: xlo(:), xla(:) Integer :: nlon,nlat,i,j Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S") print*, nlon,nlat,nx,ny !Do i=1,Ny; W(i) = SQRT(UU(i,i)); Enddo Do i=1,Ny; W(i) = UU(i,i); Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"U") print*, nlon,nlat,nx,ny Do i=1,Nx U(i,1:Ny) = UU(i,1:Ny) !*W(i) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"V") print*, nlon,nlat,nx,ny Do i=1,Ny V(i,1:Ny) = UU(i,1:Ny) !*W(i) Enddo Deallocate (xlo,xla,UU) End Subroutine Read_UVW_MatLab !----------------------------------------------------------------- Subroutine Read_UVW (Fu,Fv,Fw) Character (Len=*) :: Fu,Fv,Fw Real, Allocatable :: UU(:,:) Real, Allocatable :: xlo(:), xla(:) Integer :: nlon,nlat,i,j Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny Do i=1,Nx U(i,1:Ny) = UU(i,1:Ny) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny Do i=1,Ny V(i,1:Ny) = UU(1:Ny,i) Enddo Deallocate (xlo,xla,UU) Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"var") print*, nlon,nlat,nx,ny W(1:Ny) = UU(1:Ny,1) Deallocate (xlo,xla,UU) End Subroutine Read_UVW !----------------------------------------------------------------- Subroutine Calc_3D_AB ! Sign inversion of U and V: !U = -U; V = -V Wk(1:Nmode) = W(1:Nmode) ! Singular values Do i=1,Nmode Uk(:,i)=U(:,i); Vk(:,i)=V(:,i) Enddo ! Singular matrix: WW=0; Do j=1,Ny; WW(j,j)=W(j); Enddo ! Full WWk=0; Do j=1,Nmode; WWk(j,j)=W(j); Enddo ! Nmode ! Timeseries of Nmode Ak = MatMul(Transpose(Uk),Xc) Bk = MatMul(Transpose(Vk),Yc) ! Correlations: Do i=1,Nmode XX(1:Nt)=Ak(i,1:Nt); YY(1:Nt)=Bk(i,1:Nt) Call Correlation (XX,YY,Nt,Corr,Prob, Miss) Corre(i)=Corr; Probs(i)=prob Enddo ! Accum (%) of singular values: Sumx=Sum(W(1:Ny)); Accum(:)=0 Do i=1,Nmode Percent(i)=W(i)/Sumx*100 Accum(i)=Sum(Percent(1:i)) Enddo print*,"Calc_3D_AB... Done" End Subroutine Calc_3D_AB !-------------------------------------------------------------------- Subroutine Homo_And_Hetero_Map ! Calculate Homogenerous correlation map = r[X(t), a_k (t)] ! and Heterogenerous correlation map = r[X(t), b_k (t)] Integer :: ilon,ilat,it,i,j Real :: ctmp, ptmp Real, Allocatable :: Ho(:,:),He(:,:),pHo(:,:),pHe(:,:),X1d(:), A1d(:), B1d(:) Allocate (X1d(Nt),A1d(Nt),B1d(Nt),Ho(Nx,Nmode),He(Nx,Nmode)) Allocate (pHo(Nx,Nmode),pHe(Nx,Nmode)) Do i=1,Nmode A1d(1:Nt)=Ak(i,1:Nt); B1d(1:Nt)=Bk(i,1:Nt) Do j=1,Nx X1d(1:Nt)=X(j,1:Nt) Call Correlation(X1d,A1d,Nt,ctmp,ptmp,Miss) Ho(j,i)=ctmp; pHo(j,i)=ptmp Call Correlation(X1d,B1d,Nt,ctmp,ptmp,Miss) He(j,i)=ctmp; pHe(j,i)=ptmp Enddo Enddo ! Convert to 3D: Allocate (Homo(Nx,Nmode),p_Ho(Nx,Nmode)) Allocate (Hete(Nx,Nmode),p_He(Nx,Nmode)) Homo=Ho; Hete=He; p_Ho=pHo; p_He=pHe End Subroutine Homo_And_Hetero_Map !-------------------------------------------------------------------- Subroutine Write_Homo_Hete_Maps_to_NC (FHo,FHe,FpHo,FpHe) Character (Len=*) :: FHo,FHe,FpHo,FpHe Real, Allocatable :: xmod(:) Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Integer :: i,j,k Call Write_GrADS_Stn_Multivar_1 (Homo,Nst,Nmode,FHo) Call Write_GrADS_Stn_Multivar_1 (Hete,Nst,Nmode,FHe) Call Write_GrADS_Stn_Multivar_1 (p_Ho,Nst,Nmode,FpHo) Call Write_GrADS_Stn_Multivar_1 (p_He,Nst,Nmode,FpHe) Call Write_to_TXT (trim(FHo)//".txt",Homo,Nst,Nmode,"Homo") Call Write_to_TXT (trim(FHe)//".txt",Hete,Nst,Nmode,"Hete") Call Write_to_TXT (trim(FpHo)//".txt",p_Ho,Nst,Nmode,"p-val_Ho") Call Write_to_TXT (trim(FpHe)//".txt",p_He,Nst,Nmode,"p-val_He") End Subroutine Write_Homo_Hete_Maps_to_NC !-------------------------------------------------------------------- Subroutine Write_to_NC (Fnu,Fnv,Fnab) Character (Len=*) :: Fnu, Fnv, Fnab ! Filename of U, V and A,B Real, Allocatable :: xmod(:), X2D(:,:), xtime(:), climidx(:) Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Integer :: i,j,k Allocate (xmod(Nmode),X2D(Nmode,Ny),xtime(Nt),climidx(Nidx)) ! print*,"Write_to_TXT >>>>>>>>>>>" Do i=1,Nmode; xmod(i)=i; Enddo Do i=1,Nt; xtime(i)=i; Enddo Do i=1,Nidx; climidx(i)=i; Enddo Call Write_GrADS_Stn_Multivar_1 (Uk,Nst,Nmode,Fnu) Call Write_to_TXT (trim(Fnu)//".txt",Uk,Nst,Nmode,"Uk") ! V matrix: Do i=1,Nmode X2D(i,:)=Vk(:,i) Enddo var_ln="MCA for climate indices" Call wr_ncdf_2d (Fnv,xmod,climidx,Nmode,Nidx, X2D, "v", var_ln, & "-", "mod", "idx", Miss) ! A & B matrix var_ln="Timeseries of MCA of drought index (Ak)" Call wr_ncdf_2d (Fnab,xmod,xtime,Nmode,Nt,Ak, "a", var_ln, & "-", "mod", "mon", Miss) var_ln="Timeseries of MCA of SST (Bk)" Call wr_ncdf_2d_apd (Fnab,Bk,"b",Nmode,Nt, var_ln, "-", & "mod", "mon", Miss) var_ln="Correlation coefficients between A & B" Call wr_ncdf_1d_apd (Fnab,Corre,"cor", Nmode, var_ln,"-","mod", Miss) var_ln="p-values of Correlation coefficients between A & B" Call wr_ncdf_1d_apd (Fnab,probs,"prob", Nmode, var_ln,"-","mod", Miss) var_ln="Percentage of contribution of singular values" Call wr_ncdf_1d_apd (Fnab,Percent,"per", Nmode, var_ln,"-","mod", Miss) var_ln="Accumulated Percent. of contribution of singular values" Call wr_ncdf_1d_apd (Fnab,Accum,"cum", Nmode, var_ln,"-","mod", Miss) End Subroutine Write_to_NC !--------------------------------------------------------------------- Subroutine Write_GrADS_Stn_Multivar_1 (X,Nst,Nmode,Grads_File) ! (Nst,NLen,ndr_idx) Integer :: imon, ivar, ist, ic, i, iidx,Nst,Nmode Character (Len=8) :: STID Character *(*) :: Grads_File Real :: X(Nst,Nmode) !Nst,Nmon,NLen,ndr_idx) ! (Nst,2,Nvar) Real :: RLAT, RLON, TIM, tmp, Var Integer :: NLEV,NFLAG !---------------------------------------------------------------------------------- Open (10,file=trim(Grads_File)//'.dat',FORM='UNFORMATTED',RECORDTYPE='STREAM') TIM = 0.0 ! Co the bien doi tu -0.5 .. +0.5; thuong dat bang 0. NFLAG = 1 ! Khong doi Do i = 1,Nmode ! Nvar chi so, moi chi so duoc xem nhu 1 timeslice NLEV = 1 ! =0 o ban ghi cuoi cung cua 1 timeslice; =1 cho tat ca cac timeslice khac ic = 0 Do ist=1, Nst ic = ic+1 STID = Stname(ist)(1:8) ! Chi lay 8 ky tu dau tien RLAT = lat(ist); RLON = lon(ist) WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Var = X(ist,i) WRITE (10) Var ! Gia tri cua bien Enddo NLEV = 0 WRITE (10) STID,RLAT,RLON,TIM,NLEV,NFLAG Enddo ! ivar Close(10) !---------------------------------------------------------------------! ! Write control file Station data ! !---------------------------------------------------------------------! open (2001,file=trim(Grads_File)//'.ctl') write (2001,'(2A)') 'dset ','^'//trim(Grads_File)//'.dat' write (2001,'(A)') 'dtype station' write (2001,'(2A)') 'stnmap ','^'//trim(Grads_File)//'.dat'//'.map' write (2001,'(A,F8.1)') 'undef ',Miss write (2001,'(A)') 'title Factor Analysis' write (2001,'(A,I3,1x,A)') 'tdef ',Nmode,' linear jan1980 1mo' write (2001,'(A,I3)') 'vars 1' ! ,Nidx write (2001,"(A)") 'uvar 0 99 Drough Chars' write (2001,'(A)') 'endvars' close (2001) End Subroutine Write_GrADS_Stn_Multivar_1 !-------------------------------------------------------------------- Subroutine Correlation(X,Y,N,Corr,prob,Miss) Integer :: N, N1, i, j Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0 Do i=1,N if (X(i)/=Miss.AND.Y(i)/=Miss) then n1=n1+1 Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i) Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i) Rxy=Rxy+X(i)*Y(i) endif Enddo if (N1>2) then Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1) Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb tmp=Rxy/Real(N1) Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy) t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) ) Prob = T_Dist (t, Real(N1-2)) else Rxy=Miss; prob=Miss endif Corr=Rxy End Subroutine Correlation !--------------------------------------------------------------------------- Subroutine Centered_Data (Normalize) Integer :: i,j Logical :: Normalize Real, Allocatable :: Std(:), X1d(:) ! Centering: Allocate (Std(Nx), X1d(Nt)) Do j=1,Nx X1d(1:Nt)=X(j,1:Nt) Xbar(j)=Ave(X1d,Nt,Miss) Std (j)=Stdev_1(X1d,Nt,Miss) Do i=1,Nt Xc(j,i)=X(j,i)-Xbar(j) if (Normalize) Xc(j,i)=Xc(j,i)/Std(j) !/Stdev_1(X(j,1:Nt),Nt,Miss) Enddo Enddo Deallocate (Std); Allocate (Std(Ny)) Do j=1,Ny X1d(1:Nt)=Y(j,1:Nt) Ybar(j)=Ave(X1d,Nt,Miss) Std (j)=Stdev_1(X1d,Nt,Miss) Do i=1,Nt Yc(j,i)=Y(j,i)-Ybar(j) if (Normalize) Yc(j,i)=Yc(j,i)/Std(j) ! Stdev_1(Y(j,1:Nt),Nt,Miss) Enddo Enddo End Subroutine Centered_Data !-------------------------------------------------------------------------------------- Subroutine Read_ClimIDX (Fname) !,varn) Character (Len=*) :: Fname !,varn Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k Allocate (Y(Nidx,Ntime)) Yr11=Yr1 !+1 Do iidx=1,Nidx print*,Ntime,iidx,trim(Clim_nm(iidx)) Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx))) ! j=0 Do iyr=Yr11,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do k=1,Nmon j=j+1; Y(iidx,j)=tmp(i,k) Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) Enddo End Subroutine Read_ClimIDX !--------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N-m+1), tmp Do j=1,N-m+1 tmp=0 Do i=j,j+m-1 tmp=tmp+X(i) Enddo Y(j)=tmp/Real(m) Enddo End Subroutine Moving_Ave !-------------------------------------------------------------------- Subroutine Read_HW_data (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nst=102, Nyrs=Yr2-Yr1+1 Integer :: it, ist, iyr, imon, ireg, id, yr Real :: X3d(Nst,Nyrs,Nmon), X1d(Nyrs*Nmon), xlon,xlat Allocate (Lon(Nst),Lat(Nst),Regs(Nst),Stname(Nst)) Open (1,file=Fname,Status="Old") Read (1,*) ! Tieu de Do ist=1,Nst Do iyr=1,Nyrs Read (1,*) id, yr, (X3d(ist,iyr,imon),imon=1,Nmon), ireg, xlon, xlat if (iyr==1) then; Regs(ist)=ireg; Lon(ist)=xlon; Lat(ist)=xlat; endif Enddo !print*,Regs(ist), Lon(ist), Lat(ist) Write(Stname(ist),"(A,I3.3)") "Stn_",ist Enddo Close (1) Landpoint=Nst Nt=Nyrs*Nmon; Ntime=Nt Allocate (X(Landpoint,Nt)) Do ist=1,Nst it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1; X(ist,it)=X3d(ist,iyr,imon) Enddo Enddo Enddo ! Anomaly: if (IsAnomaly) then Do ist=1,Nst X1d(:)=X(ist,:) Call Anomaly (X1d,Nt) X(ist,:)=X1d(:) Enddo endif End Subroutine Read_HW_data !-------------------------------------------------------------------------------------- Subroutine Write_to_TXT (Fname,X,N,M,Varname) Character (Len=*) :: Fname,Varname Integer :: N,M,i,j Real :: X(N,M) print*,"Write_to_TXT X ",N,M,"--->", trim(Fname) Open (3,file=Fname,Status="Unknown") Write (3,"(A)") "Varname= "//trim(Varname) Do i=1,N Write (3,"(I4,50F10.4)") i,(X(i,j),j=1,M) print*,i Enddo Close (3) print*,"Writing ",trim(Varname)," ... Done !" End Subroutine Write_to_TXT !-------------------------------------------------------------------------------------- Subroutine Anomaly (X,N) Integer :: N, it, j, imon, ic Real :: X(N), Xtb(Nmon), tmp Do imon=1,Nmon tmp=0; ic=0 Do it=1,N j=mod(it,Nmon); if (j==0) j=Nmon if (j==imon) then; ic=ic+1; tmp=tmp+X(it); endif Enddo Xtb(imon)=tmp/ic Enddo Do it=1,N j=mod(it,Nmon); if (j==0) j=Nmon X(it)=X(it)-Xtb(j) Enddo End Subroutine Anomaly !-------------------------------------------------------------------------------------- Logical Function Ar1_Empty(X,N) Integer :: N, i Real :: X(N) Logical :: Empty Empty=.false. Do i=1,N if (X(i)==Miss) then Empty=.true.; exit endif Enddo Ar1_Empty=Empty End Function Ar1_Empty !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Stdev_1.f90" !Include "Sub_MCA.f90" !-------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90"