Use NetCDF Real, parameter :: Miss=-99.0 Integer, parameter :: Nmon=5, Ngrd_point=180, Ntime_mx=180, Nlead=45, Nlw=9, Nw=5 Real, parameter :: Lon1=104.250, Lon2=108.0, Lat1=16.00, Lat2=20.00, Del=0.125 Integer, parameter :: Nlon=Int((Lon2-Lon1)/Del)+1, Nlat=Int((Lat2-Lat1)/Del)+1 Character (Len=50), Parameter :: S2S_file(Nmon)=(/"data_s2s_text_08.txt", "data_s2s_text_09.txt", & "data_s2s_text_10.txt","data_s2s_text_11.txt","data_s2s_text_12.txt"/) Character (Len=50), Parameter :: TRMM_file(Nmon)=(/"data_trmm_text_08.txt","data_trmm_text_09.txt", & "data_trmm_text_10.txt","data_trmm_text_11.txt","data_trmm_text_12.txt"/) Integer :: Ntime Integer :: IO,i1,i2,j,k,it, imon, igrd, ilt, iw Real :: No, xlon, xlat Character (Len=255) :: Fname1, Fname2 Character (Len=10) :: Date Real, dimension (Ngrd_point,Ntime_mx,Nlead,Nmon) :: EC, TRMM ! SL du bao va quan trac Real, dimension (Ngrd_point,Ntime_mx,Nlw,Nmon) :: EC_w, TRMM_w Real, dimension (Ngrd_point,Nlw,Nmon) :: RME0, RMAE0, Cor0, Pval0 Real, dimension (Ngrd_point) :: Lon0, Lat0 Real, dimension (Nmon,Nlw,Ntime_mx) :: Spatial_Cor, Spatial_pv !Real, dimension (Nlon,Nlat) :: Mask Real :: Lon(Nlon), Lat(Nlat) Real, dimension (Nlon,Nlat,Nlw,Nmon) :: RME, RMAE, Cor, Pval Real :: X1d(Nlead), X(Ntime_mx), Y(Ntime_mx) ! Mang lam viec Real :: X1(Ngrd_point), Y1(Ngrd_point) Integer :: Date_FC(Nmon,Ntime_mx), Nspl(Nmon) ! So lan du bao trong thang Real :: tmp1, tmp2, LonMin,LonMax,LatMin,LatMax Logical :: Status !--------------------------------------------------------------------------------- LonMin=999;LonMax=-999;LatMin=999;LatMax=-999; Date_FC=Miss Do imon=1,Nmon Fname1="B4_data/"//trim(S2S_file(imon)) Fname2="B4_data/"//trim(TRMM_file(imon)) Open (1,file=Fname1, status="old"); Open (2,file=Fname2, status="old") ! Read ECMWF IO=0; it=0; Loop1: Do while (IO==0) it=it+1 Do igrd=1,Ngrd_point Read (1,*,IOStat=IO) No, xlon, xlat, Date, (X1d(ilt),ilt=1,Nlead) ! (EC(igrd,it,ilt,imon),ilt=1,Nlead) ! if (xlonLonMax) LonMax=xlon ! if (xlatLatMax) LatMax=xlat if (IO==0) then EC(igrd,it,1:Nlead,imon)=X1d(1:Nlead) if (igrd==1) Date_FC(imon,it)=IDate (Date) ! Convert Date string to Date numbering if (imon==1.AND.it==1) then Lon0(igrd)=xlon; Lat0(igrd)=xlat !print*,Lon(igrd),Lat(igrd) endif else Cycle Loop1 endif Enddo ! print*,LonMin,LonMax,LatMin,LatMax; pause Enddo Loop1 Nspl(imon)=it-1 print*,it; !print*,"Date_FC",Date_FC(imon,1:it-1) ; pause ! Read TRMM IO=0; it=0 Loop2: Do while (IO==0) it=it+1 Do igrd=1,Ngrd_point Read (2,*,IOStat=IO) No, xlon, xlat, Date, (X1d(ilt),ilt=1,Nlead) ! (EC(igrd,it,ilt,imon),ilt=1,Nlead) if (IO==0) then TRMM(igrd,it,1:Nlead,imon)=X1d(1:Nlead) else Cycle Loop2 endif Enddo Enddo Loop2 print*,it; !pause Close (1); Close (2) Enddo ! Convert to weekly (5 days) data: Do imon=1,Nmon Do igrd=1,Ngrd_point Do it=1,Nspl(imon) Do iw=1,Nlw i1=Nw*(iw-1)+1; i2=i1+(Nw-1) tmp1=0; tmp2=0 !print*,i1,i2 Do i=i1,i2 tmp1=tmp1+EC(igrd,it,i,imon); tmp2=tmp2+TRMM(igrd,it,i,imon) Enddo EC_w(igrd,it,iw,imon)=tmp1/Nw; TRMM_w(grd,it,iw,imon)=tmp2/Nw Enddo Enddo Enddo Enddo ! Calculating RME, RMAE, Corr, Pval: Do imon=1,Nmon Do igrd=1,Ngrd_point Do iw=1,Nlw Ntime=Nspl(imon) Do it=1,Ntime X(it)=EC_w(igrd,it,iw,imon); Y(it)=TRMM_w(grd,it,iw,imon) Enddo RME0 (igrd,iw,imon)= Relative_ME(X,Y,Ntime,Miss) RMAE0(igrd,iw,imon)=Relative_MAE(X,Y,Ntime,Miss) Call Correlation(X,Y,Ntime,tmp1,tmp2,Miss) Cor0(igrd,iw,imon)=tmp1; Pval0(igrd,iw,imon)=tmp2 Enddo Enddo Enddo ! Calculating spatial correlation: (Nmon,Nlw,Ntime_mx) :: Spatial_Cor, Spatial_pv Spatial_Cor=Miss; Spatial_pv=Miss Do imon=1,Nmon Ntime=Nspl(imon) Do iw=1,Nlw Do it=1,Ntime X1(1:Ngrd_point)=EC_w (1:Ngrd_point,it,iw,imon) Y1(1:Ngrd_point)=TRMM_w(1:Ngrd_point,it,iw,imon) Call Correlation(X1,Y1,Ngrd_point,tmp1,tmp2,Miss) Spatial_Cor(imon,iw,it)=tmp1 Spatial_pv (imon,iw,it)=tmp2 Enddo Enddo Enddo Call Convert_to_4D Call Write_to_NC ("KQ.nc") Call Write_spatial_cor_to_NC ("Spatial_Corr.nc") !-------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------- Subroutine Write_spatial_cor_to_NC (Fname) Character (Len=*) :: Fname Real :: Lead(Nlw), Xmon(Nmon) Real, Allocatable :: xtime(:), XX(:,:) Character (Len=100) :: t_unit, var_nm, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Do ilt=1,Nlw; Lead(ilt)=ilt; Enddo Allocate (xtime(Ntime_mx), XX(Nlw,Ntime_mx)); xtime=Miss; XX=Miss Do it=1,Ntime_mx; xtime(it)=it; Enddo Do imon=1,Nmon Ntime=Nspl(imon) !Allocate (xtime(Ntime_mx), XX(Nlw,Ntime_mx)); xtime=Miss; XX=Miss !xtime(1:Ntime)=Date_FC(imon,1:Ntime) XX(1:Nlw,1:Ntime)=Spatial_Cor(imon,1:Nlw,1:Ntime) if (imon==1) then write (var_nm,"(A,i2.2)") "cor_",imon+7; var_ln="Spatial correlation"; Call wr_ncdf_2d (Fname,Lead,xtime,Nlw,Ntime_mx,XX,var_nm, & var_ln,"-","leadtime","date_fcs",Miss) write (var_nm,"(A,i2.2)") "pvl_",imon+7; var_ln="P-value of Spatial correlation" Call wr_ncdf_2d_apd (Fname,XX,var_nm,Nlw,Ntime_mx, & var_ln,"-","leadtime","date_fcs",Miss) else write (var_nm,"(A,i2.2)") "cor_",imon+7; var_ln="Spatial correlation"; Call wr_ncdf_2d_apd (Fname,XX,var_nm,Nlw,Ntime_mx, & var_ln,"-","leadtime","date_fcs",Miss) write (var_nm,"(A,i2.2)") "pvl_",imon+7; var_ln="P-value of Spatial correlation" Call wr_ncdf_2d_apd (Fname,XX,var_nm,Nlw,Ntime_mx, & var_ln,"-","leadtime","date_fcs",Miss) endif Enddo End Subroutine Write_spatial_cor_to_NC !-------------------------------------------------------------------- Subroutine Write_to_NC (Fname) Character (Len=*) :: Fname Character (Len=100) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real :: Lead(Nlw), Xmon(Nmon) Integer :: ilt, imon Do ilt=1,Nlw; Lead(ilt)=ilt; Enddo Do imon=1,Nmon; Xmon(imon)=imon+7; Enddo t_unit="days since 2000-01-01 "; lon_nm="lon"; lat_nm="lat"; z_nm="week"; t_nm="mon" var_ln="Sai so trung binh tuong doi"; var_unit="-" Call wr_ncdf_4d (Fname,lon,lat,lead,xmon,nlon,nlat,nlw,nmon,RME,"rme", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) var_ln="Sai so tuong doi trung binh tuyet doi"; var_unit="-" Call wr_ncdf_4d_apd (Fname,RMAE,"rmae", nlon,nlat,nlw,nmon, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) var_ln="He so tuong quan"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Cor,"cor",nlon,nlat,nlw,nmon, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) var_ln="P-value cua He so tuong quan"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Pval,"pval",nlon,nlat,nlw,nmon, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) End Subroutine Write_to_NC !-------------------------------------------------------------------- Subroutine Convert_to_4D Integer :: ilon, ilat, i, j, igrd, imon,iw Real :: xlon,xlat RME=Miss; RMAE=Miss; Cor=Miss; Pval=Miss Do imon=1,Nmon Do iw=1,Nlw Do ilon=1,Nlon xlon=Lon1+(ilon-1)*Del if (imon==1.AND.iw==1) Lon(ilon)=xlon Do ilat=1,Nlat xlat=Lat1+(ilat-1)*Del if (imon==1.AND.iw==1.AND.ilon==1) Lat(ilat)=xlat Do igrd=1,Ngrd_point if (Lon0(igrd)==xlon.AND.Lat0(igrd)==xlat) then RME (ilon,ilat,iw,imon)=RME0 (igrd,iw,imon) RMAE(ilon,ilat,iw,imon)=RMAE0(igrd,iw,imon) Cor (ilon,ilat,iw,imon)=Cor0 (igrd,iw,imon) Pval(ilon,ilat,iw,imon)=Pval0(igrd,iw,imon) endif Enddo Enddo Enddo Enddo Enddo End Subroutine Convert_to_4D !-------------------------------------------------------------------- Real Function Relative_ME (F,O,N,Miss) Integer :: N,i,N1 Real :: F(N), O(N), Miss, tmp tmp=0; N1=0 Do i=1,N if (O(i)/=0) then tmp=tmp+(F(i)-O(i))/O(i); N1=N1+1 endif Enddo if (N1>0) then tmp=tmp/N1 else tmp=Miss endif Relative_ME=tmp End Function Relative_ME !-------------------------------------------------------------------- Real Function Relative_MAE (F,O,N,Miss) Integer :: N,i,N1 Real :: F(N), O(N), Miss, tmp tmp=0; N1=0 Do i=1,N if (O(i)/=0) then tmp=tmp+ABS((F(i)-O(i))/O(i)); N1=N1+1 endif Enddo if (N1>0) then tmp=tmp/N1 else tmp=Miss endif Relative_MAE=tmp End Function Relative_MAE !-------------------------------------------------------------------- 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 !-------------------------------------------------------------------- Integer Function IDate (date) Character (Len=*) :: date Integer :: Lenst, i, j, mydate Character (Len=10) St Lenst=Len(date) St=""; j=0 Do i=1,Lenst if (date(i:i)/=" ".AND.date(i:i)/="-") then; j=j+1; St(j:j)=date(i:i); endif Enddo read (St,"(I8)") mydate IDate=mydate End Function IDate !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" END !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Student_Prob.f90"