! 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=1960, Yr2=2019, Nmon=12 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=89.0, lon2_sea=147.0, lat1_sea=-11.0, lat2_sea=27.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(:) Integer :: Nlon, Nlat, Nlev, Ntime, ilon, ilat Logical :: Status Real, allocatable :: qu_all(:,:,:), qv_all(:,:,:) ! Nlon,Nlat,Ntime Real, allocatable :: q_flux(:,:) ! Ntime !--------------------------------------------------------------------------------- Character (Len=255) :: Fname !--------------------------------------------------------------------------------- Fname="/home/tanpv/Progs/2021/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/2021/NC/Moisture_flux.nc" Call Moisture_Flux (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(5) Real :: lon1,lat1,lon2,lat2 Real :: Flux, x1,x2,y1,y2,dx,dy,dqx,dqy Integer :: Ntt, 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,5), 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; tmp(5)=5 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 ! q_flux(it,1): Monthly timeseries ! q_flux(it,2): Monthly moving average timeseriess ! q_flux(it,3): Monthly Anomaly ! q_flux(it,4): Climatology Anomaly (=Monthly moving average - Mean climatology ! q_flux(it,5): Moving average of Monthly Anomaly 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(M2+1:Ntime-M2) ! Y1d(1:Nt) ! Monthly moving timeseries X1d(1:Ntime) = q_flux(1:Ntime,1) Call Anomaly (X1d,Ntime) ! Anomaly timeseries q_flux(1:Ntime,3) = X1d(1:Ntime) Deallocate (Y1d); Allocate (y1d(Nt)) y1d(1:Nt)=q_flux(M2+1:Ntime-M2,2) tmp1 = Ave(y1d,Nt,Miss) y1d(1:Nt) = y1d(1:Nt) - tmp1 ! Climatology Anomaly q_flux(M2+1:Ntime-M2,4) = y1d(1:Nt) Deallocate (Y1d); Allocate (y1d(Ntime)) X1d(1:Ntime)=q_flux(1:Ntime,3) Call Moving_Ave (X1d,Ntime,Move_Ave,Y1d) q_flux(M2+1:Ntime-M2,5) = Y1d(M2+1:Ntime-M2) ! Moving average of Monthly Anomaly ! 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) ! Ytb(imon)=tmp1/real(Nc1) !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) ! endif ! if (X1d(it)/=Miss) then ! q_flux(it,4)=X1d(it)-Ytb(imon1) ! 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,5,q_flux,"qflux", & "Moisture fluxes: 1:Monthly, 2: Mov. Ave, 3: Ano, 4: Cli Ano", & "kg/(m^2s)","time","tmp",Miss) End Subroutine Moisture_Flux !--------------------------------------------------------------------------------- Real Function Distance(lon1,lat1,lon2,lat2) Real, parameter :: E_radi=6.371E+06 ! in Meters Real :: lon1,lat1,lon2,lat2 ! In degree Real :: phi1,lda1,phi2,lda2 ! In radian Real :: Pi,C,a,delphi,dellda,a1,a2,a3 Pi=4.0*Atan(1.0) phi1=lat1*Pi/180.0; phi2=lat2*Pi/180.0 lda1=lon1*Pi/180.0; lda2=lon2*Pi/180.0 delphi=(phi2-phi1)/2.0; dellda=(lda2-lda1)/2.0 a1=sin(delphi)*sin(delphi) a2=cos(phi1)*cos(phi2) a3=sin(dellda)*sin(dellda) a = a1 + a2*a3 c = 2*atan2(Sqrt(a),Sqrt(1.0-a)) Distance = E_radi*c End Function Distance !------------------------------------------------------------------------------ Subroutine Anomaly (X,Ntime) Integer :: Ntime Real :: X(Ntime), Y(Ntime), TB(Nmon) Integer :: ic, it, imon, imon1 Real :: tmp if (mod(Ntime,Nmon)/=0) print*,"Something wrong !" Do imon=1,Nmon tmp=0; ic=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then if (X(it)/=Miss) then ic=ic+1; tmp=tmp+X(it) endif endif Enddo if (ic>0) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif Enddo Do it=1,Ntime imon=mod(it,Nmon); if (imon==0) imon=Nmon Y(it)=X(it)-TB(imon) Enddo X(:) = Y(:) End Subroutine Anomaly !--------------------------------------------------------------------------------- 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=1,N-m+1 tmp=0 Do i=j,j+m-1 tmp=tmp+X(i) Enddo Y(j+k)=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/Ave.f90" !Include "/home/tanpv/Mylibs/Distance.f90" END !----------------------------------------------------------------------------------!