! 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 ! Real, Parameter :: Slon1=50, Slon2=280, Slat1=-40, Slat2=70 ! "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 ! Composite maps of Moisture Transport: Real, allocatable :: q_dry(:,:,:,:),u_dry(:,:,:,:),v_dry(:,:,:,:) ! Q,U,V (Nlon,Nlat,Nlev,Ntime) Real, allocatable :: q_wet(:,:,:,:),u_wet(:,:,:,:),v_wet(:,:,:,:) ! Q,U,V (Nlon,Nlat,Nlev,Ntime) Real, allocatable :: q_all(:,:,:,:),u_all(:,:,:,:),v_all(:,:,:,:) Real, allocatable :: qu_dry(:,:), qu_wet(:,:), qu_all(:,:,:) ! Nlon,Nlat,Ntime Real, allocatable :: qv_dry(:,:), qv_wet(:,:), qv_all(:,:,:) Real, allocatable :: qu_w_d(:,:), qv_w_d(:,:) ! Mean all wet and dry Real, allocatable :: qu_dry_w(:,:), qu_dry_s(:,:) ! Winter, Summer Real, allocatable :: qv_dry_w(:,:), qv_dry_s(:,:) Real, allocatable :: qu_wet_w(:,:), qu_wet_s(:,:) ! Winter, Summer Real, allocatable :: qv_wet_w(:,:), qv_wet_s(:,:) Real, allocatable :: qu_all_w(:,:), qu_all_s(:,:) ! Winter, Summer Real, allocatable :: qv_all_w(:,:), qv_all_s(:,:) ! Winter, Summer Real, allocatable :: q_flux(:,:) ! Ntime Real, allocatable :: slp_win(:,:), slp_sum(:,:), slp_all(:,:,:) Real, allocatable :: u85_win(:,:), v85_win(:,:), u85_sum(:,:), v85_sum(:,:) !--------------------------------------------------------------------------------- Real, allocatable :: Areal_mean_Dr(:) ! (Ntime) !--------------------------------------------------------------------------------- Character (Len=255) :: Fname !--------------------------------------------------------------------------------- Ntime=(Yr2-Yr1+1)*Nmon ! Bo nam dau tien neu khong +1 Call Read_quv_data Fname="/home/tanpv/Progs/2021/NC/Moisture_Trans_moving_ave.nc" !Fname="/home/tanpv/Progs/2021/NC/" & ! //"Moisture_Trans_very_large.nc" Call Moisture_Transport_maps (Fname) CONTAINS !================================================================================= Subroutine Moisture_Transport_maps (Fout) Character (Len=*) :: Fout Character (Len=255) :: Fname, Dir Integer, Parameter :: SWm=11, SSm=5 ! Thang dau mua Dong, mua He Real, Allocatable :: X2d(:,:), xreg(:), time(:), xscal(:), X3d(:,:,:) Real, Allocatable :: Pval(:,:) Real, parameter :: Dr_thres=-0.5, W_thres=0.5 Integer :: Nregs, Nt, Nt1, Nt2, nscal Integer :: ist, i, j, ilon, ilat, ipoint, it Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, imon1, ilev Real :: T_Dist, tmp1, tmp2, ntmp1, ntmp2, Sx, Sy Real :: tslp1, tslp2, tu1, tu2, tv1, tv2 !Fname="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/" & ! //"Areal_mean_and_Trend_APH_spei_12.nc" !Call rd_ncdf_2d (Fname,xreg,time,Nregs,Nt,X2d,"dridx") Dir="/work/projects/G09/SEACLID/GPCC_CRU_SEA/" Fname=trim(Dir)//"SPEI_SPI_GPCC_CRU_Reg_59_19_New.nc" Call rd_ncdf_3d (Fname, xreg, xscal, time, Nregs, nscal, Nt,X3d, "spe12") Allocate (Areal_mean_Dr(Ntime)) Areal_mean_Dr(1:Ntime) = X3d(Nregs,nscal,Nmon+1:Nt) ! Chi lay trung binh toan SEA Deallocate (X3d,xreg,time,xscal) Allocate (qu_dry(Nlon,Nlat),qu_wet(Nlon,Nlat)) Allocate (qv_dry(Nlon,Nlat),qv_wet(Nlon,Nlat)) Allocate (slp_win(Nlon,Nlat),slp_sum(Nlon,Nlat)) ! SLP Allocate (u85_win(Nlon,Nlat),u85_sum(Nlon,Nlat)) ! U850 Allocate (v85_win(Nlon,Nlat),v85_sum(Nlon,Nlat)) Allocate (qu_dry_w(Nlon,Nlat),qu_dry_s(Nlon,Nlat)) Allocate (qv_dry_w(Nlon,Nlat),qv_dry_s(Nlon,Nlat)) Allocate (qu_wet_w(Nlon,Nlat),qu_wet_s(Nlon,Nlat)) Allocate (qv_wet_w(Nlon,Nlat),qv_wet_s(Nlon,Nlat)) ! Clomatological seasons: Allocate (qu_all_w(Nlon,Nlat),qu_all_s(Nlon,Nlat)) Allocate (qv_all_w(Nlon,Nlat),qv_all_s(Nlon,Nlat)) Allocate (qu_w_d(Nlon,Nlat),qv_w_d(Nlon,Nlat)) Nt=Ntime ! Mean All Do ilon=1,Nlon Do ilat=1,Nlat tmp1=0; tmp2=0 Do it=1,Nt tmp1 =tmp1+ qu_all(ilon,ilat,it) tmp2 =tmp2+ qv_all(ilon,ilat,it) Enddo qu_w_d(ilon,ilat)=tmp1/Nt qv_w_d(ilon,ilat)=tmp2/Nt Enddo Enddo ! Wet and dry Do ilon=1,Nlon Do ilat=1,Nlat it1=0; it2=0; tmp1=0; tmp2=0; ntmp1=0; ntmp2=0 !tslp1=0; tslp2=0 Do it=1,Nt if (Areal_mean_Dr(it)<=Dr_thres) then it1=it1+1 tmp1 =tmp1+ qu_all(ilon,ilat,it) ntmp1=ntmp1+qv_all(ilon,ilat,it) ! tslp1=tslp1+slp_all(ilon,ilat,it) ! SLP else if (Areal_mean_Dr(it)>=W_thres) then it2=it2+1 tmp2 =tmp2+qu_all(ilon,ilat,it) ntmp2=ntmp2+qv_all(ilon,ilat,it) ! tslp2=tslp2+slp_all(ilon,ilat,it) ! SLP endif Enddo qu_dry(ilon,ilat)= tmp1/it1; qu_wet(ilon,ilat)= tmp2/it2 qv_dry(ilon,ilat)=ntmp1/it1; qv_wet(ilon,ilat)=ntmp2/it2 !slp_dry(ilon,ilat)=tslp1/it1; slp_wet(ilon,ilat)=tslp2/it2 ! SLP Enddo Enddo ! Seasons for wet and dry: Do ilon=1,Nlon Do ilat=1,Nlat it1=0; it2=0; tmp1=0; tmp2=0; ntmp1=0; ntmp2=0 Do it=1,Nt if (Areal_mean_Dr(it)<=Dr_thres) then ! Dry imon=mod(it,Nmon); if (imon==0) imon=Nmon if (imon>=SWm.OR.imon=SSm.AND.imon=W_thres) then !>Dr_thres) then ! Wet imon=mod(it,Nmon); if (imon==0) imon=Nmon if (imon>=SWm.OR.imon=SSm.AND.imon=SWm.OR.imon=SSm.AND.imon0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_dry, "qv_dry", Nlon,Nlat, & "qv in case Dr_idx<0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_wet, "qv_wet", Nlon,Nlat, & "qv in case Dr_idx>0.5", "g/(m.s)", "lon", "lat", Miss) ! Wet and dry Seasons: Call wr_ncdf_2d_apd (Fout,qu_dry_w, "qu_dry_w", Nlon,Nlat, & "qu_w in case Dr_idx<0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_dry_w, "qv_dry_w", Nlon,Nlat, & "qv_w in case Dr_idx<0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qu_dry_s, "qu_dry_s", Nlon,Nlat, & "qu_s in case Dr_idx<0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_dry_s, "qv_dry_s", Nlon,Nlat, & "qv_s in case Dr_idx<0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qu_wet_w, "qu_wet_w", Nlon,Nlat, & "qu_w in case Dr_idx>0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_wet_w, "qv_wet_w", Nlon,Nlat, & "qu_w in case Dr_idx>0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qu_wet_s, "qu_wet_s", Nlon,Nlat, & "qu_s in case Dr_idx>0.5", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_wet_s, "qv_wet_s", Nlon,Nlat, & "qv_s in case Dr_idx>0.5", "g/(m.s)", "lon", "lat", Miss) ! Climatological seasons: Call wr_ncdf_2d_apd (Fout,qu_all_w, "qu_all_w", Nlon,Nlat, & "qu in winter", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_all_w, "qv_all_w", Nlon,Nlat, & "qv in winter", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qu_all_s, "qu_all_s", Nlon,Nlat, & "qu in summer", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_all_s, "qv_all_s", Nlon,Nlat, & "qv in summer", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,slp_win, "slp_w", Nlon,Nlat, & "slp in winter", "mb", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,slp_sum, "slp_s", Nlon,Nlat, & "slp in summer", "mb", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,u85_win, "u85_w", Nlon,Nlat, & "u850 in winter", "m/s", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,v85_win, "v85_w", Nlon,Nlat, & "v850 in winter", "m/s", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,u85_sum, "u85_s", Nlon,Nlat, & "u850 in summer", "m/s", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,v85_sum, "v85_s", Nlon,Nlat, & "v850 in summer", "m/s", "lon", "lat", Miss) ! All: Call wr_ncdf_2d_apd (Fout,qu_w_d, "qu_all", Nlon,Nlat, & "qu in all cases", "g/(m.s)", "lon", "lat", Miss) Call wr_ncdf_2d_apd (Fout,qv_w_d, "qv_all", Nlon,Nlat, & "qv in all cases", "g/(m.s)", "lon", "lat", Miss) End Subroutine Moisture_Transport_maps !--------------------------------------------------------------------------------- Real Function P_value(Xtb, Ytb, Sx, Sy, Nx, Ny) ! Testing whether Xmean and Ymean are different Real :: t, Xtb, Ytb, Sx, Sy, Nx, Ny, df, Sx2, Sy2 Real :: T_Dist, A, tmp1, tmp2, tmp3, tmp4, pval df=nx+ny-2 if (Xtb/=Miss .AND. Ytb/=Miss .AND. Sx>0 .AND. Sy>0 .AND. Nx>1 .AND. Ny>1) then Sx2=Sx*Sx; Sy2=Sy*Sy ! S^2 tmp1=1.0/real(nx)+1.0/real(ny); tmp1=1.0/SQRT(tmp1) ! 1/Sqrt(1/n1 + 1/n2) tmp2=SQRT(df) tmp3=SQRT((nx-1)*Sx2 + (ny-1)*Sy2) tmp4=tmp3/tmp2 A=tmp1/tmp4 t=ABS(A*(Xtb-Ytb)) pval=T_Dist (t, df) else pval=Miss endif P_value = pval End Function P_value !-------------------------------------------------------------------- Subroutine Read_quv_data Character (Len=255) :: Fname Real, Parameter :: g=9.8, f=100.0 ! m/s2 !Integer, Parameter :: Yr_sta=1891 ! Nam bat dau chuoi SL SST Real, Allocatable :: Q(:,:,:,:),U(:,:,:,:),V(:,:,:,:), SLP(:,:,:) Real, Allocatable :: tlon(:), tlat(:), tlev(:), time(:) Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlo, nla, nle Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, ilev, Nlev1 Character (Len=100) :: var_nm, t_unit, var_ln, var_unit Real :: Delp, tmpu, tmpv, Delu, Delv, Delq Real, Allocatable :: X1d(:), Y1d(:) Integer :: it11, it22, M_ave, Md2, Nave print*, "Reading quv ..." Fname="/work/projects/G09/SEACLID/NCEP/shum.mon.mean.nc" Call rd_ncdf_4d (Fname,tlon,tlat,tlev,time,nlo,nla,nle,nt, & Q,"shum") Nlev=nle ! Lay so muc cua Q Deallocate (tlon,tlat,tlev,time) Fname="/work/projects/G09/SEACLID/NCEP/uwnd.mon.mean.nc" Call rd_ncdf_4d (Fname,tlon,tlat,tlev,time,nlo,nla,nle,nt, & U,"uwnd") Deallocate (tlon,tlat,tlev,time) Fname="/work/projects/G09/SEACLID/NCEP/vwnd.mon.mean.nc" Call rd_ncdf_4d (Fname,tlon,tlat,tlev,time,nlo,nla,nle,nt, & V,"vwnd") Deallocate (tlon,tlat,time) Fname="/work/projects/G09/SEACLID/NCEP/slp.mon.mean.nc" Call rd_ncdf_3d (Fname,tlon,tlat,time,nlo,nla,nt,SLP,"slp") Where (IsNaN(Q)); Q=Miss; Endwhere Where (IsNaN(U)); U=Miss; Endwhere Where (IsNaN(V)); V=Miss; Endwhere Where (Q==Miss_quv); Q=Miss; Endwhere Where (U==Miss_quv); U=Miss; Endwhere Where (V==Miss_quv); V=Miss; Endwhere i=1; Do while (tlon(i)Slat1); j=j+1; Enddo; ilat1=j j=1; Do while (tlat(j)>Slat2); j=j+1; Enddo; ilat2=j it1=(Yr1-Yr_sta)*Nmon+1 ! + Nmon ! Bo nam dau tien neu +Nmon it2=it1+Ntime-1 ! For moving average: M_ave=13; Md2=M_ave/2 it11=it1-Md2; it22=it2+Md2; Nave=it22-it11+1 Allocate (X1d(Nave), Y1d(Ntime)) !print*, "it1, it2, Ntime: ",it1, it2, Ntime !print*, "it11, it22, Nave, Ntime1: ",it11, it22, Nave,Nave-M_ave+1 !pause Nlon=ilon2-ilon1+1; Nlat=ilat1-ilat2+1 ! Chu y: ilat1>ilat2 vi Latitude giam tu North to South Allocate (Lon(Nlon), Lat(Nlat), Lev(Nlev), xtime(Ntime)) Allocate (q_all(Nlon,Nlat,Nlev,Ntime)) Allocate (u_all(Nlon,Nlat,Nlev,Ntime)) Allocate (v_all(Nlon,Nlat,Nlev,Ntime)) Allocate (slp_all(Nlon,Nlat,Ntime)) Lev(1:Nlev)=tlev(1:Nlev) xtime(1:Ntime)=time(it1:it2) i=0 Do ilon=1,Nlo !n if (ilon>=ilon1.AND.ilon<=ilon2) then i=i+1; Lon(i)=tlon(ilon) j=Nlat+1 Do ilat=1,Nla !t if (ilat>=ilat2.AND.ilat<=ilat1) then j=j-1; Lat(j)=tlat(ilat) X1d(1:Nave)=SLP(ilon,ilat,it11:it22) Call Moving_Ave (X1d,Nave,M_ave,Y1d) slp_all(i,j,1:Ntime)=Y1d(1:Ntime) !SLP(ilon,ilat,it1:it2) ! Sea level pressure Do ilev=1,Nlev !u_all(i,j,ilev,1:Ntime)=U(ilon,ilat,ilev,it1:it2) !v_all(i,j,ilev,1:Ntime)=V(ilon,ilat,ilev,it1:it2) !q_all(i,j,ilev,1:Ntime)=Q(ilon,ilat,ilev,it1:it2) X1d(1:Nave)=U(ilon,ilat,ilev,it11:it22) Call Moving_Ave (X1d,Nave,M_ave,Y1d) u_all(i,j,ilev,1:Ntime)=Y1d(1:Ntime) X1d(1:Nave)=V(ilon,ilat,ilev,it11:it22) Call Moving_Ave (X1d,Nave,M_ave,Y1d) v_all(i,j,ilev,1:Ntime)=Y1d(1:Ntime) X1d(1:Nave)=Q(ilon,ilat,ilev,it11:it22) Call Moving_Ave (X1d,Nave,M_ave,Y1d) q_all(i,j,ilev,1:Ntime)=Y1d(1:Ntime) Enddo endif Enddo endif Enddo ! qu_all(:,:,:), qv_all ! Nlon,Nlat,Ntime Allocate (qu_all(Nlon,Nlat,Ntime),qv_all(Nlon,Nlat,Ntime)) ! level = 1000, 925, 850, 700, 600, 500, 400, 300 ; !Nlev1=Nlev-4 ! Skip levels of 400 and 300 mb Nlev1=6 ! 1000, 925, 850, 700, 600, 500 print*,"Lev=", Nlev1,Lev(1:Nlev1) pause Do it=1,Ntime Do ilon=1,Nlon Do ilat=1,Nlat tmpu=0; tmpv=0 Do ilev=1,Nlev1-1 ! Nlev-1 Delp=(Lev(ilev)-Lev(ilev+1))*f Delu=(u_all(ilon,ilat,ilev,it)+u_all(ilon,ilat,ilev+1,it))/2 Delv=(v_all(ilon,ilat,ilev,it)+v_all(ilon,ilat,ilev+1,it))/2 Delq=(q_all(ilon,ilat,ilev,it)+q_all(ilon,ilat,ilev+1,it))/2 tmpu=tmpu+Delu*Delq*Delp tmpv=tmpv+Delv*Delq*Delp Enddo qu_all(ilon,ilat,it)=tmpu/g qv_all(ilon,ilat,it)=tmpv/g Enddo Enddo Enddo !Fname="/home/tanpv/Progs/2021/NC/" & ! //"quv.mon.mean_very_large.nc" Fname="/home/tanpv/Progs/2021/NC/" & //"quv.mon.mean_moving_ave.nc" t_unit="months since 1960-01-01 "; var_ln="U-wind"; var_unit="m/s" Call wr_ncdf_4d (Fname,lon,lat,lev,xtime,nlon,nlat,nlev,Ntime, & u_all,"u",t_unit, var_ln, var_unit, "lon", "lat","lev","time",Miss) var_ln="V-wind"; var_unit="m/s" Call wr_ncdf_4d_apd (Fname,v_all,"v", nlon,nlat,nlev,Ntime, & var_ln,var_unit,"lon", "lat","lev","time",Miss) var_ln="Specific humidity"; var_unit="g/kg" Call wr_ncdf_4d_apd (Fname,q_all,"q", nlon,nlat,nlev,Ntime, & var_ln,var_unit,"lon", "lat","lev","time",Miss) var_ln="Zonal moisture transport"; var_unit="kg/(ms)" Call wr_ncdf_3d_apd (Fname,qu_all,"qu",nlon,nlat,Ntime, var_ln, & var_unit, "lon", "lat", "time", Miss) var_ln="Meridional moisture transport"; var_unit="kg/(ms)" Call wr_ncdf_3d_apd (Fname,qv_all,"qv",nlon,nlat,Ntime, var_ln, & var_unit, "lon", "lat", "time", Miss) Deallocate (U,V,Q,tlat,tlon,tlev,time) End Subroutine Read_quv_data !--------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y) Integer :: N,m, i,j,k ! m must be an odd number Real :: X(N), Y(N-m+1), tmp Do j=1,N-m+1 tmp=0 Do i=j,j+m-1 tmp=tmp+X(i) Enddo 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" END !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/Student_Prob.f90"