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 = 4 ! mod,cli,ann,mlr Character (Len=3), Parameter :: Cases(Ncase)=(/"mod","cli","ann","mlr"/) Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Pr_All(:,:,:,:,:,:), T2_All(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase Real, Allocatable :: Pr_1 (:,:,:,:,:), T2_1(:,:,:,:,:) ! Nlon,Nlat,Nmem,Ntime,Ncase - Lay LT=1 noi voi cac LT khac Real, Allocatable :: Pr_Ope(:,:,:,:,:), T2_Ope(:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ncase Real, Allocatable :: xmem(:), xlead(:), xtime(:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname Character (Len=255) :: Inp_Dir, Dir_Ope, Mask_File 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/ Inp_Dir, Model_nm, Mask_File, Dir_Ope, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Nmodel Namelist /Model_name/ Model_nm Namelist /Time_Infor/ YOpe, MOpe !------------------------------------------------------------------------------ Open (1, file="Namelist_Combine_for_Ope.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) ! Yr1=1982, Yr2=2025 Allocate (Model_nm(Nmodel)) Read (1, Model_name) Read (1, Time_Infor) Close (1) Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*Nmon Call Read_Regs print*,"Parameters: Nlat,Nlon,Ntime=",Nlat,Nlon,Ntime Do imod=1,Nmodel Call Reinit Call Read_last_data (Model_nm(imod)) Call Read_Ope_Data (Model_nm(imod)) print*,"Reading data done ..." !------------------------------------------------------------------------------ print*,"Parameters: Nlat,Nlon,Ntime=",Nlat,Nlon,Ntime ! Ntime=Ntime1 ! write (Fname,"(A,I1,A)") trim(Out_Dir)//"SPEI_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" ! Call Write_Dr_IDXs_to_file (Fname,ilead,Model_nm(imod)) ! Enddo Call Combine Call Write_for_drought (Model_nm(imod)) Call Write_for_drought_Lead_1 (Model_nm(imod)) Enddo CONTAINS !------------------------------------------------------------------------------ Subroutine Reinit if (Allocated(T2_All)) Deallocate (T2_All) if (Allocated(Pr_All)) Deallocate (Pr_All) if (Allocated(T2_1)) Deallocate (T2_1) if (Allocated(Pr_1)) Deallocate (Pr_1) if (Allocated(Pr_Ope)) Deallocate (Pr_Ope) if (Allocated(T2_Ope)) Deallocate (T2_Ope) Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase)) Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase)) Allocate (Pr_1(Nlon,Nlat,Nmem,Ntime,Ncase)) Allocate (T2_1(Nlon,Nlat,Nmem,Ntime,Ncase)) ! Nlon,Nlat,Nmem,Ntime,Ncase Allocate (Pr_Ope(Nlon,Nlat,Nmem,Nlead,Ncase)) Allocate (T2_Ope(Nlon,Nlat,Nmem,Nlead,Ncase)) ! Nlon,Nlat,Nmem,Nlead,Ncase T2_All=Miss; Pr_All=Miss; Pr_1=Miss; T2_1=Miss; Pr_Ope=Miss; T2_Ope=Miss End Subroutine Reinit !-------------------------------------------------------------------------------- Subroutine Write_for_drought (Model_nm) Character (Len=*) :: Model_nm Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real, Allocatable, Dimension (:) :: zlev, xtime INteger :: i, Nvar1 Allocate (zlev(Nmem),xtime(Ntime)) Do i=1,Nmem; zlev(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo t_unit="months since 1982-01-01 "; var_unit="-" lon_nm="lon"; lat_nm="lat"; t_nm="month"; z_nm="mem" ! Pr_All(:,:,imem,ilead,it,icas) Allocate (X4d(Nlon,Nlat,Nmem,Ntime)) Do ilead=1,NLead ! T2 write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" !write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" Do ivar=1,Ncase Do imem=1,Nmem Do it=1,Ntime X4d(:,:,imem,it)=T2_All(:,:,imem,ilead,it,ivar) Enddo Enddo var_ln="Monthly mean values of "//trim(TVar)//". z=1..6 member, z=7 Ens. t=1-1982..12-2025" var_nm=trim(Cases(ivar)) if (ivar==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nmem,Ntime,X4d,var_nm, & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nlon,Nlat,Nmem,Ntime, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo Enddo ! Do ilead=1,NLead ! Pr write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(RVar)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" !write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" Do ivar=1,Ncase Do imem=1,Nmem Do it=1,Ntime X4d(:,:,imem,it)=Pr_All(:,:,imem,ilead,it,ivar) Enddo Enddo var_ln="Monthly mean values of "//trim(RVar)//". z=1..6 member, z=7 Ens. t=1-1982..12-2025" var_nm=trim(Cases(ivar)) if (ivar==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nmem,Ntime,X4d,var_nm, & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nlon,Nlat,Nmem,Ntime, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo Enddo Deallocate (zlev,xtime,X4d) End Subroutine Write_for_drought !-------------------------------------------------------------------------------- Subroutine Write_for_drought_Lead_1 (Model_nm) Character (Len=*) :: Model_nm Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real, Allocatable, Dimension (:) :: zlev, xtime INteger :: i, Nvar1 Allocate (zlev(Nmem),xtime(Ntime)) Do i=1,Nmem; zlev(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo t_unit="months since 1982-01-01 "; var_unit="-" lon_nm="lon"; lat_nm="lat"; t_nm="month"; z_nm="mem" ! Pr_All(:,:,imem,ilead,it,icas) Allocate (X4d(Nlon,Nlat,Nmem,Ntime)) !Do ilead=1,NLead ! T2 write (Fname,"(A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm)//".nc" !write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" Do ivar=1,Ncase Do imem=1,Nmem Do it=1,Ntime X4d(:,:,imem,it)=T2_1(:,:,imem,it,ivar) Enddo Enddo var_ln="Monthly mean values of "//trim(TVar)//". z=1..6 member, z=7 Ens. t=1-1982..12-2025" var_nm=trim(Cases(ivar)) if (ivar==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nmem,Ntime,X4d,var_nm, & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nlon,Nlat,Nmem,Ntime, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo !Enddo ! !Do ilead=1,NLead ! Pr write (Fname,"(A)") trim(Inp_Dir)//trim(RVar)//"_"//trim(Model_nm)//".nc" !write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" Do ivar=1,Ncase Do imem=1,Nmem Do it=1,Ntime X4d(:,:,imem,it)=Pr_1(:,:,imem,it,ivar) Enddo Enddo var_ln="Monthly mean values of "//trim(RVar)//". z=1..6 member, z=7 Ens. t=1-1982..12-2025" var_nm=trim(Cases(ivar)) if (ivar==1) then Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nmem,Ntime,X4d,var_nm, & t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss) else Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nlon,Nlat,Nmem,Ntime, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo !Enddo Deallocate (zlev,xtime,X4d) End Subroutine Write_for_drought_Lead_1 !------------------------------------------------------------------------------ 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 Combine Integer :: ilon,ilat,ilead,imem,it,icas, it1, it2, Nt0, Nt1 ! Pr_All(:,:,imem,ilead,it,icas) Pr_Ope(Nlon,Nlat,Nmem,Nlead,Ncase) Nt0=(YOpe-Yr1)*Nmon ! (2021-1982)*12 = 468 = 12/2020 ung voi Initial=12/2020 Nt1=Nt0+MOpe ! Operational Initial month Do ilead=1,Nlead it1=Nt1+ilead if (ilead==1) then Pr_1(:,:,:,:,:)=Pr_All(:,:,:,ilead,:,:) T2_1(:,:,:,:,:)=T2_All(:,:,:,ilead,:,:) endif Do ilon=1,Nlon Do ilat=1,Nlat Do imem=1,Nmem Do icas=1,Ncase Pr_All(ilon,ilat,imem,ilead,it1,icas)=Pr_Ope(ilon,ilat,imem,ilead,icas) T2_All(ilon,ilat,imem,ilead,it1,icas)=T2_Ope(ilon,ilat,imem,ilead,icas) Pr_1 (ilon,ilat,imem,Nt1+ilead,icas)=Pr_Ope(ilon,ilat,imem,ilead,icas) T2_1 (ilon,ilat,imem,Nt1+ilead,icas)=T2_Ope(ilon,ilat,imem,ilead,icas) Enddo Enddo Enddo Enddo Enddo End Subroutine Combine !------------------------------------------------------------------------------ Subroutine Read_Ope_Data (Model_nm) Character (Len=*) :: Model_nm Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:),xtime(:) write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"BC_Ope_",YOpe,MOpe,"_"//trim(TVar)//"_"//trim(Model_nm)//".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)=X4d(:,:,:,:) Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo write (Fname,"(A,I4.4,I2.2,A)") trim(Dir_Ope)//"BC_Ope_",YOpe,MOpe,"_"//trim(RVar)//"_"//trim(Model_nm)//".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)=X4d(:,:,:,:) Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo End Subroutine Read_Ope_Data !------------------------------------------------------------------------------ Subroutine Read_last_data (Model_nm) Character (Len=*) :: Model_nm Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:),xtime(:) !it1=Nmon+1 ! Tinh tu 1/1983 - pr_RegCM_lead_01.nc !it2=it1+Ntime-1 ! Tinh den 12/2021 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 ilead=1,Nlead write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//trim(TVar)//"_"//trim(Model_nm)//"_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))) print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt !if (icas==1.AND.ilead==1) then ! Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase),time(Ntime)); ! T2_All=Miss; time(1:Ntime)=xtime(1:Ntime) ! Luu lai de write out !endif Do imem=1,Nmem Do it=1,Nt ! Ntime T2_All(:,:,imem,ilead,it,icas)=X4d(:,:,imem,it) !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)//"_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 (icas==1.AND.ilead==1) then; Allocate (Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase)); Pr_All=Miss; endif Do imem=1,Nmem Do it=1,Nt !Ntime Pr_All(:,:,imem,ilead,it,icas)=X4d(:,:,imem,it) !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 Print*,"Nt=",Nt if (Allocated(xtime)) Deallocate (xtime) Allocate (xtime(Ntime)) !pause;print*,time !xtime(:)=time(:) ! restore xtime !Deallocate (time) Print*,"Reading data ... Done" End Subroutine Read_last_data !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" !-------------------------------------------------------------------- END ! Program