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 :: Nseason=4 ! 4 seasons Character (Len=5), Parameter :: Season(Nseason)=(/"12-02","03-05","06-08","09-11"/) !Integer, Parameter :: Yr1=2003, Yr2=2022, Nfcs=18, Nyrs=(Yr2-Yr1+1), Nspl=Nyrs*Nfcs ! 2 thang x 9 lan x 20 nam Integer :: Nfc_per_year(Nseason) ! Obtained from Namelist Integer :: Yr1, Yr2, Nyrs, Nspl, Lw, Nws 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=250) :: Dir_home, Dir_Mask, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_Veri, Dir_Corrected Character (Len=12) :: date0, timeres, idate0, idate1 Character (Len=5), Allocatable :: hdate(:) 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 :: MEe(:,:,:), MAEe(:,:,:), CORe(:,:,:), Prbe(:,:,:) ! (Lon,Lat,LT) !---------------------------------------------------------------------------------------------- Real, Allocatable :: RMEreg(:), RMAEreg(:), CORreg(:) ! (LT) Integer :: imem, ifcs, iy, im, id, NLd, iobs !---------------------------------------------------------------------------------------------- Namelist /Dir_info/ Dir_home,Dir_Mask,DirInEC,DirInObs,Dir_ANN_Save,Dir_Corrected !Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug Namelist /Model_info/ Yr1, Yr2, Lw, Nfc_per_year !---------------------------------------------------------------------------------------------- Open (1, file="Namelist_S2S_ecmwf.txt",Status="old") Read (1, Dir_info) !Read (1, ann_info) Read (1, Model_info) ! Allocate (Layer_size(Nlayer)) ! Read (1, layer_info) Close (1) Nws=NLead-Lw+1 Nyrs=(Yr2-Yr1+1) Call rd_ncdf_2d (trim(Dir_Mask),Lon,Lat,Nlon,Nlat,Regs,"vn") Where (IsNaN(Regs)) Regs=Miss Do isea=1,Nseason ! Calculation for each season among 4 seasons Fname=trim(Dir_home)//"hdate_"//trim(Season(isea))//".txt" Open (1,file=trim(Fname), status="old") Nfcs = Nfc_per_year(isea) ! Number of forecast in each year Nspl = Nfcs*Nyrs ! Total forecast over entire period if (Allocated (hdate)) Deallocate (hdate) Allocate (hdate(Nfcs)) ! Days of forecast Do i=1,Nfcs read (1,*) hdate(i) Enddo Close (1) Call Read_data ! Verification of original ecmwf s2s forecast Call Calc_weekly_data Call Very_Score ! Relative Errors Fname=trim(Dir_Corrected)//"Veri_Score_"//trim(Season(isea))//"_org.nc" Call Write_Very_Scores (Fname) Call Very_Score_1 ! Absolute Errors Fname=trim(Dir_Corrected)//"Veri_Score_"//trim(Season(isea))//"_ME_MAE_org.nc" Call Write_Very_Scores_1 (Fname) Enddo 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)) if (Allocated (R_ME)) Deallocate (R_ME); if (Allocated (R_MAE)) Deallocate (R_MAE) if (Allocated (COR)) Deallocate (COR); if (Allocated (Prb)) Deallocate (Prb) 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 if (Allocated (Ob)) Deallocate (Ob); if (Allocated (Fc)) Deallocate (Fc) 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 Very_Score_1 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)) if (Allocated (MEe)) Deallocate (MEe); if (Allocated (MAEe)) Deallocate (MAEe) if (Allocated (CORe)) Deallocate (CORe); if (Allocated (Prbe)) Deallocate (Prbe) Allocate (MEe(Nlon,Nlat,NLead), MAEe(Nlon,Nlat,NLead), CORe(Nlon,Nlat,NLead),Prbe(Nlon,Nlat,NLead)) MEe=Miss; MAEe=Miss; CORe=Miss; Prbe=Miss if (Allocated (Ob)) Deallocate (Ob); if (Allocated (Fc)) Deallocate (Fc) 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) MEe (ilon,ilat,ilead) = ME (Fc,Ob,Nspl,Miss) !*100. MAEe(ilon,ilat,ilead) = MAE(Fc,Ob,Nspl,Miss) !*100. Enddo endif Enddo Enddo End Subroutine Very_Score_1 !---------------------------------------------------------------------------------------------- Subroutine Write_Very_Scores_1 (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" if (Allocated (lead_day)) Deallocate (lead_day) Allocate (lead_day(Nlead)) Do ilead=1,Nlead; lead_day(ilead)=ilead; Enddo ! Relative Mean Errors: var_ln="Mean Error" Call wr_ncdf_3d (Fname,Lon,Lat,lead_day,Nlon,Nlat,Nlead,MEe,"me",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Mean Absolute Error" Call wr_ncdf_3d_apd (Fname,MAEe,"mae",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_Very_Scores_1 !---------------------------------------------------------------------------------------------- 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" if (Allocated (lead_day)) Deallocate (lead_day) 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: if (Allocated (EC)) Deallocate (EC) 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: if (Allocated (OBS)) Deallocate (OBS) 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) if (Allocated (EC_w)) Deallocate (EC_w) if (Allocated (Obs_w)) Deallocate (Obs_w) 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" !---------------------------------------------------------------------------