Use netCDF ! Chuong trinh tinh he so tuong quan giua So luong TC tren VES va cac truong KQ-DD cua CFS_Rfc ! Chuoi thoi gian tu 1982-2009 ! Lead time lay tu 1-6 thang ! TC co the la monthly hoac seasonally (tc, tc3m,...). Co the thay TC bang NTC ! Cac truong co the la monthly hoac seasonally. Neu la seasonally thi lay trung binh truot ! Mien tinh la khu vuc TBD-ADD duoc xac dinh boi Lon1, Lon2, Lat1, Lat2 Integer, Parameter :: Nmon=12, Yr1=1982, Yr2=2009, Nyrs=(Yr2-Yr1)+1 Integer, Parameter :: Ntime=Nyrs*Nmon, LT=6 Real, Parameter :: Miss=-99.0, missing_value = 9.999e+20 Real, Parameter :: Dlon=1.0, Dlat=1.0, Lon1=60.0, Lon2=290.0, Lat1=-20.0, lat2=60.0 Integer :: Nst,Nlon,Nlat, ilon,ilat,ist, ilt Real, Allocatable :: sst(:,:,:,:,:) ! (Nlon,Nlat,Nyrs,Nmon,LT) Real, Allocatable :: TC(:,:) ! (Nyrs,Nmon) Real, Allocatable :: xlon(:),xlat(:),xyrs(:),xmon(:) ! Nlon,Nlat,Nyrs,Nmon Character (Len=100) :: Fname, Dir Character (Len=10) :: vname, vartc Character (Len=1) :: Lst !------------------------------------------------------------------------------------------- Real, Allocatable :: Corre(:,:,:,:), Pval(:,:,:,:) ! (Nlon,Nlat,Nmon,LT) Real, Allocatable :: X2D(:,:) ! (Nlon,Nlat) !------------------------------------------------------------------------------------------- Integer, Parameter :: Nvar=9 Integer :: lev, ivar, Mov_Ave ! AO fields: Character (Len=5) :: Vars(Nvar)= & (/"sst","prmsl","vor","vws","hgt","vvel","rh5","rh7","h85"/) !------------------------------------------------------------------------------------------- Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc" vartc="tc3m" Mov_Ave=0 ! =0 neu khong tinh trung binh truot (cho monthly) ! /=0 Tinh trung binh mua (3 hoac 6 thang) ! Reading TC data: Call Read_TC (Fname,trim(vartc)) ! Reading AO fields: Dir="/home/tanpv/Data/CFS_Rfc/ensmean_all/" !--------------------------------------------------------------------------------- Do ivar=1,Nvar vname=trim(Vars(ivar)) print*,"Calculating for ",vname if (ivar <= 4) then ! Reading 3D-arrays (lon,lat,time) Call Read_var (vname,Mov_Ave) else if (ivar<=9) then ! Reading 4D-arrays (lon,lat,lev,time) if (trim(Vars(ivar))=="hgt") then lev=500 Call Read_var_4D (vname,lev,Mov_Ave) else if (trim(Vars(ivar))=="vvel") then lev=850 Call Read_var_4D (vname,lev,Mov_Ave) else if (trim(Vars(ivar))=="rh5") then lev=500 Call Read_var_4D ("rh",lev,Mov_Ave) else if (trim(Vars(ivar))=="rh7") then lev=700 Call Read_var_4D ("rh",lev,Mov_Ave) else if (trim(Vars(ivar))=="h85") then lev=850 Call Read_var_4D ("hgt",lev,Mov_Ave) endif ! Call Read_var_4D (vname,lev,Mov_Ave) !else if (ivar==8) then ! Reading 4D-arrays, multi levs (700,600,500) ! Call Read_var_4D_mlv ("rh",Mov_Ave) !else if (ivar==8) then ! h850 ! lev=850 ! Call Read_var_4D ("hgt",lev,Mov_Ave) endif !Do ilt=1,LT ! Do iyr=1,Nyrs ! print"(I5,12F10.3)",iyr+Yr1-1,sst(20,20,iyr,1:Nmon,ilt) ! Enddo ! pause !Enddo Fname="Corr_"//trim(vartc)//"_"//trim(vname)//"_Rfc.nc" Call Calc_TQ (Fname,vname) Enddo !--------------------------------------------------------------------------------- CONTAINS !------------------------------------------------------------------------------------------- Subroutine Read_TC (Fname,vname) ! vname: tc, tc3m,... Character (Len=*) :: Fname, vname Real, Allocatable :: xyear(:),xmon(:), xtc(:,:) Integer :: Ny, Nm, Nt, it, ix, iy, imon, iyr Call rd_ncdf_2d (Fname,xyear,xmon,Ny,Nm,xtc,trim(vname)) Allocate (TC(Nyrs,Nmon)) iy=0 Do iyr=Yr1,Yr2 ! Chi lay chuoi thoi gian tu Yr1-Yr2 Do it=1,Ny if (int(xyear(it))==iyr) then iy=iy+1; TC(iy,1:Nmon)=xtc(it,1:Nmon) endif Enddo !print"(I5,12F5.1)",iy,(TC(iy,1:Nmon)) Enddo ! print for testing if need: Open (3,file=trim(vartc)//".txt",status="unknown") write (3,"(A4,12I5)")"Yr ",(imon,imon=1,Nmon) Do iyr=1,Nyrs write (3,"(I4.4,12F5.0)") iyr+Yr1-1, (TC(iyr,imon),imon=1,Nmon) Enddo Close (3) End Subroutine Read_TC !------------------------------------------------------------------------------------------- Subroutine Calc_TQ (Fname, vname) Character (Len=*) :: Fname, vname Integer :: Lag1, lag,imtc,imvar,ytc1,ytc2,yva1,yva2 Integer :: iyr, imon,ilon,ilat Character (Len=50) :: var_nm, t_unit, var_ln, var_unit, & lon_nm, lat_nm, z_nm, t_nm Real :: Corr, Prob Real, Allocatable :: X(:), Y(:), xmon(:), xlag(:) Allocate (X(Nyrs),Y(Nyrs)) Allocate (Corre(1:Nlon,1:Nlat,1:Nmon,1:LT), Pval(1:Nlon,1:Nlat,1:Nmon,1:LT)) Do imtc=1,Nmon ! TC mon X(:)=Miss; X(1:Nyrs)=TC(1:Nyrs,imtc) Do lag=1,LT ! Leadtimes !X(:)=Miss; X(1:Nyrs)=TC(1:Nyrs,imtc) Do ilat=1,Nlat Do ilon=1,Nlon Y(:)=Miss; Do iyr=1,Nyrs Y(iyr)=sst(ilon,ilat,iyr,imtc,lag) Enddo Call Correlation (X, Y, Nyrs, Corr, Prob, Miss) Corre(ilon,ilat,imtc,lag)=Corr ! Correlation coefficients Pval (ilon,ilat,imtc,lag)=Prob ! P-values (t-test) Enddo Enddo Enddo Enddo Allocate (xmon(Nmon),xlag(LT)) Do lag=1,LT; xlag(lag)=lag; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo lon_nm="lon"; lat_nm="lat"; z_nm="mon"; t_nm="leadtime"; var_unit="-"; t_unit="months since 1981-01-01 " var_nm="cor"; var_ln="Lag Correlation Coefficients - TC-"//trim(vname) Call wr_ncdf_4d (Fname,xlon,xlat,xmon,xlag,nlon,nlat,nmon,LT,Corre,var_nm, & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) var_nm="pval"; var_ln="P-values of t-test" Call wr_ncdf_4d_apd (Fname,Pval, var_nm, nlon,nlat,nmon,LT, & var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) Deallocate (X,Y,Corre,Pval,xmon,xlag) Deallocate (xlon, xlat, sst) End Subroutine Calc_TQ !------------------------------------------------------------------------------------------- Subroutine Read_var_4D_mlv (vname,Mov_Ave) Character (Len=*) :: vname Integer :: Nt, Nx, Ny, Nlev, it, ix, iy, imon, iyr, ilt, lev,ilev Integer :: ilv1, ilv2, ilv3 Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:), lon(:),lat(:),xt(:) Real, Allocatable :: lonx1(:),latx1(:),lev1(:),xt1(:),XX(:,:,:,:) Real, Allocatable :: X1(:,:,:,:),X2(:,:,:,:),X3(:,:,:,:) Real, Allocatable :: xave(:) Character (Len=250) :: Fname Integer :: Mov_Ave ilv1=Level(700); ilv2=Level(600); ilv3=Level(300) Do ilt=1,LT write (Lst,"(I1)") ilt Fname=trim(Dir)//Lst//"/"//trim(vname)//".1982-2009.lead"//Lst//".nc" Call rd_ncdf_4d (trim(Fname),lonx1,latx1,lev1,xt1,nx,ny,nlev,nt,X,trim(vname)) if (ilt==1) then Allocate (X1(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt)) Allocate (X2(nx,ny,Ntime,LT),X3(nx,ny,Ntime,LT)) Allocate (XX(nx,ny,Ntime,LT)) lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:) endif X1(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv1,1:Nt) ! Chi lay muc ilev X2(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv2,1:Nt) X3(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilv3,1:Nt) Deallocate (lonx1,latx1,lev1,xt1,X) Enddo XX=(X1+X2+X3)/3.0 ! Trung binh 3 truong Deallocate (X1, X2, X3) Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ! Mien TBD-ADD ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime)) ! Lat, lon for plotting Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime)) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Do ilt=1,LT xtmp(:,:,:)=Miss ! missing_value xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(ix:ix+Nlon,iy:iy+Nlat,1:Nt,ilt) Where (xtmp==missing_value) xtmp=Miss ! Convert to my missing values if (Mov_Ave>0) then ! Trung binh truot Do ilat=1,Nlat Do ilon=1,Nlon xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime) Call Moving_Ave (xave,Ntime,Mov_Ave) xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime) Enddo Enddo endif it=0 Do iyr=1,Nyrs ! Resharp to 5D-array (lon,lat,year,mon,leadtime) Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo Enddo !Where (sst==missing_value) sst=Miss Deallocate (xtmp,lon,lat,xt,XX) End Subroutine Read_var_4D_mlv !------------------------------------------------------------------------------------------- Subroutine Read_var_4D (vname,lev,Mov_Ave) Character (Len=*) :: vname Integer :: Nt, Nx, Ny, Nlev, it, ix, iy, imon, iyr, ilt, lev,ilev Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:), lon(:),lat(:),xt(:) Real, Allocatable :: XX(:,:,:,:), lonx1(:),latx1(:),lev1(:),xt1(:) Real, Allocatable :: xave(:) Character (Len=250) :: Fname Integer :: Mov_Ave ilev=Level(lev) Do ilt=1,LT write (Lst,"(I1)") ilt Fname=trim(Dir)//Lst//"/"//trim(vname)//".1982-2009.lead"//Lst//".nc" Call rd_ncdf_4d (trim(Fname),lonx1,latx1,lev1,xt1,nx,ny,nlev,nt,X,trim(vname)) if (ilt==1) then Allocate (XX(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt)) lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:) endif XX(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,ilev,1:Nt) ! Chi lay muc ilev Deallocate (lonx1,latx1,lev1,xt1,X) Enddo Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ! Mien TBD-ADD ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime)) ! Lat, lon for plotting Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime)) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Do ilt=1,LT xtmp(:,:,:)=Miss ! missing_value xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(ix:ix+Nlon,iy:iy+Nlat,1:Nt,ilt) Where (xtmp==missing_value) xtmp=Miss ! Convert to my missing values if (Mov_Ave>0) then ! Trung binh truot Do ilat=1,Nlat Do ilon=1,Nlon xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime) Call Moving_Ave (xave,Ntime,Mov_Ave) xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime) Enddo Enddo endif it=0 Do iyr=1,Nyrs ! Resharp to 5D-array (lon,lat,year,mon,leadtime) Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo Enddo !Where (sst==missing_value) sst=Miss Deallocate (xtmp,lon,lat,xt,XX,xave) End Subroutine Read_var_4D !------------------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m) ! Moving average with m steps and filling up all first X's elements Integer :: N,m,i,j,mj Real :: X(N), Y(N), tmp Y(:)=Miss Do i=1,N-m+1 tmp=0; mj=0 Do j=i,i+m-1 if (X(j)/=Miss) then mj=mj+1; tmp=tmp+X(j) endif Enddo if (mj==m) Y(i)=tmp/mj Enddo X(:)=Y(:) End Subroutine Moving_Ave !------------------------------------------------------------------------------------------- Integer Function Level (lev) Integer :: lev,ilev if (lev==1000) then ilev=1 elseif (lev==925) then ilev=2 elseif (lev==850) then ilev=3 elseif (lev==700) then ilev=4 elseif (lev==600) then ilev=5 elseif (lev==500) then ilev=6 elseif (lev==300) then ilev=7 elseif (lev==200) then ilev=8 elseif (lev==100) then ilev=9 endif Level=ilev End Function Level !------------------------------------------------------------------------------------------- Subroutine Read_var (vname,Mov_Ave) Character (Len=*) :: vname Integer :: Nt, Nx, Ny, it, ix, iy, imon, iyr, ilt Real, Allocatable :: xtmp(:,:,:), X(:,:,:), lon(:),lat(:),xt(:) Real, Allocatable :: XX(:,:,:,:), lonx1(:),latx1(:),xt1(:) Character (Len=250) :: Fname Real, Allocatable :: xave(:) Integer :: Mov_Ave Real :: tmp Do ilt=1,LT write (Lst,"(I1)") ilt Fname=trim(Dir)//Lst//"/"//trim(vname)//".1982-2009.lead"//Lst//".nc" Call rd_ncdf_3d (trim(Fname),lonx1,latx1,xt1,nx,ny,nt, X, trim(vname)) if (ilt==1) then Allocate (XX(nx,ny,Ntime,LT),lon(nx),lat(ny),xt(nt)) lon(:)=lonx1(:); lat(:)=latx1(:); xt(:)=xt1(:) endif XX(1:Nx,1:Ny,1:Nt,ilt)=X(1:Nx,1:Ny,1:Nt) Deallocate (lonx1,latx1,xt1,X) Enddo Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1 ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1 Allocate (xlon(Nlon),xlat(Nlat),xtmp(Nlon,Nlat,Ntime)) Allocate (sst(Nlon,Nlat,Nyrs,Nmon,LT),xave(Ntime)) Allocate (X2D(Nlon,Nlat)) xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat) Do ilt=1,LT xtmp(:,:,:)=Miss ! missing_value xtmp(1:Nlon,1:Nlat,ilt+1:Ntime) = XX(ix:ix+Nlon,iy:iy+Nlat,1:Nt,ilt) Where (xtmp==missing_value) xtmp=Miss if (Mov_Ave>0) then ! Trung binh truot Do ilat=1,Nlat Do ilon=1,Nlon xave(1:Ntime)=xtmp(ilon,ilat,1:Ntime) Call Moving_Ave (xave,Ntime,Mov_Ave) xtmp(ilon,ilat,1:Ntime)=xave(1:Ntime) Enddo Enddo endif it=0 Do iyr=1,Nyrs Do imon=1,Nmon it=it+1 sst(1:Nlon,1:Nlat,iyr,imon,ilt)=xtmp(1:Nlon,1:Nlat,it) Enddo Enddo Enddo if (trim(vname)=="sst") then Where (sst/=Miss) sst=sst-273.15 !missing_value) sst=sst-273.15 else if (trim(vname)=="prmsl") then Where (sst/=Miss) sst=sst/100.0 endif ! Test with SST: !if (trim(vname)=="sst") then ! Call Gradient !endif print*,nx,ny,ntime,Nlon,Nlat,ix,iy print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat) print*,ix,ix+Nlon,iy,iy+Nlat Deallocate (xtmp,lon,lat,xt,XX,X2D) !,xtmp1,xtmp2,ymean) End Subroutine Read_var !------------------------------------------------------------------------------------------- Subroutine Gradient Integer :: ilon,ilat,iyr,imon,ilt Real, Allocatable :: xtmp1(:,:,:,:,:),xtmp2(:,:,:,:,:) Real :: tmp Allocate (xtmp1(Nlon,Nlat,Nyrs,Nmon,LT)) Allocate (xtmp2(Nlon,Nlat,Nyrs,Nmon,LT)) xtmp1=sst; xtmp2=Miss Do ilt=1,LT Do iyr=1,Nyrs Do imon=1,Nmon Do ilon=1,Nlon; Do ilat=1,Nlat X2D(ilon,ilat)=xtmp1(ilon,ilat,iyr,imon,ilt) Enddo; Enddo tmp=Sub_mean (imon) !tmp=xtmp1(69,34,iyr,imon,ilt) if (tmp /= Miss) then Do ilon=1,Nlon; Do ilat=1,Nlat if (xtmp1(ilon,ilat,iyr,imon,ilt)/=Miss) & xtmp2(ilon,ilat,iyr,imon,ilt) = & xtmp1(ilon,ilat,iyr,imon,ilt)-tmp Enddo; Enddo endif Enddo Enddo Enddo sst=xtmp2 Deallocate (xtmp1,xtmp2) End Subroutine Gradient !------------------------------------------------------------------------------------------- Real Function Sub_mean (imon) Character (Len=3) :: SubDomain Integer :: ilon,ilat, ilo1,ilo2,ila1,ila2 Integer :: imon,iyr,ic Real :: tmp, xmean ! Lo1-Lo2 La1-La2 i1-i2 j1-j2 ! PLP: 130-140, 8- 18 (71-81, 29-39) ! VES: 110-118, 10- 18 (51-59, 31-39) ! CHN: 130-140, 20- 30 (71-81, 41-51) ! IND: 80- 90, -5- 5 (21-31, 16- 26) Select Case (imon) Case (1,6,7) SubDomain="PLP" !ilo1=71; ilo2=81; ila1=29; ila2=39 ilo1=68; ilo2=69; ila1=30; ila2=31 Case (2,3,4,5,11,12) SubDomain="PLP" !"VES" !ilo1=51; ilo2=59; ila1=31; ila2=39 ilo1=55; ilo2=56; ila1=34; ila2=35 Case (8) SubDomain="PLP" !"CHN" !ilo1=71; ilo2=81; ila1=41; ila2=51 !ilo1=68; ilo2=70; ila1=41; ila2=43 ilo1=25; ilo2=26; ila1=16; ila2=17 ! IND Case (9,10) SubDomain="PLP" !"IND" !ilo1=21; ilo2=31; ila1=16; ila2=26 ilo1=25; ilo2=26; ila1=16; ila2=17 Case DEFAULT SubDomain="XXX" Sub_mean=0; Return End Select !ilo1=70; ilo2=75; ila1=30; ila2=35 tmp=0.0; ic=0 Do ilat=ila1,ila2 Do ilon=ilo1,ilo2 if (X2D(ilon,ilat) /= Miss) then ic=ic+1; tmp=tmp+X2D(ilon,ilat) endif Enddo Enddo if (ic>0) then xmean=tmp/real(ic) else xmean=Miss endif !print*,ilo1,ilo2,ila1,ila2,xmean; pause Sub_mean=xmean End Function Sub_mean !------------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90" !------------------------------------------------------------------------------------------- END Include "/home/tanpv/Mylibs/Student_Prob.f90" !-------------------------------------------------------------------------------------------