Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99.0 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_All(:,:,:,:,:,:,:), T2_All(:,:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel Real, Allocatable :: Pr_Obs(:,:,:), T2_Obs(:,:,:) ! Nlon,Nlat,Ntime Real, Allocatable :: Pr_Obs_Cli(:,:,:), T2_Obs_Cli(:,:,:) ! Nlon,Nlat,Nmon Real, Allocatable :: Pr_ens(:,:,:,:), T2_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime Real, Allocatable :: Pr_ens_Cli(:,:,:,:), T2_ens_Cli(:,:,:,:) ! Nlon,Nlat,Nlead,Nmon Real, Allocatable :: xmem(:), xlead(:), xtime(:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File Character (Len=3) :: TVar, RVar Character (Len=10),Allocatable :: Model_nm(:), IDXname(:) Logical :: status, Empty Integer :: Yr1, Yr2, Nmodel, NYear, maxyrs, Nidx 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, xspi !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Nmodel Namelist /Model_name/ Model_nm !------------------------------------------------------------------------------ Open (1, file="Namelist_Pr_T2_Clim.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) Close (1) Ntime = (Yr2-Yr1+1)*Nmon; NYear=Yr2-Yr1+1 ! Chi tinh cho giai doan 1983-2020 Call Read_Regs print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime Call Read_all_data ! Read all data of "cli","ann","mlr" 1983 to present Call Calc_ens_model ! Ensemble over Nspl=Nmodel*Ncase*Nmem corrected members. Currently, Nspl=2*3*7=42 Call Calc_Obs_Monthly_Cli ! Observed Climatological means Call Calc_Model_Monthly_Cli ! Model Ensemble Climatological means Fname=trim(Out_Dir)//"Monthly_Clim.nc" Call Write_Pr_Temp_Cli (Fname) CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_Model_Monthly_Cli Integer :: ilon,ilat,it,imon,imon1,Nyr,iyr Real :: xtmp(Nyear),xpr(Nyear) ! Pr_ens_Cli(:,:,:,:), T2_ens_Cli(:,:,:,:) ! Nlon,Nlat,Nlead,Nmon ! Pr_ens(Nlon,Nlat,Nlead,Ntime),T2_ens(Nlon,Nlat,Nlead,Ntime) Allocate (Pr_ens_Cli(Nlon,Nlat,Nlead,Nmon),T2_ens_Cli(Nlon,Nlat,Nlead,Nmon)) Pr_ens_Cli=Miss; T2_ens_Cli=Miss Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do imon=1,Nmon iyr=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then iyr=iyr+1; xtmp(iyr)=T2_ens(ilon,ilat,ilead,it); xpr (iyr)=Pr_ens(ilon,ilat,ilead,it) endif Enddo T2_ens_Cli(ilon,ilat,ilead,imon)=Ave(xtmp,iyr,Miss) Pr_ens_Cli(ilon,ilat,ilead,imon)=Ave(xpr, iyr,Miss) Enddo endif Enddo Enddo Enddo End Subroutine Calc_Model_Monthly_Cli !------------------------------------------------------------------------------ Subroutine Calc_Obs_Monthly_Cli Integer :: ilon,ilat,it,imon,imon1,Nyr,iyr Real :: xtmp(Nyear),xpr(Nyear) ! Pr_Obs_Cli(:,:,:), T2_Obs_Cli(:,:,:) ! Nlon,Nlat,Nmon Allocate (Pr_Obs_Cli(Nlon,Nlat,Nmon),T2_Obs_Cli(Nlon,Nlat,Nmon)) Pr_Obs_Cli=Miss; T2_Obs_Cli=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do imon=1,Nmon iyr=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then iyr=iyr+1; xtmp(iyr)=T2_Obs(ilon,ilat,it); xpr(iyr)=Pr_Obs(ilon,ilat,it) endif Enddo T2_Obs_Cli(ilon,ilat,imon)=Ave(xtmp,iyr,Miss) Pr_Obs_Cli(ilon,ilat,imon)=Ave(xpr, iyr,Miss) Enddo endif Enddo Enddo End Subroutine Calc_Obs_Monthly_Cli !------------------------------------------------------------------------------ Subroutine Write_Pr_Temp_Cli (Fname) Character (Len=*) :: Fname Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg Real :: xlead(Nlead),xmon(Nmon) Real, Allocatable :: X4d(:,:,:,:),xcase(:) ! Real, Allocatable :: Pr_Obs_Cli(:,:,:), T2_Obs_Cli(:,:,:) ! Nlon,Nlat,Nmon ! Real, Allocatable :: Pr_ens_Cli(:,:,:,:), T2_ens_Cli(:,:,:,:) ! Nlon,Nlat,Nlead,Nmon !write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Monthly_Clim.nc" t_unit="months since 1983-01-01 " lon_nm="lon"; lat_nm="lat"; z_nm="lead"; t_nm="month" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do imon=1,Nmon; xmon(imon)=imon; Enddo Var_ln="Climatological Monthly Ensemble Mean of Model Rainfall averaged over 7 members"// & " of Cli and ANN bias correction of RegCM and clWRF for 6 leadtimes (z)" var_unit="mm/month" Call wr_ncdf_4d (Fname,Lon,Lat,xlead,xmon, Nlon,Nlat,Nlead,Nmon,Pr_ens_Cli,"pr_m", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Climatological Monthly Ensemble Mean of Model 2m-Temperature averaged over 7 members"// & " of Cli and ANN bias correction of RegCM and clWRF for 6 leadtimes (z)" var_unit="C degree" Call wr_ncdf_4d_apd (Fname,T2_ens_Cli,"t2_m",Nlon,Nlat,Nlead,Nmon, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Climatological Monthly Ensemble Mean of Observed Rainfall" var_unit="mm/month" Call wr_ncdf_3d_apd (Fname,Pr_Obs_Cli,"pr_o",Nlon,Nlat,Nmon,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Var_ln="Climatological Monthly Ensemble Mean of Observed 2m-Temperature" var_unit="C degree" Call wr_ncdf_3d_apd (Fname,T2_Obs_Cli,"t2_o",Nlon,Nlat,Nmon,Var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) End Subroutine Write_Pr_Temp_Cli !------------------------------------------------------------------------------ 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 !-------------------------------------------------------------------------------- Integer Function Get_Reg_Points (ireg) Integer :: ireg,ilon,ilat,Npoint Npoint=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==ireg) Npoint=Npoint+1 Enddo Enddo Get_Reg_Points = Npoint End Function Get_Reg_Points !------------------------------------------------------------------------------ Subroutine Calc_ens_model Integer :: ispl,imod,ilon,ilat,ilead,icas,imem, Nspl Real, Allocatable :: xtmp(:),xpr(:) Nspl=Nmodel*Ncase*Nmem Allocate (xtmp(Nspl),xpr(Nspl)) Allocate (Pr_ens(Nlon,Nlat,Nlead,Ntime),T2_ens(Nlon,Nlat,Nlead,Ntime)) ! Pr_ens(:,:,:,:), T2_ens(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime Pr_ens=Miss; T2_ens=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do ilead=1,Nlead Do it=1,Ntime ispl=0 ! = Nmodel x Nmem x Ncase = 2 x 7 x 3 = 42 Do imod=1,Nmodel Do imem=1,Nmem Do icas=1,Ncase ispl=ispl+1; xpr (ispl)=Pr_All(ilon,ilat,imem,ilead,it,icas,imod) xtmp(ispl)=T2_All(ilon,ilat,imem,ilead,it,icas,imod) Enddo Enddo Enddo Pr_ens(ilon,ilat,ilead,it)=Ave(xpr, ispl,Miss) T2_ens(ilon,ilat,ilead,it)=Ave(xtmp,ispl,Miss) Enddo Enddo endif Enddo Enddo End Subroutine Calc_ens_model !------------------------------------------------------------------------------ Subroutine Read_all_data ! (Model_nm) !Character (Len=*) :: Model_nm Integer :: imod,ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2, Yr0 Character (Len=255) :: Fname, Dir_obs Real, Allocatable :: X3d(:,:,:), X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:) !,xtime(:) it1=Nmon+1 ! Tinh tu 1/1983 pr_clWRF_lead_01.nc Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)) Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)) T2_All=Miss; Pr_All=Miss; if (Allocated(X4d)) Deallocate (X4d); if (Allocated(xmem)) Deallocate (xmem) if (Allocated(xlon)) Deallocate (xlon); if (Allocated(xlat)) Deallocate (xlat) if (Allocated(time)) Deallocate (time); if (Allocated(xtime)) Deallocate (xtime) Do imod=1,Nmodel Do ilead=1,Nlead write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" print*,"Reading ",trim(Fname) Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))) Do imem=1,Nmem Do it=1,Ntime T2_All(:,:,imem,ilead,it,icas,imod)=X4d(:,:,imem,Nmon+it) ! Luu vao T2_All voi time tu Nmon+1 = it1 (Tinh tu 1/1983) Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(RVar)//"_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))) print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt !if (imod==1.AND.icas==1.AND.ilead==1) then; Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel)); Pr_All=Miss; endif Do imem=1,Nmem Do it=1,Ntime Pr_All(:,:,imem,ilead,it,icas,imod)=X4d(:,:,imem,Nmon+it) ! Luu vao Pr_All voi time tu Nmon+1 = it1 (Tinh tu 1/1983) Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo Print*,"Reading data ...LT=",ilead Enddo Enddo Print*,"Nt=",Nt ! Read OBS data: Dir_obs="/home/tanpv/Data/OBS_Vn_GRD/0.25/Monthly/" ! Pr_Obs(:,:,:), T2_Obs(:,:,:) ! Nlon,Nlat,Ntime Allocate (Pr_Obs(Nlon,Nlat,Ntime),T2_Obs(Nlon,Nlat,Ntime)) Pr_Obs=Miss; T2_Obs=Miss Fname=trim(Dir_obs)//"R.mon.1980-2021_0.25.nc" Yr0=1980 it1=(Yr1-Yr0)*Nmon+1; it2=it1+Ntime-1 ! it1=1/1983, it2=12/2020 Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"r") Pr_Obs(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) Fname=trim(Dir_obs)//"T2m.mon.1961-2021_0.25.nc" Yr0=1961; it1=(Yr1-Yr0)*Nmon+1; it2=it1+Ntime-1 Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"t2m") T2_Obs(:,:,1:Ntime)=X3d(:,:,it1:it2) Deallocate (xlon,xlat,xtime,X3d) Print*,"Reading data ... Done" End Subroutine Read_all_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" !Include "/home/tanpv/Mylibs/Sub_SPI_One_Stn_2023.f90" !Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_2023.f90" !Include "/home/tanpv/Mylibs/PET_thornthwaite.f90" !Include "/home/tanpv/Mylibs/Sub_SPEI_One_Stn_PET_2023.f90" !Include "/home/tanpv/Mylibs/Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI !Include "/home/tanpv/Mylibs/Sub_SPEI.f90" ! Chuong trinh tinh chi so han SPEI Include "/home/tanpv/Mylibs/Sub_Percentile.f90" !-------------------------------------------------------------------- END ! Program