Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss = -99. Integer, Parameter :: Nmon = 12 Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6 Integer, Parameter :: Ncase = 3 ! cli,ann Character (Len=3), Parameter :: Cases(Ncase)=(/"cli","ann","mlr"/) Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Pr_Ope(:,:,:,:,:,:), T2_Ope(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel Real, Allocatable :: Pr_Ope_ens(:,:,:), T2_Ope_ens(:,:,:) ! Nlon,Nlat,Nlead Real, Allocatable :: Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:) ! Nlon,Nlat,Nlead = *Cli_Obs + *Ano Real, Allocatable :: Pr_Ano(:,:,:), T2_Ano(:,:,:) ! Nlon,Nlat,Nlead Real, Allocatable :: Pr_Cli_ens(:,:,:,:), T2_Cli_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Nmon Real, Allocatable :: Pr_Cli(:,:,:,:,:), T2_Cli(:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Nmodel Real, Allocatable :: Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:) ! Nlon,Nlat,Nmon Real, Allocatable :: xmem(:), xlead(:), xtime(:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname Character (Len=255) :: Cli_Dir, Dir_Ope, Mask_File, Dir_Out Character (Len=3) :: TVar, RVar Character (Len=10),Allocatable :: Model_nm(:) Logical :: status, Empty Integer :: Yr1, Yr2, YOpe, MOpe, Nmodel, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat Integer :: Nt1 !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d,Xt !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Cli_Dir, Mask_File, Dir_Ope, Dir_Out, TVar, RVar Namelist /Data_Info/ Nmodel Namelist /Model_name/ Model_nm Namelist /Time_Infor/ YOpe, MOpe !------------------------------------------------------------------------------ Open (1, file="Namelist_Anomaly_Pr_T2_Ope.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Model_nm(Nmodel)) Read (1, Model_name) Read (1, Time_Infor) Close (1) Call Read_Regs print*,"Parameters: Nlat,Nlon=",Nlat,Nlon Call Read_Ope_Data ! Read data of "cli","ann","mlr" from Ope run at YOpe, MOpe !Fname=trim(Cli_Dir)//"Clim_Mean_Bias_"//RegCM_pr.nc" Call Read_Clim_Data ! Read Climatological mean of "pr", "t2", "pr_o", "t2_o" Call Do_ens_Ano_Ope ! Ensemble mean over Nspl = Nmodel*Nmem*Ncase = 2*7*3 = 42 corrected members Fname=trim(Dir_Out)//"Forecast_Ope_Pr_T2_and_Anomaly.nc" Call Write_Forecast_Mean_Ano_Pr_T2 (Fname) CONTAINS !------------------------------------------------------------------------------ Subroutine Do_ens_Ano_Ope Integer :: imod,imem,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Integer :: Nspl,islp Real, Allocatable :: xtmp(:), xpr(:) Real :: xfcs, pope, pcli, tope, tcli, pobs, tobs, pmod, tmod ! Pr_Cli(:,:,:,:), T2_Cli(:,:,:,:) ! Nlon,Nlat,Nlead,Nmon ! Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:) ! Nlon,Nlat,Nlead ! Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:) ! Nlon,Nlat,Nmon Nspl=Nmodel*Nmem*Ncase Allocate (xtmp(Nspl),xpr(Nspl)) Allocate (Pr_Ope_ens(Nlon,Nlat,Nlead),T2_Ope_ens(Nlon,Nlat,Nlead)) Allocate (Pr_Ope_fcs(Nlon,Nlat,Nlead),T2_Ope_fcs(Nlon,Nlat,Nlead)) Pr_Ope_ens=Miss; T2_Ope_ens=Miss; Pr_Ope_fcs=Miss; T2_Ope_fcs=Miss Allocate (Pr_Ano(Nlon,Nlat,Nlead),T2_Ano(Nlon,Nlat,Nlead)) Pr_Ano=Miss; T2_Ano=Miss Do ilead=1,Nlead imon=MOpe+ilead; if (imon>Nmon) imon=imon-Nmon ! Month of the year Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then pobs=Pr_Cli_Obs(ilon,ilat,imon); tobs=T2_Cli_Obs(ilon,ilat,imon) ! Obs islp=0 Do imod=1,Nmodel ! =2 pcli=Pr_Cli(ilon,ilat,ilead,imon,imod) ! Model Clim tcli=T2_Cli(ilon,ilat,ilead,imon,imod) Do imem=1,Nmem ! =7 Do icas=1,Ncase ! =3 islp=islp+1 pmod = Pr_Ope(ilon,ilat,imem,ilead,icas,imod) ! Ope model output tmod = T2_Ope(ilon,ilat,imem,ilead,icas,imod) xpr (islp) = Fcs (pmod,pcli,pobs) ! Fcs = (model-Cli)+Obs xtmp(islp) = Fcs (tmod,tcli,tobs) !xpr (islp)=Pr_Ope(ilon,ilat,imem,ilead,icas,imod) !xtmp(islp)=T2_Ope(ilon,ilat,imem,ilead,icas,imod) Enddo Enddo Enddo ! Fcs ensemble mean: Pr_Ope_ens(ilon,ilat,ilead)=Ave(xpr, islp,Miss) T2_Ope_ens(ilon,ilat,ilead)=Ave(xtmp,islp,Miss) !print*,Pr_Ope_ens(ilon,ilat,ilead),T2_Ope_ens(ilon,ilat,ilead) pope=Pr_Ope_ens(ilon,ilat,ilead); tope=T2_Ope_ens(ilon,ilat,ilead); pcli=0; tcli=0 Do imod=1,Nmodel ! Ens of Clim nean over all models pcli=pcli+Pr_Cli(ilon,ilat,ilead,imon,imod); tcli=tcli+T2_Cli(ilon,ilat,ilead,imon,imod); Enddo pcli=pcli/Nmodel; tcli=tcli/Nmodel !print"(10F6.1)",pobs,pcli,pope,tobs,tcli,tope ! Fcs Anomaly: if (pope/=Miss.AND.pcli/=Miss) Pr_Ano(ilon,ilat,ilead)=pope-pcli ! otherwise have been set to Miss if (tope/=Miss.AND.tcli/=Miss) T2_Ano(ilon,ilat,ilead)=tope-tcli !Pr_Ano(ilon,ilat,ilead)=Pr_Ope_ens(ilon,ilat,ilead)-Pr_Cli(ilon,ilat,ilead,imon) ! Ano = Ens_fcs - Mod_Cli !T2_Ano(ilon,ilat,ilead)=T2_Ope_ens(ilon,ilat,ilead)-T2_Cli(ilon,ilat,ilead,imon) ! Operational fcs: !if (pobs/=Miss.AND.Pr_Ano(ilon,ilat,ilead)/=Miss) Pr_Ope_fcs(ilon,ilat,ilead)=pobs+Pr_Ano(ilon,ilat,ilead) ! Fcs = Obs_Cli + Ano !if (tobs/=Miss.AND.T2_Ano(ilon,ilat,ilead)/=Miss) T2_Ope_fcs(ilon,ilat,ilead)=tobs+T2_Ano(ilon,ilat,ilead) Pr_Ope_fcs = Pr_Ope_ens; T2_Ope_fcs = T2_Ope_ens print"(10F6.1)",Pr_Ano(ilon,ilat,ilead),Pr_Ope_fcs(ilon,ilat,ilead),T2_Ano(ilon,ilat,ilead),T2_Ope_fcs(ilon,ilat,ilead) !Pr_Ope_fcs(ilon,ilat,ilead)=Pr_Cli_Obs(ilon,ilat,imon)+Pr_Ano(ilon,ilat,ilead) ! Fcs = Obs_Cli + Ano !T2_Ope_fcs(ilon,ilat,ilead)=T2_Cli_Obs(ilon,ilat,imon)+T2_Ano(ilon,ilat,ilead) ! Relative anomaly for rainfall: if (pobs>0.AND.Pr_Ano(ilon,ilat,ilead)/=Miss) Pr_Ano(ilon,ilat,ilead)=Pr_Ano(ilon,ilat,ilead)/pobs !Pr_Ano(ilon,ilat,ilead)=Pr_Ano(ilon,ilat,ilead)/Pr_Cli_Obs(ilon,ilat,imon) ! Relative Ano = Ano/Obs_Cli endif Enddo Enddo Enddo End Subroutine Do_ens_Ano_Ope !------------------------------------------------------------------------------ Real Function Fcs (Xmodel,Xcli,Xobs) Real :: Xmodel,Xcli,Xobs, Resl, Bias if (Xmodel/=Miss.AND.Xcli/=Miss.AND.Xobs/=Miss) then Bias = Xmodel - Xcli Resl = Xobs + Bias else Resl = Miss endif Fcs = Resl End Function Fcs !------------------------------------------------------------------------------ Subroutine Write_Forecast_Mean_Ano_Pr_T2 (Fname) Character (Len=*) :: Fname Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real, Allocatable, Dimension (:) :: xtime INteger :: i, Nvar1 ! Allocate (Pr_Ope_ens(Nlon,Nlat,Nlead),T2_Ope_ens(Nlon,Nlat,Nlead)) ! Pr_Ope_ens=Miss; T2_Ope_ens=Miss ! Allocate (Pr_Ano(Nlon,Nlat,Nlead),T2_Ano(Nlon,Nlat,Nlead)) Allocate (xtime(Nlead)) Do i=1,Nlead; xtime(i)=i; Enddo Write(t_unit,"(A,I4.4,A,I2.2,A)") "months since ",YOpe,"-",MOpe+1,"-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="lead"; !Fname=trim(Dir_Out)//"Forecast_Ope_Pr_T2_and_Anomaly.nc" var_ln="Model rainfall ensemble forecast" var_unit="mm/month" Call wr_ncdf_3d (Fname,Lon,Lat,xtime,Nlon,Nlat,Nlead,Pr_Ope_ens,"pr",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Model 2-m temperature ensemble forecast" var_unit="C degree" Call wr_ncdf_3d_apd (Fname,T2_Ope_ens,"t2",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Relative Anomaly of the Model rainfall forecast" var_unit="-" Call wr_ncdf_3d_apd (Fname,Pr_Ano,"pr_ano",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Anomaly of the Model 2-m temperature forecast" var_unit="C degree" Call wr_ncdf_3d_apd (Fname,T2_Ano,"t2_ano",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Pr_Ope_fcs(:,:,:), T2_Ope_fcs(:,:,:) ! Nlon,Nlat,Nlead var_ln="Model rainfall forecast obtained by adding anomaly forecast to observation climatology" var_unit="mm/month" Call wr_ncdf_3d_apd (Fname,Pr_Ope_fcs,"pr_fcs",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) var_ln="Model 2-m temperature forecast obtained by adding anomaly forecast to observation climatology" var_unit="C degree" Call wr_ncdf_3d_apd (Fname,T2_Ope_fcs,"t2_fcs",Nlon,Nlat,Nlead, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Deallocate (xtime) End Subroutine Write_Forecast_Mean_Ano_Pr_T2 !-------------------------------------------------------------------------------- Subroutine Read_Regs Integer :: ilon,ilat, ireg, ic,j Real :: tmp Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg") Nst=0; Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Nst=Nst+1 endif Enddo Enddo !print*,"Nst=",Nst; pause End Subroutine Read_Regs !------------------------------------------------------------------------------ Subroutine Read_Clim_Data !(Fname) !Character (Len=*) :: Fname Integer :: Nlo,Nla,Nle,Nt,imod,it Real, Allocatable :: xlon(:),xlat(:),xled(:),xtime(:),X4d(:,:,:,:) ! Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:) ! Nlon,Nlat,Nmon Clim_Mean_Bias_RegCM_pr.nc ! Pr_Cli(:,:,:,:,:), T2_Cli(:,:,:,:,:) ! Nlon,Nlat,Nlead,Nmon,Nmodel ! Pr_Cli_Obs(:,:,:), T2_Cli_Obs(:,:,:) ! Nlon,Nlat,Nmon Allocate (Pr_Cli(Nlon,Nlat,Nlead,Nmon,Nmodel), T2_Cli(Nlon,Nlat,Nlead,Nmon,Nmodel)) Allocate (Pr_Cli_Obs(Nlon,Nlat,Nmon), T2_Cli_Obs(Nlon,Nlat,Nmon)) Do imod=1,Nmodel ! Cac mang nay tinh tu BC_ANN_1982-2020_Feb_2024.f90 Fname=trim(Cli_Dir)//"Clim_Mean_Bias_"//trim(Model_nm(imod))//"_pr.nc" ! Chi can doc tu 1 file *_pr.nc hoac *_t2.nc Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,X4d,"pr") !_m") Pr_Cli(:,:,:,:,imod)=X4d(:,:,:,:) Deallocate (xlon,xlat,xled,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,X4d,"t2") !_m") T2_Cli(:,:,:,:,imod)=X4d(:,:,:,:) Deallocate (xlon,xlat,xled,xtime,X4d) if (imod==1) then Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,X4d,"pr_o") Do it=1,Nt Pr_Cli_Obs(:,:,it)=X4d(:,:,imod,it) Enddo !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,Pr_Cli_Obs,"pr_o") !print*,Pr_Cli_Obs; pause Deallocate (xlon,xlat,xled,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xled,xtime,Nlo,Nla,Nle,Nt,X4d,"t2_o") Do it=1,Nt T2_Cli_Obs(:,:,it)=X4d(:,:,imod,it) Enddo !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,T2_Cli_Obs,"t2_o") !Deallocate (xlon,xlat,xtime) Deallocate (xlon,xlat,xled,xtime,X4d) endif Enddo End Subroutine Read_Clim_Data !------------------------------------------------------------------------------ Subroutine Read_Ope_Data Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:),xtime(:) ! Pr_Ope(:,:,:,:,:,:), T2_Ope(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel Allocate (T2_Ope(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel)) Allocate (Pr_Ope(Nlon,Nlat,Nmem,Nlead,Ncase,Nmodel)) Do imod=1,Nmodel write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"BC_Ope_",YOpe,MOpe,"_"//trim(TVar)//"_"//trim(Model_nm(imod))//".nc" Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))) T2_Ope(:,:,:,:,icas,imod)=X4d(:,:,:,:) !print*,T2_Ope(:,:,:,:,icas,imod) Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo Enddo Do imod=1,Nmodel write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"BC_Ope_",YOpe,MOpe,"_"//trim(RVar)//"_"//trim(Model_nm(imod))//".nc" Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))) Pr_Ope(:,:,:,:,icas,imod)=X4d(:,:,:,:) Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo Enddo End Subroutine Read_Ope_Data !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" !-------------------------------------------------------------------- END ! Program