Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6 Integer, Parameter :: Ncase = 4 ! obs,mod,cli,ann Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","mod","cli","ann"/) !Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","cli","ann"/) Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI Real, Allocatable :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:) Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Pr_All(:,:,:,:,:,:), T2_All(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Nlead,Ntime,Ncase Real, Allocatable :: xmem(:), xlead(:), xtime(:) Real, Allocatable :: SPEI_All(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Ntime,Ncase,Nlen ! Notes: ! 1) Dat_Mon(ivar,ist,iyr,imw), where ivar=1:Temp (F), ivar=2:Rainfall (inch); ist=1..Nst; iyr=YBas1..YBas2; imw=1..Nmon (or Nweek) ! Dat_Mon used to calculating PDSI ! 2) Temp_Inp (Temp in C deg), Rain_Inp (Rainfall in mm/month) (1:Nlon1,1:Nlat1,1:Ntime), Ntime=(YBas2-YBas1+1)*Nmon ! Temp_Inp and Rain_Inp to be used to write out as input data ! 3) Prec (mm/month), Temp (C) (1:Ntime,1:Nst) ! 4) DatBL (mm/month), DatTMP (C) (1:Ntime,1:Nst): Used to calculating SPI, SPEI !------------------------------------------------------------------------------ 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(:) Logical :: status, Empty Integer :: Yr1, Yr2, Yr_fc, Mn_fc, 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, Out_Dir, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Yr_fc, Mn_fc, Nmodel, NLen Namelist /Model_name/ Model_nm Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_Dr_Obs_Mod.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Len_SPI(NLen)); Allocate (Model_nm(Nmodel)) Read (1, Model_name) Read (1, Dr_Timescale) Close (1) Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Nt = Ntime - Nmon ! Tru SL nam cuoi cung (vi chi co cac thang du bao ung voi cac LT) Call Read_Regs print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime Do imod=1,Nmodel Call Reinit Call Read_all_data (Model_nm(imod)) Call Data_Preparation print*,"Reading data done ..." !------------------------------------------------------------------------------ print*,"Beginning: Ntime, Nst=",Ntime, Nst print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime ! For calculating SPI, SPEI: Ntime1=Ntime !!!!!!!!!! Do ilead=1,Nlead Nt1 = Nt + ilead SPEI_All=Miss Do imem=1,Nmem Do icas=1,Ncase ist=0; DatBL=Miss; DatTMP=Miss; Ntime=Ntime1 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then ist=ist+1; DatBL (1:Ntime,ist) = Pr_All(ilon,ilat,imem,ilead,1:Ntime,icas) !Rain_Inp(ilon,ilat,1:Ntime) DatTMP(1:Ntime,ist) = T2_All(ilon,ilat,imem,ilead,1:Ntime,icas) !Temp_Inp(ilon,ilat,1:Ntime) Lat(ist)=xlat(ilat) !X1d(1:Ntime)=DatBL (1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1 !X1d(1:Ntime)=DatTMP(1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1 !X1d(1:Ntime)=DatTMP(1:Ntime,ist); Nt1=Length_series (X1d,Ntime,Miss) !if (Nt1>Ntime1/3) then; Ntime=Nt1; Nt=Nt1; endif endif Enddo Enddo !pause Ntime=Nt1 if (Allocated(SPI)) Deallocate (SPI); Allocate (SPI(Nst,Ntime,NLen)) SPI(:,:,:) = amssng Call Cal_SPEI !Call Cal_SPI !if (icas==1.AND.imem==3) then; print*,SPI,"LT=",ilead; pause; endif ist=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then ist=ist+1 Do ilen=1,Nlen SPEI_All(ilon,ilat,imem,1:Ntime,icas,ilen)=SPI(ist,1:Ntime,ilen) Enddo !X1d(1:Ntime)=SPI(ist,1:Ntime,1); Nt1=Length_series (X1d,Ntime,Miss); print"(A,8I6)","ilead,imem,icas: ",ilead,imem,icas,Nt1 endif Enddo Enddo !pause Enddo Enddo !iyr=Nt1/Nmon; imon=Nt1-iyr*Nmon; iyr=Yr1+iyr-1; !print*,"Years: ",Yr1,"-",iyr,"/",imon; pause Ntime=Ntime1 write (Fname,"(A,I1,A)") trim(Out_Dir)//"SPEI_Obs_Fcs_"//trim(Model_nm(imod))//"_lead_",ilead,".nc" Call Write_Dr_IDXs_to_file (Fname,ilead,Model_nm(imod)) Enddo Enddo CONTAINS !------------------------------------------------------------------------------ Subroutine Reinit if (Allocated(T2_All)) Deallocate (T2_All) if (Allocated(Pr_All)) Deallocate (Pr_All) if (Allocated (Rain)) Deallocate (Rain) if (Allocated (DatBL)) Deallocate (DatBL) if (Allocated (DatTMP)) Deallocate (DatTMP) if (Allocated (beta_spi)) Deallocate (beta_spi) if (Allocated (gamm)) Deallocate (gamm) if (Allocated (pzero)) Deallocate (pzero) if (Allocated (SPEI_All)) Deallocate (SPEI_All) if (Allocated (Lat)) Deallocate (Lat) if (Allocated (X1d)) Deallocate (X1d) Allocate (Rain(Ntime),DatBL(Ntime,Nst),DatTMP(Ntime,Nst)) Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth)) !Allocate (SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen)) Allocate (SPEI_All(Nlon,Nlat,Nmem,Ntime,Ncase,Nlen)) Allocate (Lat(Nst),X1d(Ntime)) End Subroutine Reinit !------------------------------------------------------------------------------ Subroutine Write_Dr_IDXs_to_file (Fname,ilead,Model_nm) Character (Len=*) :: Fname,Model_nm Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: it, imem, ilead, ireg, icas, ilen, imod Real :: xmem(Nmem) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="mem"; t_nm="mon" Do imem=1,Nmem; xmem(imem)=imem; Enddo Allocate (X4d(Nlon,Nlat,Nmem,Ntime)) Do icas=1,Ncase Do ilen=1,Nlen Do imem=1,Nmem Do ilon=1,Nlon Do ilat=1,Nlat Do it=1,Ntime X4d(ilon,ilat,imem,it)=SPEI_All(ilon,ilat,imem,it,icas,ilen) Enddo Enddo Enddo Enddo write(Var_ln,"(A,I2.2,A,I1)") "SPEI for "//trim(Cases(icas))//" timescale = ",ilen," months of model "//trim(Model_nm)//" with Leadtime=",ilead write(Var_nm,"(A,I2.2)") trim(Cases(icas))//"_",Len_SPI(ilen) if (icas==1.AND.ilen==1) then Call wr_ncdf_4d (Fname, xlon,xlat,xmem,xtime, Nlon,Nlat,Nmem,Ntime,X4d,Var_nm, & t_unit,Var_ln,var_unit,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 (X4d) End Subroutine Write_Dr_IDXs_to_file !------------------------------------------------------------------------------ Integer Function Length_series (X,Ninp,Miss) Integer :: Ninp,i,j Real :: X(Ninp), Miss j=0 Do i=1,Ninp !,1,-1 if (X(i)/=Miss) then j=j+1 else exit endif Enddo !print*,Ninp,j Length_series=j !Ninp-j End Function Length_series !------------------------------------------------------------------------------ Subroutine Read_Regs Integer :: ilon,ilat, ireg, ic,j Real :: tmp Call rd_ncdf_2d (Mask_File,xlon,xlat,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 Data_Preparation Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Integer :: Nt1, Nt2, Nyr1, Nyr2 Nyr1 = (Yr_fc-Yr1) ! So nam su dung SL OBS tinh tu Yr1 = (2011-1983) = 28 yrs Nt1 = Nyr1*Nmon + Mn_fc - 1 ! So thang su dung SL OBS cho LT=1 tinh tu 01/Yr1 = 28 * 12 + 3 = 339 Do ilead=1,Nlead it2=Nt1+ilead-1 ! So thu tu thang tinh tu 01/Yr1 (voi LT=1: it2=339) Do imem=1,Nmem Do ilon=1,Nlon Do ilat=1,Nlat Do icas=2,Ncase T2_All(ilon,ilat,imem,ilead,1:it2,icas)=T2_All(ilon,ilat,imem,ilead,1:it2,1) ! Gan SL du bao boi SL Obs Pr_All(ilon,ilat,imem,ilead,1:it2,icas)=Pr_All(ilon,ilat,imem,ilead,1:it2,1) Enddo Enddo Enddo Enddo Enddo End Subroutine Data_Preparation !------------------------------------------------------------------------------ Subroutine Read_all_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 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,Ntime T2_All(:,:,imem,ilead,it,icas)=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)//"_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,Ntime Pr_All(:,:,imem,ilead,it,icas)=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 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_all_data !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.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 !-------------------------------------------------------------------- END ! Program