Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Nreg=7, Nmon=12, IMiss=-99, NLead=5 Integer, Parameter :: Yr1=1982, Yr2=2010, Ntime0=(Yr2-Yr1+1)*Nmon, Ntime=Ntime0+NLead Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255) :: Dir Character (Len=10) :: Model_Name Real, Allocatable :: lon(:), lat(:), Regs(:,:) !---------------------------------------------------------------------------------------------- ! Input data: Real, Allocatable, dimension(:,:,:,:) :: Pr, T2m, Tx, Tm ! (Nlon,Nlat,Nlead,Ntime) Real, Allocatable, dimension(:,:,:) :: Pr_O, T2m_O, Tx_O, Tm_O ! (Nlon,Nlat,Ntime) Real, Allocatable :: xlon(:), xlat(:), xtime(:) Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) Integer :: Nlon, Nlat, nlo, nla !---------------------------------------------------------------------------------------------- ! Model Climatology: Real, Allocatable, dimension(:,:,:,:) :: Pr_M_cli, T2m_M_cli, Tx_M_cli, Tm_M_cli ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:) :: Pr_O_cli, T2m_O_cli, Tx_O_cli, Tm_O_cli ! (Nlon,Nlat,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Pr_BS, T2m_BS, Tx_BS, Tm_BS ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- ! Differences between Model and Obs: Real, Allocatable, dimension(:,:,:,:) :: Pr_Df, T2_Df, Tx_Df, Tm_Df ! (Nlon,Nlat,Nlead,Ntime) !---------------------------------------------------------------------------------------------- ! Verification Scores: Real, Allocatable, dimension(:,:,:,:) :: Pr_ME, Pr_MAE, Pr_RME, Pr_RMAE, Pr_COR, Pr_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: T2_ME, T2_MAE, T2_COR, T2_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tx_ME, Tx_MAE, Tx_COR, Tx_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tm_ME, Tm_MAE, Tm_COR, Tm_pval ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:) ! (Nreg,Mon,LT,Mem) Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:) ! (Nreg,Mon,Nws,Mem) Integer :: ilon,ilat,it,imon Character (Len=1),Allocatable :: IRegs(:,:) Fname="/home/tanpv/Data/562_VNU/Mask_VN_0.1.nc" Call rd_ncdf_2d (Fname,lon,lat,Nlon,Nlat,Regs,"reg") ! 7 sub-regs mask !Allocate (IRegs(Nlon,Nlat)); Where (Regs==Miss); IRegs=" "; elsewhere; IRegs="*"; End where !do ilat=Nlat,1,-1 ! print"(200A1)",(IRegs(ilon,ilat),ilon=1,Nlon) !enddo !pause Call Read_OBS ! Read obs data of Pr, T2m, Tx, Tm ! Read model data: Dir="/home/tanpv/Data/562_VNU/Rfc_VN/" Model_Name="RegCM" Call Read_Fcs (Dir,Model_Name) ! Calculating Climatology Call Calc_OBS_Cli Call Calc_Model_Cli Call Calc_Bias_cli Fname="../NC/Obs_Model_BS_Cli_"//trim(Model_Name)//".nc" Call Write_Model_Obs_Cli_Bias (Fname) Call Verification Fname="../NC/Veri_"//trim(Model_Name)//".nc" Call Write_Verification (Fname) CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Write_Verification (Fname) Character (Len=*) :: Fname Character (Len=100) :: var_ln, t_unit, var_unit, lon_nm, lat_nm, z_nm, t_nm Real :: Xlead(Nlead), Xmon(Nmon) Integer :: id, imon ! wr_ncdf_4d (filenm,lon,lat,z,t,nlon,nlat,nz,nt,var,var_nm, & ! t_unit, var_long_name, var_unit, lon_nm, lat_nm, z_nm, t_nm, flag) ! wr_ncdf_4d_apd (filenm,var,var_nm, nlon,nlat,nz,nt, & ! var_long_name,var_unit, lon_nm, lat_nm, z_nm, t_nm,flag) Do id=1,Nlead; Xlead(id)=id; Enddo Do imon=1,Nmon; Xmon(imon)=imon; Enddo t_unit="months since 1982-01-01 "; lon_nm="lon"; lat_nm="lat"; z_nm="mon"; t_nm="month" var_ln="Mean Errors (ME) of Precipitation"; var_unit="mm" Call wr_ncdf_4d (Fname,Lon,Lat,Xlead,Xmon,Nlon,Nlat,Nlead,Nmon,Pr_ME,"pr_me",t_unit,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Absolute Errors (MAE) of Precipitation"; var_unit="mm" Call wr_ncdf_4d_apd (Fname,Pr_MAE,"pr_mae",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Relative Mean Errors (RME) of Precipitation"; var_unit="%" Call wr_ncdf_4d_apd (Fname,Pr_RME,"pr_rme",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Relative Mean Absolute Errors (RMAE) of Precipitation"; var_unit="%" Call wr_ncdf_4d_apd (Fname,Pr_RMAE,"pr_rmae",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Correlation between the model and observational precipitation"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Pr_COR,"pr_cor",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="p-values of Correlation between the model and observational precipitation"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Pr_pval,"pr_pval",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Errors (ME) of 2m-Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,T2_ME,"t2_me",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Absolute Errors (RMAE) of 2m-Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,T2_MAE,"t2_mae",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Correlation between the model and observational 2m-Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,T2_COR,"t2_cor",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="p-values of Correlation between the model and observational 2m-Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,T2_pval,"t2_pval",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Errors (ME) of Maximum Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,Tx_ME,"tx_me",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Absolute Errors (RMAE) of Maximum Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,Tx_MAE,"tx_mae",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Correlation between the model and observational Maximum Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Tx_COR,"tx_cor",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="p-values of Correlation between the model and observational Maximum Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Tx_pval,"tx_pval",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Errors (ME) of Minimum Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,Tm_ME,"tm_me",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Mean Absolute Errors (RMAE) of Minimum Temperature"; var_unit="C" Call wr_ncdf_4d_apd (Fname,Tm_MAE,"tm_mae",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="Correlation between the model and observational Minimum Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Tm_COR,"tm_cor",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) var_ln="p-values of Correlation between the model and observational Minimum Temperature"; var_unit="-" Call wr_ncdf_4d_apd (Fname,Tm_pval,"tm_pval",Nlon,Nlat,Nlead,Nmon,var_ln,var_unit,lon_nm,lat_nm,z_nm,t_nm,Miss) End Subroutine Write_Verification !---------------------------------------------------------------------------------------------- Subroutine Verification Integer :: ilat,ilon,id,it,imon,j,ipr,it2,itx,itm, N Real, Allocatable :: F_pr(:),O_pr(:),F_t2(:),O_t2(:),F_tx(:),O_tx(:),F_tm(:),O_tm(:) Real :: Corr, prob N=Ntime/Nmon+1 Allocate (F_pr(N),O_pr(N),F_t2(N),O_t2(N),F_tx(N),O_tx(N),F_tm(N),O_tm(N)) Allocate (Pr_ME(Nlon,Nlat,Nlead,Nmon),Pr_MAE(Nlon,Nlat,Nlead,Nmon),Pr_COR(Nlon,Nlat,Nlead,Nmon)) Allocate (Pr_RME(Nlon,Nlat,Nlead,Nmon),Pr_RMAE(Nlon,Nlat,Nlead,Nmon)) Allocate (T2_ME(Nlon,Nlat,Nlead,Nmon),T2_MAE(Nlon,Nlat,Nlead,Nmon),T2_COR(Nlon,Nlat,Nlead,Nmon)) Allocate (Tx_ME(Nlon,Nlat,Nlead,Nmon),Tx_MAE(Nlon,Nlat,Nlead,Nmon),Tx_COR(Nlon,Nlat,Nlead,Nmon)) Allocate (Tm_ME(Nlon,Nlat,Nlead,Nmon),Tm_MAE(Nlon,Nlat,Nlead,Nmon),Tm_COR(Nlon,Nlat,Nlead,Nmon)) Allocate (Pr_pval(Nlon,Nlat,Nlead,Nmon),T2_pval(Nlon,Nlat,Nlead,Nmon),Tx_pval(Nlon,Nlat,Nlead,Nmon),Tm_pval(Nlon,Nlat,Nlead,Nmon)) Pr_ME=Miss; Pr_MAE=Miss; Pr_COR=Miss; Pr_pval=Miss; Pr_RME=Miss; Pr_RMAE=Miss T2_ME=Miss; T2_MAE=Miss; T2_COR=Miss; T2_pval=Miss; Tx_ME=Miss; Tx_MAE=Miss; Tx_COR=Miss; Tx_pval=Miss; Tm_ME=Miss; Tm_MAE=Miss; Tm_COR=Miss; Tm_pval=Miss; Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)/=Miss) then Do id=1,Nlead Do imon=1,Nmon ipr=0; it2=0; itx=o; itm=0; F_pr=Miss; O_pr=Miss; F_t2=Miss; O_t2=Miss; F_tx=Miss; O_tx=Miss; F_tm=Miss; O_tm=Miss Do it=1,Ntime j=mod(it,Nmon); if (j==0) j=Nmon if (j==imon) then if (Pr(ilon,ilat,id,it)/=Miss.AND.Pr_O (ilon,ilat,it)/=Miss) then ipr=ipr+1; F_pr(ipr)=Pr(ilon,ilat,id,it); O_pr(ipr)=Pr_O (ilon,ilat,it) endif if (T2m(ilon,ilat,id,it)/=Miss.AND.T2m_O(ilon,ilat,it)/=Miss) then it2=it2+1; F_t2(it2)=T2m(ilon,ilat,id,it); O_t2(it2)=T2m_O(ilon,ilat,it) endif if (Tx(ilon,ilat,id,it)/=Miss.AND.Tx_O (ilon,ilat,it)/=Miss) then itx=itx+1; F_tx(itx)=Tx(ilon,ilat,id,it); O_tx(itx)=Tx_O (ilon,ilat,it) endif if (Tm(ilon,ilat,id,it)/=Miss.AND.Tm_O (ilon,ilat,it)/=Miss) then itm=itm+1; F_tm(itx)=Tm(ilon,ilat,id,it); O_tm(itx)=Tm_O (ilon,ilat,it) endif endif Enddo !print"(30F6.1)",F_pr !print"(30F6.1)",O_pr !print"(30F6.1)",(F_pr-O_pr)/O_pr !print* Pr_ME(ilon,ilat,id,imon) =ME (F_pr,O_pr,ipr); Pr_MAE (ilon,ilat,id,imon)=MAE (F_pr,O_pr,ipr) Pr_RME(ilon,ilat,id,imon)=RME(F_pr,O_pr,ipr); Pr_RMAE(ilon,ilat,id,imon)=RMAE(F_pr,O_pr,ipr) Call Correlation(F_pr,O_pr,ipr,Corr,prob,Miss); Pr_COR(ilon,ilat,id,imon)=Corr; Pr_pval(ilon,ilat,id,imon)=prob T2_ME(ilon,ilat,id,imon)=ME(F_t2,O_t2,it2); T2_MAE(ilon,ilat,id,imon)=MAE(F_t2,O_t2,it2) Call Correlation(F_t2,O_t2,it2,Corr,prob,Miss); T2_COR(ilon,ilat,id,imon)=Corr; T2_pval(ilon,ilat,id,imon)=prob Tx_ME(ilon,ilat,id,imon)=ME(F_tx,O_tx,itx); Tx_MAE(ilon,ilat,id,imon)=MAE(F_tx,O_tx,itx) Call Correlation(F_tx,O_tx,itx,Corr,prob,Miss); Tx_COR(ilon,ilat,id,imon)=Corr; Tx_pval(ilon,ilat,id,imon)=prob Tm_ME(ilon,ilat,id,imon)=ME(F_tm,O_tm,itm); Tm_MAE(ilon,ilat,id,imon)=MAE(F_tm,O_tm,itm) Call Correlation(F_tm,O_tm,itm,Corr,prob,Miss); Tm_COR(ilon,ilat,id,imon)=Corr; T2_pval(ilon,ilat,id,imon)=prob Enddo !print"(30F6.1)",Pr_RMAE(ilon,ilat,id,:) !pause Enddo !pause endif Enddo Enddo End Subroutine Verification !---------------------------------------------------------------------------------------------- Subroutine Calc_Differences ! Model minus Obs Integer :: ilat,ilon,id,it Allocate (Pr_Df(Nlon,Nlat,Nlead,Ntime), T2_Df(Nlon,Nlat,Nlead,Ntime), Tx_Df(Nlon,Nlat,Nlead,Ntime), Tm_Df(Nlon,Nlat,Nlead,Ntime)) Pr_Df=Miss; T2_Df=Miss; Tx_Df=Miss; Tm_Df=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)/=Miss) then Do id=1,Nlead Do it=1,Ntime if (Pr(ilon,ilat,id,it) /=Miss.AND.Pr_O (ilon,ilat,it)/=Miss) Pr_Df(ilon,ilat,id,it)=Pr(ilon,ilat,id,it) -Pr_O(ilon,ilat,it) if (T2m(ilon,ilat,id,it)/=Miss.AND.T2m_O(ilon,ilat,it)/=Miss) T2_Df(ilon,ilat,id,it)=T2m(ilon,ilat,id,it)-T2m_O(ilon,ilat,it) if (Tx(ilon,ilat,id,it) /=Miss.AND.Tx_O (ilon,ilat,it)/=Miss) Tx_Df(ilon,ilat,id,it)=Tx(ilon,ilat,id,it) -Tx_O(ilon,ilat,it) if (Tm(ilon,ilat,id,it) /=Miss.AND.Tm_O (ilon,ilat,it)/=Miss) Tm_Df(ilon,ilat,id,it)=Tm(ilon,ilat,id,it) -Tm_O(ilon,ilat,it) Enddo Enddo endif Enddo Enddo End Subroutine Calc_Differences !---------------------------------------------------------------------------------------------- Subroutine Write_Model_Obs_Cli_Bias (Fname) Integer :: imon, id Character (Len=*) :: Fname Character (Len=100) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Real :: xmon(Nmon) ! wr_ncdf_3d_apd (filenm,var,var_nm,nlon,nlat,nt, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag) lon_nm="lon"; lat_nm="lat"; t_nm="mon"; t_unit="months since 1982-01-01 " Do imon=1,Nmon; xmon(imon)=imon; Enddo Allocate (X3d(Nlon,Nlat,Nmon)) ! Obs cli: var_ln="Monthly Prec obs climatology"; var_unit="mm/month" Call wr_ncdf_3d(Fname,lon,lat,xmon,Nlon,Nlat,Nmon,Pr_O_cli,"pr_o",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Monthly T2m obs climatology"; var_unit="Deg. C" Call wr_ncdf_3d_apd (Fname,T2m_O_cli,"t2m_o",Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Monthly Tx obs climatology"; var_unit="Deg. C" Call wr_ncdf_3d_apd(Fname,Tx_O_cli,"tx_o",Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Monthly Tm obs climatology"; var_unit="Deg. C" Call wr_ncdf_3d_apd(Fname,Tm_O_cli,"tm_o",Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) ! Model cli: Do id=1,Nlead X3d(:,:,1:Nmon)=Pr_M_cli(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly model Prec climatology for Leadtime=",id," month"; var_unit="mm/month" write (var_nm,"(A,I1.1)") "pr_m_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=T2m_M_cli(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly model T2m climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "t2m_m_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=Tx_M_cli(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly model Tx climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "tx_m_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=Tm_M_cli(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly model Tm climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "tm_m_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Enddo ! Model Bias: Do id=1,Nlead X3d(:,:,1:Nmon)=Pr_BS(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly bias model Prec climatology for Leadtime=",id," month"; var_unit="mm/month" write (var_nm,"(A,I1.1)") "pr_bs_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=T2m_BS(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly bias model T2m climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "t2m_bs_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=Tx_BS(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly bias model Tx climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "tx_bs_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) X3d(:,:,1:Nmon)=Tm_BS(:,:,id,1:Nmon) write (var_ln,"(A,I1.1,A)") "Monthly bias model Tm climatology for Leadtime=",id," month"; var_unit="Deg. C" write (var_nm,"(A,I1.1)") "tm_bs_",id Call wr_ncdf_3d_apd(Fname,X3d,var_nm,Nlon,Nlat,Nmon,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Enddo End Subroutine Write_Model_Obs_Cli_Bias !---------------------------------------------------------------------------------------------- Subroutine Calc_Bias_cli ! Bias climatology Integer :: ilat,ilon,id,imon Allocate (Pr_BS(Nlon,Nlat,Nlead,Nmon), T2m_BS(Nlon,Nlat,Nlead,Nmon), Tx_BS(Nlon,Nlat,Nlead,Nmon), Tm_BS(Nlon,Nlat,Nlead,Nmon)) Pr_BS=Miss; T2m_BS=Miss; Tx_BS=Miss; Tm_BS=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)/=Miss) then Do id=1,Nlead Do imon=1,Nmon Pr_BS (ilon,ilat,id,imon)=Pr_M_cli (ilon,ilat,id,imon)-Pr_O_cli (ilon,ilat,imon) T2m_BS(ilon,ilat,id,imon)=T2m_M_cli(ilon,ilat,id,imon)-T2m_O_cli(ilon,ilat,imon) Tx_BS (ilon,ilat,id,imon)=Tx_M_cli (ilon,ilat,id,imon)-Tx_O_cli (ilon,ilat,imon) Tm_BS (ilon,ilat,id,imon)=Tm_M_cli (ilon,ilat,id,imon)-Tm_O_cli (ilon,ilat,imon) Enddo Enddo endif Enddo Enddo End Subroutine Calc_Bias_cli !---------------------------------------------------------------------------------------------- Subroutine Calc_Model_Cli Integer :: ilat,ilon,id,imon Real :: X1(Ntime), Y(Nmon) Allocate (Pr_M_cli(Nlon,Nlat,Nlead,Nmon), T2m_M_cli(Nlon,Nlat,Nlead,Nmon), Tx_M_cli(Nlon,Nlat,Nlead,Nmon), Tm_M_cli(Nlon,Nlat,Nlead,Nmon)) Pr_M_cli=Miss; T2m_M_cli=Miss; Tx_M_cli=Miss; Tm_M_cli=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)/=Miss) then Do id=1,Nlead X1(1:Ntime)=Pr (ilon,ilat,id,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Pr_M_cli (ilon,ilat,id,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=T2m(ilon,ilat,id,1:Ntime); Call Calc_Cli(X1,Ntime,Y); T2m_M_cli(ilon,ilat,id,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=Tx (ilon,ilat,id,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Tx_M_cli (ilon,ilat,id,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=Tm (ilon,ilat,id,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Tm_M_cli (ilon,ilat,id,1:Nmon)=Y(1:Nmon) Enddo endif Enddo Enddo End Subroutine Calc_Model_Cli !---------------------------------------------------------------------------------------------- Subroutine Calc_OBS_Cli Integer :: ilat,ilon,it,imon Real :: X1(Ntime), Y(Nmon) Allocate (Pr_O_cli(Nlon,Nlat,Nmon), T2m_O_cli(Nlon,Nlat,Nmon), Tx_O_cli(Nlon,Nlat,Nmon), Tm_O_cli(Nlon,Nlat,Nmon)) Pr_O_cli=Miss; T2m_O_cli=Miss; Tx_O_cli=Miss; Tm_O_cli=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)/=Miss) then X1(1:Ntime)=Pr_O (ilon,ilat,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Pr_O_cli (ilon,ilat,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=T2m_O(ilon,ilat,1:Ntime); Call Calc_Cli(X1,Ntime,Y); T2m_O_cli(ilon,ilat,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=Tx_O (ilon,ilat,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Tx_O_cli (ilon,ilat,1:Nmon)=Y(1:Nmon) X1(1:Ntime)=Tm_O (ilon,ilat,1:Ntime); Call Calc_Cli(X1,Ntime,Y); Tm_O_cli (ilon,ilat,1:Nmon)=Y(1:Nmon) endif Enddo Enddo End Subroutine Calc_OBS_Cli !---------------------------------------------------------------------------------------------- Subroutine Calc_Cli (X,N,Y) Integer :: N, i, j, k, imon Real :: X(N), Y(Nmon) Do imon=1,Nmon Y(imon)=0.0; k=0 Do i=1,N j=mod(i,Nmon); if (j==0) j=Nmon if (X(i)/=Miss.AND.j==imon) then k=k+1; Y(imon)=Y(imon)+X(i) endif Enddo if (k/=0) then Y(imon)=Y(imon)/Real(k) else Y(imon)=Miss endif Enddo End Subroutine Calc_Cli !---------------------------------------------------------------------------------------------- Subroutine Read_Fcs (Dir,Model_Name) ! Read forecast data Character (Len=*) :: Dir,Model_Name Integer :: Nla, Nlo, Nt, Yr0, it1, it2, ilead Character (Len=255) Fname !Dir="/home/tanpv/Data/562_VNU/Rfc_VN/" Allocate (Pr(Nlon,Nlat,Nlead,Ntime), T2m(Nlon,Nlat,Nlead,Ntime)) Allocate (Tx(Nlon,Nlat,Nlead,Ntime), Tm(Nlon,Nlat,Nlead,Ntime)) ! (Nlon,Nlat,Nlead,Ntime) Pr=Miss; T2m=Miss; Tx=Miss; Tm=Miss Do ilead=1,NLead write(Fname,"(4A,I1.1,3A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_baseline_VN.nc" !print*,trim(Fname) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"pr") it1=ilead+1; it2=it1+Nt-1 ! Ntime-1 Pr(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tas") T2m(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmax") Tx(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmin") Tm(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Enddo End Subroutine Read_Fcs !---------------------------------------------------------------------------------------------- Subroutine Read_OBS ! Read observed gridded data of Pr, T2m, Tx, Tm Integer :: j, it, ilon, ilat, Nla, Nlo, Nt, Yr0, it1, it2, imon Character (Len=255) Fname Yr0=1980; it1=(Yr1-Yr0)*Nmon+1 ! Lay tu 1/1982 it2=it1+Ntime-1 ! Lay den 5/2011 !print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause Allocate (Pr_O(Nlon,Nlat,Ntime), T2m_O(Nlon,Nlat,Ntime), Tx_O(Nlon,Nlat,Ntime), Tm_O(Nlon,Nlat,Ntime)) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/R.1980.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"rain") Pr_O(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/T2m.1961.2018.mm.nc" Yr0=1961; it1=(Yr1-Yr0)*Nmon+1; it2=it1+Ntime-1 !print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"t2m") T2m_O(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tx.1961.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tx") Tx_O(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tm.1961.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tm") Tm_O(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) End Subroutine Read_OBS !---------------------------------------------------------------------------------------------- Real Function ME(F,O,N) Integer :: N, i, ic Real :: F(N), O(N), tmp tmp=0.0; ic=0 Do i=1,N if (F(i)/=Miss.AND.O(i)/=Miss) then ic=ic+1; tmp=tmp+(F(i)-O(i)) endif Enddo if (ic>N/2) then tmp=tmp/Real(ic) else tmp=Miss endif ME=tmp End Function ME !---------------------------------------------------------------------------------------------- Real Function RME(F,O,N) ! For rainfall Integer :: N, i, ic Real :: F(N), O(N), tmp tmp=0.0; ic=0 Do i=1,N if (O(i)==0.0) O(i)=0.01 ! Gan Rainfall boi gia tri gia (0.01 mm) if (F(i)/=Miss.AND.O(i)/=Miss) then ic=ic+1; tmp=tmp+(F(i)-O(i))/O(i) endif Enddo if (ic>N/2) then tmp=tmp/Real(ic)*100.0 else tmp=Miss endif RME=tmp End Function RME !---------------------------------------------------------------------------------------------- Real Function MAE(F,O,N) Integer :: N, i, ic Real :: F(N), O(N), tmp tmp=0.0; ic=0 Do i=1,N if (F(i)/=Miss.AND.O(i)/=Miss) then ic=ic+1; tmp=tmp+ABS(F(i)-O(i)) endif Enddo if (ic>N/2) then tmp=tmp/Real(ic) else tmp=Miss endif MAE=tmp End Function MAE !---------------------------------------------------------------------------------------------- Real Function RMAE(F,O,N) Integer :: N, i, ic Real :: F(N), O(N), tmp tmp=0.0; ic=0 Do i=1,N if (O(i)==0.0) O(i)=0.01 ! Gan Rainfall boi gia tri gia (0.01 mm) if (F(i)/=Miss.AND.O(i)/=Miss) then ic=ic+1; tmp=tmp+ABS( F(i)/O(i) - 1.0 ) endif Enddo if (ic>N/2) then tmp=tmp/Real(ic)*100.0 else tmp=Miss endif RMAE=tmp End Function RMAE !---------------------------------------------------------------------------------------------- Subroutine wr_ncdf_4d_new (filenm,lon,lat,z,t,nlon,nlat,nz,nt,var,var_nm, & t_unit, var_long_name, var_unit, lon_nm, lat_nm, z_nm, t_nm, flag) character (len = *), intent (in) :: filenm, var_nm, t_unit, var_long_name, var_unit character (len = *) :: lon_nm, lat_nm, z_nm, t_nm integer, parameter :: NDIMS = 4 integer :: ncid, var_id, dims_id(NDIMS) integer :: lon_id, lat_id, z_id, t_id integer, intent (in) :: nlon, nlat, nz, nt real, dimension (nlon,nlat,nz,nt), intent (in) :: var real, intent (in) :: lon(nlon), lat(nlat), flag !, z(nz), t(nt), flag integer, intent (in) :: z(nz), t(nt) print*," Writing netCDF file: ",filenm,"..." call hulk ( nf90_create (filenm, nf90_write, ncid) ) call hulk ( nf90_def_dim (ncid, lon_nm, nlon, dims_id(1)) ) call hulk ( nf90_def_dim (ncid, lat_nm, nlat, dims_id(2)) ) call hulk ( nf90_def_dim (ncid, z_nm, nz, dims_id(3)) ) call hulk ( nf90_def_dim (ncid, t_nm, nt, dims_id(4)) ) call hulk ( nf90_def_var (ncid,lon_nm,nf90_float,dims_id(1),lon_id) ) call hulk ( nf90_def_var (ncid,lat_nm,nf90_float,dims_id(2),lat_id) ) call hulk ( nf90_def_var (ncid,z_nm,NF90_INT,dims_id(3),z_id) ) call hulk ( nf90_def_var (ncid,t_nm,NF90_INT,dims_id(4),t_id) ) call hulk ( nf90_def_var (ncid,var_nm,nf90_float,dims_id(1:4),var_id) ) call hulk ( nf90_put_att (ncid,lon_id,"long_name","Longitude") ) call hulk ( nf90_put_att (ncid,lon_id,"units","degree_east") ) call hulk ( nf90_put_att (ncid,lon_id,"axis","X") ) call hulk ( nf90_put_att (ncid,lat_id,"long_name","Latitude") ) call hulk ( nf90_put_att (ncid,lat_id,"units","degree_north") ) call hulk ( nf90_put_att (ncid,lat_id,"axis","Y") ) call hulk ( nf90_put_att (ncid,z_id, "long_name","z") ) call hulk ( nf90_put_att (ncid,z_id, "units","Lev") ) call hulk ( nf90_put_att (ncid,z_id, "axis","Z") ) call hulk ( nf90_put_att (ncid,t_id,"long_name","Time") ) call hulk ( nf90_put_att (ncid,t_id,"units",trim(t_unit)//" 00:00:00") ) call hulk ( nf90_put_att (ncid,t_id,"calendar","proleptic_gregorian") ) call hulk ( nf90_put_att (ncid,t_id,"axis","T") ) call hulk ( nf90_put_att (ncid,var_id,"long_name",trim(var_long_name))) call hulk ( nf90_put_att (ncid,var_id,"units", trim(var_unit))) call hulk ( nf90_put_att (ncid,var_id,"missing_value",flag) ) call hulk ( nf90_enddef (ncid) ) call hulk ( nf90_put_var (ncid,lon_id,lon) ) call hulk ( nf90_put_var (ncid,lat_id,lat) ) call hulk ( nf90_put_var (ncid,z_id, z) ) call hulk ( nf90_put_var (ncid,t_id, t) ) call hulk ( nf90_put_var (ncid,var_id,var) ) call hulk ( nf90_close (ncid) ) print*," Completed writing file ",filenm End Subroutine wr_ncdf_4d_new !---------------------------------------------------------------------------------------------- Subroutine rd_ncdf_4d_new (filenm,lon,lat,lev,t,nlon,nlat,nlev,nt,var,var_nm) ! use netcdf character (len = *), intent (in) :: filenm, var_nm character (len = 15) :: lon_nm, lat_nm, lev_nm, t_nm integer, parameter :: NDIMS = 4 integer :: ncid, var_id, dims_id(NDIMS) integer :: lon_id, lat_id, lev_id, t_id integer, intent (in out) :: nlon, nlat, nlev,nt real, allocatable, dimension (:,:,:,:), intent (in out) :: var real, allocatable, dimension (:), intent (in out) :: lon, lat !, lev, t integer, allocatable, dimension (:), intent (in out) :: lev, t print*," Reading netCDF file: ",filenm,"..." status = nf90_open (filenm,nf90_Nowrite,ncid) status = nf90_inq_varid (ncid,var_nm,var_id) status = nf90_inquire_variable (ncid,var_id,dimids = dims_id(1:4)) status = nf90_inquire_dimension (ncid, dims_id(1),name = lon_nm, len = nlon) status = nf90_inquire_dimension (ncid, dims_id(2),name = lat_nm, len = nlat) status = nf90_inquire_dimension (ncid, dims_id(3),name = lev_nm, len = nlev) status = nf90_inquire_dimension (ncid, dims_id(4),name = t_nm, len = nt) print*," - Longitude's name = ",lon_nm,"; length = ", nlon print*," - Latitude's name = ",lat_nm,"; length = ", nlat print*," - Level's name = ",lev_nm,"; length = ", nlev print*," - Time's name = ",t_nm,"; length = ", nt ALLOCATE (lon(nlon),lat(nlat),lev(nlev),t(nt),var(nlon,nlat,nlev,nt)) status = nf90_inq_varid (ncid,lon_nm,lon_id) status = nf90_get_var (ncid,lon_id,lon) status = nf90_inq_varid (ncid,lat_nm,lat_id) status = nf90_get_var (ncid,lat_id,lat) status = nf90_inq_varid (ncid,lev_nm,lev_id) status = nf90_get_var (ncid,lev_id,lev) status = nf90_inq_varid (ncid,t_nm,t_id) status = nf90_get_var (ncid,t_id,t) status = nf90_get_var(ncid,var_id,var) status = nf90_close(ncid) print*," Completed reading file ",filenm return End Subroutine rd_ncdf_4d_new !---------------------------------------------------------------------------------- Subroutine Correlation(X,Y,N,Corr,prob,Miss) Integer :: N, N1, i, j Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0 Do i=1,N if (X(i)/=Miss.AND.Y(i)/=Miss) then n1=n1+1 Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i) Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i) Rxy=Rxy+X(i)*Y(i) endif Enddo if (N1>2) then Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1) Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb tmp=Rxy/Real(N1) Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy) t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) ) Prob = T_Dist (t, Real(N1-2)) else Rxy=Miss; prob=Miss endif Corr=Rxy End Subroutine Correlation !--------------------------------------------------------------------------------- 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 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 !----------------------------------------------------------------------- Real function percentile(x, length, per, RMiss) integer length real x(length), per, RMiss real xtos(length),bb,cc,oout integer nn, i nn=0 do i=1, length if(x(i) /= RMiss)then nn=nn+1 xtos(nn)=x(i) endif enddo if(nn.eq.0) then oout=RMiss else call sort(nn,xtos) bb=nn*per+per/3.+1/3. cc=real(int(bb)) if(int(cc).ge.nn) then oout=xtos(nn) else oout=xtos(int(cc))+(bb-cc)*(xtos(int(cc)+1)-xtos(int(cc))) endif endif percentile = oout end function percentile !----------------------------------------------------------------------- SUBROUTINE sort(n,arr) INTEGER n,M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,ir,j,jstack,k,l,istack(NSTACK) REAL a,temp jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 12 j=l+1,ir a=arr(j) do 11 i=j-1,1,-1 if(arr(i).le.a)goto 2 arr(i+1)=arr(i) 11 continue i=0 2 arr(i+1)=a 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp if(arr(l+1).gt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp endif if(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp endif if(arr(l+1).gt.arr(l))then temp=arr(l+1) arr(l+1)=arr(l) arr(l)=temp endif i=l+1 j=ir a=arr(l) 3 continue i=i+1 if(arr(i).lt.a)goto 3 4 continue j=j-1 if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp goto 3 5 arr(l)=arr(j) arr(j)=a jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in sort' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END SUBROUTINE sort !--------------------------------------------------------------------------------- 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/DayOfMon.f90" END !--------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/Student_Prob.f90"