! Cau lenh doi ngay: ! date -u -d " 1-1-1 00:00:00 17696376 hours " +%Y%m%d%H Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Yr1=1976, Yr2=2005, Nmon=12 Real, Parameter :: Phi1=20, Phi2=23, Lda1=103, Lda2=108 ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon_1(:), xlat_1(:), xtime_1(:) Integer :: Nlon_1, Nlat_1, Ntime_1 Logical :: Status Real, allocatable, dimension(:,:,:) :: olr_1, slp_1, q2_1, rh2_1, t2_1 Real, allocatable, dimension(:,:,:) :: rainc_1, rainnc_1, u10m_1, v10m_1 Real, allocatable, dimension(:,:,:) :: olr, slp, q2, rh2, t2 Real, allocatable, dimension(:,:,:) :: rainc, rainnc, u10m, v10m Real, allocatable, dimension(:,:,:,:) :: dat_all ! (Nlon,Nlat,Ntime,Nvar) Real, allocatable :: xlon(:), xlat(:), xtime(:) Integer :: Nlon, Nlat, Ntime, Nvar !--------------------------------------------------------------------------------- Real, allocatable, dimension(:,:,:) :: PC_slp, PC_q2, PC_rh2, PC_t2 Real, allocatable, dimension(:,:) :: EOF_slp, EOF_q2, EOF_rh2, EOF_t2 !--------------------------------------------------------------------------------- Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:) Real*8, Allocatable :: EOFs(:,:), PCs(:,:) Real*8, Allocatable :: Eigenvalues(:), ECumulate(:) Integer :: Nt, Nx ! Kich thuoc thuc cua mang X Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K Character (Len=250) :: Fname, pca_out, eof_out Real*8 :: Factor Namelist /Model_Info/ pca_out, eof_out, & switch, standard, PC_Stard, Nmode !--------------------------------------------------------------------------------- Open (1,file="Namelist_PCA_New.txt",Status="old") Read (1, Model_Info) Close (1) !--------------------------------------------------------------------------------- Fname="Data/clwrf_25km_2019101600_6hourly.nc" Call Read_Data_wrf (Fname) Call Find_corners (ilon1,ilat1,ilon2,ilat2) print *, "ilon1,ilat1,ilon2,ilat2=", ilon1,ilat1,ilon2,ilat2 !--------------------------------------------------------------------------------- Nt=Ntime; Nvar=4 ! (slp, q2, rh2, t2) Nx=Nlon*Nlat*Nvar Call Subdomain_data !--------------------------------------------------------------------------------- if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif Factor = Sqrt(Real(Nt)) Allocate (xdat(Nt,Nx), ydat(Nt,Nx)) Allocate (Xtb(Nx), Sx(Nx)) ! TB, DLC, ma tran Covar hoac Corre Allocate (Eigenvalues(Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode)) Call Create_dat_for_PCA ! Tinh cac EOFs, PCs Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, & Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) !Allocate (PC_slp(Nlon,Nlat,Nmode), PC_q2(Nlon,Nlat,Nmode)) !Allocate (PC_rh2(Nlon,Nlat,Nmode), PC_t2(Nlon,Nlat,Nmode)) Call Write_EOFs (eof_out) Call Write_PCs (pca_out) !--------------------------------------------------------------------------------- CONTAINS !--------------------------------------------------------------------------------- Subroutine Write_PCs (F_pcs) Character (Len=*) :: F_pcs Character (Len=100) :: var_ln,t_unit Character (Len=3) :: Varnm(Nvar) Real, Allocatable :: PCA(:,:,:), xmode(:) Integer :: ilon,ilat,ivar,imod,ist Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode)) Do imod=1,Nmode; xmode(imod)=imod; Enddo Varnm(1)="slp"; Varnm(2)="q2"; Varnm(3)="rh2"; Varnm(4)="t2" t_unit="hours since 1-1-1 00:00:00" ist=0 Do ivar=1, Nvar Do ilat=1,Nlat Do ilon=1,Nlon ist=ist+1 PCA(ilon,ilat,1:Nmode)=PCs(ist,1:Nmode) Enddo Enddo var_ln="Spatial distribution of PCAs of "//trim(Varnm(ivar)) if (ivar==1) then Call wr_ncdf_3d(F_pcs,xlon,xlat,xmode,nlon,nlat,Nmode, & PCA,Varnm(ivar),trim(t_unit),trim(var_ln), & "-", "lon", "lat","mode", Miss) else Call wr_ncdf_3d_apd(F_pcs,PCA,Varnm(ivar),nlon,nlat,Nmode, & trim(var_ln), "-", "lon", "lat", "mode", Miss) endif Enddo End Subroutine Write_PCs !--------------------------------------------------------------------------------- Subroutine Write_EOFs (F_eof) Character (Len=*) :: F_eof Character (Len=100) :: var_ln, t_unit Real, Allocatable :: xmode(:), xtmp(:), EOF1(:,:,:) integer :: imod, it, i Real :: xlo(1), tmp t_unit="hours since 1-1-1 00:00:00" Allocate (xmode(Nmode), Xtmp(Nmode), EOF1(1,Nmode,Nt)) ! EOF1(Nt,Nmode)) Do imod=1,Nmode; xmode(imod)=imod; Enddo !Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(1:Nt,i); Enddo Do i=1, Nmode; EOF1(1,i,1:Nt)=EOFs(1:Nt,i); Enddo var_ln="Timeseries of eofs of slp, q2, rh2, t2" !Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1,"eof", & ! trim(var_ln), "-", "hour","mode", Miss) xlo(1)=1 Call wr_ncdf_3d (F_eof,xlo,xmode,xtime,1,Nmode,Nt, EOF1, & "eof", trim(t_unit),trim(var_ln), & "-", "lo", "mode","hour", Miss) Xtmp(1:Nmode) = Eigenvalues(1:Nmode) Call wr_ncdf_1d_apd (F_eof, Xtmp,"eig",Nmode,"Eigen values","-", "mode", Miss) tmp=0.0 Do i=1,Mcount tmp=tmp+Eigenvalues(i) if (I<=Nmode) then ECumulate(i)=0 Do j=1,I ECumulate(i)=ECumulate(i)+Eigenvalues(j) Enddo endif Enddo Do i=1,Nmode ECumulate(i)=ECumulate(i)/tmp; Xtmp(i)=Eigenvalues(i)/tmp Enddo Call wr_ncdf_1d_apd (F_eof,Xtmp,"var",Nmode, & "Percent Contribution","-", "mode", Miss) Xtmp(1:Nmode) = ECumulate(1:Nmode) Call wr_ncdf_1d_apd (F_eof,Xtmp,"cum",Nmode, & "Cumuluted Eigenvalues","-", "mode", Miss) End Subroutine Write_EOFs !--------------------------------------------------------------------------------- Subroutine Create_dat_for_PCA Integer :: ilon, ilat, it, ivar, ist ist=0 Do ivar=1, Nvar Do ilat=1,Nlat Do ilon=1,Nlon ist=ist+1 xdat(1:Nt,ist)=dat_all(ilon,ilat,1:Nt,ivar) Enddo Enddo Enddo print*,"ist,Nt,Nx=",ist,Nt,Nx End Subroutine Create_dat_for_PCA !--------------------------------------------------------------------------------- Subroutine Subdomain_data Integer :: ilon, ilat, it, ilo,ila !Ntime=Ntime_1-1 ! Bo timeslice dau tien !Allocate (slp(Nlon,Nlat,Ntime),q2(Nlon,Nlat,Ntime), rh2(Nlon,Nlat,Ntime), t2(Nlon,Nlat,Ntime)) print*,"Nlon,Nlat,Ntime,Nvar=",Nlon,Nlat,Ntime,Nvar Allocate (dat_all(Nlon,Nlat,Ntime,Nvar)) !print*,"Nlon,Nlat,Ntime,Nvar=",Nlon,Nlat,Ntime,Nvar ila=0 Do ilat=ilat1,ilat2 ila=ila+1; ilo=0 Do ilon=ilon1,ilon2 ilo=ilo+1 !slp(ilo,ila,1:Ntime)=slp_1(ilon,ilat,2:Ntime_1) !q2 (ilo,ila,1:Ntime)=q2_1 (ilon,ilat,2:Ntime_1) !rh2(ilo,ila,1:Ntime)=rh2_1(ilon,ilat,2:Ntime_1) !t2 (ilo,ila,1:Ntime)=t2_1 (ilon,ilat,2:Ntime_1) dat_all(ilo,ila,1:Ntime,1)=slp_1(ilon,ilat,2:Ntime_1) dat_all(ilo,ila,1:Ntime,2)=q2_1 (ilon,ilat,2:Ntime_1) dat_all(ilo,ila,1:Ntime,3)=rh2_1(ilon,ilat,2:Ntime_1) dat_all(ilo,ila,1:Ntime,4)=t2_1 (ilon,ilat,2:Ntime_1) Enddo Enddo End Subroutine Subdomain_data !--------------------------------------------------------------------------------- Subroutine Read_Data_wrf (Fname) Character (Len=*) :: Fname Integer :: ilon, ilat, it !Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, olr_1, "olr") !Deallocate (xlon_1, xlat_1, xtime_1) Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, slp_1, "slp") Deallocate (xlon_1, xlat_1, xtime_1) Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, q2_1, "q2") Deallocate (xlon_1, xlat_1, xtime_1) Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, rh2_1, "rh2") Deallocate (xlon_1, xlat_1, xtime_1) Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, t2_1, "t2") !Deallocate (xlon_1, xlat_1, xtime_1) !Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, rainc_1, "rainc") !Deallocate (xlon_1, xlat_1, xtime_1) !Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, rainnc_1, "rainnc") !Deallocate (xlon_1, xlat_1, xtime_1) !Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, u10m_1, "u10m") !Deallocate (xlon_1, xlat_1, xtime_1) !Call rd_ncdf_3d (trim(Fname), xlon_1, xlat_1, xtime_1, nlon_1, nlat_1, Ntime_1, v10m_1, "v10m") Ntime=Ntime_1-1 ! Bo timeslice dau tien Allocate (xtime(Ntime)) xtime(1:Ntime)=xtime_1(2:Ntime_1) ! Mang gia tri thoi gian print*, "nlon_1, nlat_1, Ntime_1=", nlon_1, nlat_1, Ntime_1 End Subroutine Read_Data_wrf !--------------------------------------------------------------------------------- Subroutine Find_corners (ilon1,ilat1,ilon2,ilat2) !Real, Parameter :: Phi1=20, Phi2=23, Lda1=103, Lda2=108 Integer :: ilon,ilat, ilo,ila ilat=0 Do ilat=ilat+1 if (xlat_1(ilat)>=Phi1) then; ilat1=ilat; exit; endif Enddo ilat=nlat_1+1 Do ilat=ilat-1 if (xlat_1(ilat)<=Phi2) then; ilat2=ilat; exit; endif Enddo ilon=0 Do ilon=ilon+1 if (xlon_1(ilon)>=Lda1) then; ilon1=ilon; exit; endif Enddo ilon=nlon_1+1 Do ilon=ilon-1 if (xlon_1(ilon)<=Lda2) then; ilon2=ilon; exit; endif Enddo ! Subdomain size: Nlon=ilon2-ilon1+1; Nlat=ilat2-ilat1+1 ! Subdomain's gridpoints: Allocate (xlon(Nlon),xlat(Nlat)) ilon=0 Do ilo=ilon1,ilon2 ilon=ilon+1; xlon(ilon)=xlon_1(ilo) Enddo ilat=0 Do ila=ilat1,ilat2 ilat=ilat+1; xlat(ilat)=xlat_1(ila) Enddo print*,"Nlon,Nlat=",Nlon,Nlat !print*,"ilon1,ilat1,lon1,lat1: ",ilon1,ilat1,xlon_1(ilon1),xlat_1(ilat1) !print*,"ilon2,ilat2,lon2,lat2: ",ilon2,ilat2,xlon_1(ilon2),xlat_1(ilat2) End Subroutine Find_corners !--------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90" END