Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, IMiss=-99, NLead=46, Lw=5, Nws=NLead-Lw+1 ! Nws=So tuan; Lw = Do dai tuan Integer, Parameter :: Yr1=2001, Yr2=2020, Nfcs=18, Nyrs=(Yr2-Yr1+1), Nspl=Nyrs*Nfcs ! 2 thang x 9 lan x 20 nam Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255) :: DirInEC, DirInObs, DirInStn Character (Len=12) :: date0, timeres, idate0, idate1 Character (Len=5) :: hdate(Nfcs) Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 Real*8 :: factor, offset, Missing Real, Allocatable :: lon(:), lat(:), Regs(:,:) Integer, Allocatable :: lead(:), fcsdate(:,:) Real, Allocatable :: EC(:,:,:,:),OBS(:,:,:,:) ! (Nlon,Nlat,Nspl,NLead) Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:) ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan Real, Allocatable :: xlon(:), xlat(:), xtime(:) Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) Integer :: Nlon, Nlat, date_fcs, Nlon0, Nlat0 ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB Integer :: Nlo, Nla, Ntime, j ! Nday from 2000-01-01 to 2019-12-31 !---------------------------------------------------------------------------------------------- ! Verification Scores: Real, Allocatable :: R_ME(:,:,:), R_MAE(:,:,:), COR(:,:,:), Prb(:,:,:) ! (Lon,Lat,LT) !---------------------------------------------------------------------------------------------- Real, Allocatable :: RMEreg(:), RMAEreg(:), CORreg(:) ! (LT) Integer :: imem, ifcs, iy, im, id, NLd, iobs Fname="/home/tanpv/Data/S2S/2024/Mask_R6_0.125_S2S.nc" Call rd_ncdf_2d (trim(Fname),lon,lat,Nlon,Nlat,Regs,"r6") ! 7 sub-regs mask Fname="/home/cordex_sea/Chi/hdate.txt" Open (1,file=trim(Fname), status="old") Do i=1,Nfcs read (1,*) hdate(i) Enddo Close (1) !print*,"From Mask: ",lat; pause !Deallocate (lon,lat) DirInEC ="/home/cordex_sea/Chi/TayNguyen/EC_TP/" DirInObs ="/home/cordex_sea/Chi/TayNguyen/CHIRPS/" Call Read_data Call Calc_weekly_data Call Very_Score Fname="/home/cordex_sea/Chi/Res/Veri_Score_R6.nc" Call Write_Very_Scores (Fname) CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Very_Score Integer :: ireg,ilon,ilat,imon,ilead,imem,ifcs,im,iobs Real,Allocatable :: Ob(:),Fc(:) !(Nspl) Real :: Corr,prob ! RME(:,:,:), RMAE(:,:,:), COR(:,:,:) ! (Lon,Lat,LT) ! EC_w(Nlon,Nlat,Nspl,Nlead),Obs_w(Nlon,Nlat,Nspl,Nlead)) Allocate (R_ME(Nlon,Nlat,NLead), R_MAE(Nlon,Nlat,NLead), COR(Nlon,Nlat,NLead),Prb(Nlon,Nlat,NLead)) R_ME=Miss; R_MAE=Miss; COR=Miss; Prb=Miss Allocate (Ob(Nspl),Fc(Nspl)) Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ilead=1,Nlead Ob(1:Nspl) = Obs_w(ilon,ilat,1:Nspl,ilead) Fc(1:Nspl) = EC_w (ilon,ilat,1:Nspl,ilead) print"(22F6.1)",Ob(10:30) !print"(22F6.1)",Fc(10:30) R_ME (ilon,ilat,ilead) = RME_1 (Fc,Ob,Nspl,Miss) !*100. R_MAE(ilon,ilat,ilead) = RMAE_1(Fc,Ob,Nspl,Miss) !*100. Call Correlation(Fc,Ob,Nspl,Corr,prob,Miss) COR (ilon,ilat,ilead) = Corr Prb (ilon,ilat,ilead) = prob Enddo R_ME (ilon,ilat,3:Nlead-2) = R_ME (ilon,ilat,3:Nlead-2)*100. R_MAE(ilon,ilat,3:Nlead-2) = R_MAE(ilon,ilat,3:Nlead-2)*100. endif Enddo Enddo End Subroutine Very_Score !---------------------------------------------------------------------------------------------- Subroutine Write_Very_Scores (Fname) Character (Len=*) :: Fname Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_nm Real, Allocatable :: X4d(:,:,:,:),xmon(:),lead_day(:),lead_week(:) Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs,Nspl,Nlead_day t_unit="days since 2001-01-01 "; var_unit="-";lon_nm="lon"; lat_nm="lat"; t_nm="week" Allocate (lead_day(Nlead)) Do ilead=1,Nlead; lead_day(ilead)=ilead; Enddo ! Relative Mean Errors: var_ln="Relative Mean Error" Call wr_ncdf_3d (Fname,Lon,Lat,lead_day,Nlon,Nlat,Nlead,R_ME,"rme",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Relative Mean Absolute Error" Call wr_ncdf_3d_apd (Fname,R_MAE,"rmae",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Correlation between forecast and observed" Call wr_ncdf_3d_apd (Fname,COR, "cor", Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="p-value of Correlation between forecast and observed" Call wr_ncdf_3d_apd (Fname,Prb, "pval",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_Very_Scores !---------------------------------------------------------------------------------------------- Subroutine Areal_Mean (X,Nlon,Nlat,NLead,Y,ireg) Integer :: Nlon,Nlat,NLead, ireg Real :: X(Nlon,Nlat,NLead), Y(NLead) ! EC(ilon,ilat,1:NLead,ifcs,imem) Integer :: ilon,ilat, ilead Real :: tmp, ntmp Y(:)=Miss Do ilead=1,NLead tmp=0; ntmp=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg.AND.X(ilon,ilat,ilead)/=Miss) then ntmp=ntmp+1; tmp=tmp+X(ilon,ilat,ilead) endif Enddo Enddo if (ntmp>0) Y(ilead)=tmp/ntmp Enddo End Subroutine Areal_Mean !---------------------------------------------------------------------------------------------- Subroutine Read_data Real, Parameter :: ec_miss=-32767.0, obs_miss=-9999.0 Integer :: imem, ifcs Real, Allocatable :: xlon(:), xlat(:), xtime(:), X3d(:,:,:) Integer :: Nlo,Nla,NLt,i,ifc,ispl,iyr ! Read precipitation forecast data: Allocate (EC(Nlon,Nlat,Nspl,Nlead)) ispl=0 Do iyr=Yr1,Yr2 Do ifc=1,Nfcs write(Fname,"(A,I4.4,A)") trim(DirInEC)//"LH_SH_OLR_TP_cf_",iyr,"-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"tp") Where (X3d==ec_miss) X3d=Miss ispl=ispl+1 EC(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo ! Read precipitation obs data: Allocate (OBS(Nlon,Nlat,Nspl,Nlead)) ispl=0 Do iyr=Yr1,Yr2 Do ifc=1,Nfcs write(Fname,"(A,I4.4,A)") trim(DirInObs)//"chirps_",iyr,"-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"precip") Where (X3d==obs_miss) X3d=Miss ispl=ispl+1 OBS(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo End Subroutine Read_data !---------------------------------------------------------------------------------------------- Subroutine Calc_weekly_data ! Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:) ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan Integer :: ilon,ilat,ilead,ispl Real :: X(Nlead), Y(Nlead) Allocate (EC_w(Nlon,Nlat,Nspl,Nlead),Obs_w(Nlon,Nlat,Nspl,Nlead)) Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ispl=1,Nspl X(1:Nlead)=EC(ilon,ilat,ispl,1:Nlead) Call Moving_Ave (X,Nlead,Lw,Y,Miss) EC_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead) !print"(22F6.1)",EC_w(ilon,ilat,ispl,10:30) X(1:Nlead)=OBS(ilon,ilat,ispl,1:Nlead) Call Moving_Ave (X,Nlead,Lw,Y,Miss) Obs_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead) Enddo endif Enddo Enddo End Subroutine Calc_weekly_data !---------------------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y,Miss) Integer :: N,m, i,j,k,m2 ! m must be an odd number Real :: X(N), Y(N), A(m), tmp, Miss Y(:) = Miss; m2=m/2 Do j=1,N-m+1 k = j+m-1 A(1:m) = X(j:k) tmp = Ave(A,m,Miss) Y(j+m2)=tmp*m ! Gan cho gia tri giua Enddo End Subroutine Moving_Ave !--------------------------------------------------------------------------------- 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/Sub_Verification.f90" Include "/home/tanpv/Mylibs/Correlation.f90" END Include "/home/tanpv/Mylibs/Student_Prob.f90" !---------------------------------------------------------------------------