! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10) Use NetCDF Implicit None Real, Parameter :: Miss=-99.0, Lat_min=-11.0 ! Cat phan Autralia Integer, parameter :: Yr1=1980, Yr2=2018, Nmon=12, Nreg=8 Integer, Parameter :: Yr_sta=1948, Nidx=24 Real, Parameter :: Slon1=60, Slon2=180, Slat1=-20, Slat2=50 ! "drawrec 89 147 -11 27 " Real, Parameter :: lon1_sea=100.0, lon2_sea=110.0, lat1_sea=8.0, lat2_sea=24.0 Real, Parameter :: Miss_quv=-9.96921e+36 !--------------------------------------------------------------------------------- Integer :: Nyear ! Thoi ky Integer :: i,j,L ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: lon(:), lat(:), lev(:), xtime(:) Real, Allocatable :: Mask_VN(:,:), lon_vn(:), lat_vn(:) Real, Allocatable :: sm_vn_reg(:,:,:) ! (Nreg,Ntime,:) Soil moisture Integer :: Nlon, Nlat, Nlev, Ntime, ilon, ilat Integer :: Nlon_vn, Nlat_vn Logical :: Status Real, allocatable :: qu_all(:,:,:), qv_all(:,:,:) ! Nlon,Nlat,Ntime Real, allocatable :: q_flux(:,:) ! Ntime !--------------------------------------------------------------------------------- Character (Len=255) :: Fname !--------------------------------------------------------------------------------- Fname="/home/tanpv/Progs/Dr_VN_2020/VN_PP/NC/quv.mon.mean.nc" Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,Ntime,qu_all,"qu") Deallocate (lon,lat,xtime) Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,Ntime,qv_all,"qv") Fname="/home/tanpv/Progs/Dr_VN_2020/VN_PP/NC/" & //"Moisture_flux.nc" Call Moisture_Flux (Fname) Fname="/home/tanpv/Progs/Dr_VN_2020/VN_PP/NC/" & //"Soil_moisture_reg.nc" Call Soil_Moisture (Fname) CONTAINS !================================================================================= Subroutine Moisture_Flux (Fname) Character (Len=*) :: Fname Integer :: it, i1,i2,j1,j2,ilon,ilat, i,j, imon, imon1, Nc, Nc1 Real :: Qx(Nlon,Nlat), Qy(Nlon,Nlat), tmp(4) Real :: lon1,lat1,lon2,lat2 Real :: Flux, x1,x2,y1,y2,dx,dy,dqx,dqy Integer :: Nx1, Nx2, Ny1, Ny2, Nx, Ny, Move_Ave, M2, Nt Real, Allocatable :: xlon(:),xlat(:), X1d(:), Y1d(:) Real :: Xtb(Nmon), Ytb(Nmon), tmp0, tmp1 Move_Ave=13 ! Moving average M2 = Move_Ave/2; Nt = Ntime-Move_Ave+1 lon1=lon1_sea; lat1=lat1_sea; lon2=lon2_sea; lat2=lat2_sea Do ilon=1,Nlon-1 if (Lon(ilon)<=lon1.AND.Lon(ilon+1)>=lon1) then i1=ilon; exit endif Enddo Do ilon=1,Nlon-1 if (Lon(ilon)<=lon2.AND.Lon(ilon+1)>=lon2) then i2=ilon+1; exit endif Enddo Do ilat=1,Nlat-1 if (Lat(ilat)<=lat1.AND.Lat(ilat+1)>=lat1) then j1=ilat; exit endif Enddo Do ilat=1,Nlat-1 if (Lat(ilat)<=lat2.AND.Lat(ilat+1)>=lat2) then j2=ilat+1; exit endif Enddo !------------------------------------------- !print*,"i1,i2,j1,j2=",i1,i2,j1,j2 Nx1=i1; Nx2=i2; Ny1=j1; Ny2=j2 Nx=Nx2-Nx1+1; Ny=Ny2-Ny1+1 Allocate (xlon(Nx),xlat(Ny)) xlon(1:Nx)=Lon(Nx1:Nx2); xlat(1:Ny)=Lat(Ny1:Ny2) Allocate (q_flux(Ntime,4), X1d(Ntime), Y1d(Ntime)) !q_flux(Ntime,1) : monthly !q_flux(Ntime,2) : monthly moving average !q_flux(Ntime,3) : monthly Anomaly q_flux = Miss; tmp(1)=1; tmp(2)=2; tmp(3)=3; tmp(4)=4 Do it=1,Ntime Qx(:,:) = qu_all(:,:,it); Qy(:,:) = qv_all(:,:,it) Flux=0 Do i=1,Nx-1 lon1=xlon(i); lon2=xlon(i+1) Do j=1,Ny-1 lat1=xlat(j); lat2=xlat(j+1) dx=Distance(lon1,lat1,lon2,lat1) dy=Distance(lon1,lat1,lon1,lat2) dqx=Qx(i+1,j)-Qx(i,j) dqy=Qy(i,j+1)-Qy(i,j) Flux=Flux+(dqy/dx-dqx/dy)*dx*dy Enddo ! pause Enddo q_flux(it,1) = Flux Enddo X1d(1:Ntime) = q_flux(1:Ntime,1) ! Monthly timeseries ! Moving average: Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d) !q_flux(M2+1:Ntime-M2,2) = Y1d(1:Nt) ! Moving average of Monthly timeseries q_flux(1:Ntime,2) = Y1d(1:Ntime) ! Anomaly: Do imon=1,Nmon Nc=0; tmp0=0; Nc1=0; tmp1=0 Do it=1,Ntime imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then if (Y1d(it)/=Miss) then Nc=Nc+1; tmp0=tmp0+Y1d(it) endif if (X1d(it)/=Miss) then Nc1=Nc1+1; tmp1=tmp1+X1d(it) endif endif Enddo Xtb(imon)=tmp0/real(Nc) ! Mean of Monthly Moving average Ytb(imon)=tmp1/real(Nc1) ! Mean of Monthly Enddo Do it=1,Ntime imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (Y1d(it)/=Miss) then q_flux(it,3)=Y1d(it)-Xtb(imon1) ! Monthly anomaly of Moving average of Monthly timeseries endif if (X1d(it)/=Miss) then q_flux(it,4)=X1d(it)-Ytb(imon1) ! Monthly anomaly endif Enddo ! Moving average: !Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d) !q_flux(M2+1:Ntime-M2,2) = Y1d(1:Nt) !Call wr_ncdf_1d (Fname, xtime, Ntime, q_flux, "qflux", & ! "Moisture fluxes", "kg/(m^2s)","time", Miss) Call wr_ncdf_2d (Fname,xtime,tmp,Ntime,4,q_flux,"qflux", & "Moisture fluxes","kg/(m^2s)","time","tmp",Miss) End Subroutine Moisture_Flux !-------------------------------------------------------------------- Subroutine Soil_Moisture (Fout) Character (Len=*) :: Fout Character (Len=255) :: Fname Real, Parameter :: SMMiss=-9.96921e+36 Real, Allocatable :: sm_vn(:,:,:), X1d(:), Y1d(:), time(:) Real, Allocatable :: xtime(:),X3d(:,:,:) Integer :: nlon,nlat, Nt, it1, it2, it, i,ireg Integer :: Move_Ave, imon,imon1, M2 Real :: nc, tmp0,tmp, nc1, tmp1, xtmp(4) Real :: Xtb(Nmon,Nreg), Ytb(Nmon,Nreg), TB Fname="/home/tanpv/Data/AWC/Land_Mask_VN.nc" Call rd_ncdf_2d (Fname,lon_vn,lat_vn,nlon_vn,nlat_vn, & Mask_vn,"all") Where (Mask_vn==2) Mask_vn=1; Where (Mask_vn==3) Mask_vn=2 Where (Mask_vn==4) Mask_vn=3; Where (Mask_vn==5) Mask_vn=4 Where (Mask_vn==6) Mask_vn=5; Where (Mask_vn==7) Mask_vn=6 Where (Mask_vn==8) Mask_vn=7 Deallocate (lon_vn,lat_vn) Fname="/home/tanpv/Data/Soil_moisture/Soil_moisture_VN_0.25.nc" Call rd_ncdf_3d (Fname,lon_vn,lat_vn,xtime,nlon_vn,nlat_vn,nt, & sm_vn,"soilw") Where (sm_vn==SMMiss) sm_vn=Miss Where (IsNaN(sm_vn)) sm_vn=Miss it1=(Yr1-Yr_sta)*Nmon+1 it2=it1+Ntime-1 Allocate (sm_vn_reg(Nreg,Ntime,4),time(Ntime)) sm_vn_reg(:,:,:) = Miss; time(1:Ntime)=xtime(it1:it2) ! 1: Monthly, 2: monthly anomaly, 3: Moving aver, 4: Moving aver anomaly it=0 Do i=it1,it2 it=it+1 Do ireg=1,Nreg-1 nc=0; tmp=0 Do ilon=1,Nlon_vn Do ilat=1,Nlat_vn if (Mask_vn(ilon,ilat)==ireg) then ! subreg if (sm_vn(ilon,ilat,i)/=Miss) then nc=nc+1; tmp=tmp+sm_vn(ilon,ilat,i) endif endif Enddo Enddo if (nc>0) sm_vn_reg(ireg,it,1)=tmp/nc Enddo nc=0; tmp=0 Do ilon=1,Nlon_vn Do ilat=1,Nlat_vn if (Mask_vn(ilon,ilat)>0) then ! VN if (sm_vn(ilon,ilat,i)/=Miss) then nc=nc+1; tmp=tmp+sm_vn(ilon,ilat,i) endif endif Enddo Enddo if (nc>0) sm_vn_reg(Nreg,it,1)=tmp/nc ! Origin monthly regional mean Enddo ! Moving Aver: Move_Ave=13 ! Moving average M2 = Move_Ave/2; Nt = Ntime-Move_Ave+1 Allocate (X1d(Ntime), Y1d(Ntime)) !print*,"Ntime=",Ntime; pause Do ireg=1,Nreg X1d(1:Ntime) = sm_vn_reg(ireg,1:Ntime,1) ! Monthly timeseries Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d) !sm_vn_reg(ireg,M2+1:Ntime-M2,2) = Y1d(1:Ntime) print*,"X1d=", X1d(1:Ntime); pause sm_vn_reg(ireg,1:Ntime,2) = Y1d(1:Ntime) Enddo ! Monthly Climatology: Do ireg=1,Nreg X1d(1:Ntime) = sm_vn_reg(ireg,1:Ntime,1) ! monthly Y1d(1:Ntime) = sm_vn_reg(ireg,1:Ntime,2) ! mov ave print*,"Y1d=sm_vn_reg(ireg,1:Ntime,2)=",Y1d(1:Ntime); pause TB = Ave(Y1d, Ntime, Miss) ! Mean value Do imon=1,Nmon Nc=0; tmp0=0; Nc1=0; tmp1=0 Do it=1,Ntime imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then if (X1d(it)/=Miss) then Nc=Nc+1; tmp0=tmp0+X1d(it) endif if (Y1d(it)/=Miss) then Nc1=Nc1+1; tmp1=tmp1+Y1d(it) endif endif Enddo Xtb(imon,ireg)=tmp0/real(Nc) ! Mean of Monthly Moving average Ytb(imon,ireg)=tmp1/real(Nc1) ! Mean of Monthly Enddo print*,"Xtb=",Xtb(:,ireg) print*,"Ytb=",Ytb(:,ireg) Enddo ! Anomaly: Do ireg=1,Nreg X1d(1:Ntime) = sm_vn_reg(ireg,1:Ntime,1) Y1d(1:Ntime) = sm_vn_reg(ireg,1:Ntime,2) Do it=1,Ntime imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon if (X1d(it)/=Miss) then sm_vn_reg(ireg,it,3)=X1d(it)-Xtb(imon1,ireg) ! Monthly anomaly of Moving average of Monthly timeseries endif if (Y1d(it)/=Miss) then sm_vn_reg(ireg,it,4)=Y1d(it)-Ytb(imon1,ireg) ! Monthly anomaly endif Enddo Enddo ! Write to file: Deallocate (X1d,Y1d) Allocate (X1d(Nreg), X3d(Nreg,4,Ntime)) Do ireg=1,Nreg; X1d(ireg)=ireg; Enddo Do i=1,4; xtmp(i)=i; Enddo Do ireg=1,Nreg Do i=1,4 X3d(ireg,i,1:Ntime)=sm_vn_reg(ireg,1:Ntime,i) Enddo Enddo Call wr_ncdf_3d (Fout,X1d,xtmp,time,Nreg,4,Ntime,X3d,"sm", & "months since 1980-01-01 ", "Soil moisture", "mm", & "reg", "cha","mon", Miss) End Subroutine Soil_Moisture !--------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N), tmp ! Y(N-m+1), tmp k=m/2; Y=Miss Do j=k+1,N-k !1,N-m+1 tmp=0 Do i=j-k,j+k !j,j+m-1 tmp=tmp+X(i) Enddo !Y(j+k)=tmp/Real(m) Y(j)=tmp/Real(m) Enddo End Subroutine Moving_Ave !--------------------------------------------------------------------------------- Logical Function Empty_Ar1(X,N) Integer :: N,i Real :: X(N) Logical :: Empty Empty=.False. Do i=1,N if (X(i)==Miss.OR.IsNaN(X(i)).OR.ABS(X(i))>1000) then Empty=.True.; Exit endif Enddo Empty_Ar1=Empty End Function Empty_Ar1 !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Distance.f90" Include "/home/tanpv/Mylibs/Ave.f90" END !----------------------------------------------------------------------------------!