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 Integer, Parameter :: Ncase_1 = 2 ! obs, ens Character (Len=3), Parameter :: Case_1(Ncase_1)=(/"obs","ens"/) 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,Nmodel Real, Allocatable :: xmem(:), xlead(:), xtime(:) Real, Allocatable :: SPEI_All(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Ntime,Ncase,Nlen !------------------------------------------------------------------------------ ! Ensemble: Real, Allocatable :: Pr_ens (:,:,:,:,:), T2_ens(:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_1 Real, Allocatable :: SPEI_ens(:,:,:,:,:) ! Nlon,Nlat,Ntime,Ncase_1,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, 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, Nmodel, NLen Namelist /Model_name/ Model_nm Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_Dr_Reg.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 Call Read_all_data Call Do_Ensemble !stop ! T2_ens(:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_1 !Do imod=1,Nmodel Call Reinit !Call Read_all_data (Model_nm(imod)) 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: SPEI_ens=Miss Ntime1=Ntime !!!!!!!!!! Do ilead=1,Nlead Nt1 = Nt + ilead Do icas=1,Ncase_1 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,imod) !Rain_Inp(ilon,ilat,1:Ntime) !DatTMP(1:Ntime,ist) = T2_All(ilon,ilat,imem,ilead,1:Ntime,icas,imod) !Temp_Inp(ilon,ilat,1:Ntime) DatBL (1:Ntime,ist) = Pr_ens(ilon,ilat,ilead,1:Ntime,icas) DatTMP(1:Ntime,ist) = T2_ens(ilon,ilat,ilead,1:Ntime,icas) Lat(ist)=xlat(ilat) 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_ens(ilon,ilat,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 Ntime=Ntime1 write (Fname,"(A,I1,A)") trim(Out_Dir)//"SPEI_ens_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_ens)) Deallocate (SPEI_ens) 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 (SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen)) ! Nlon,Nlat,Ntime,Ncase_1,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 :: zlev(Nlen) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="timescale"; t_nm="time" Do ilen=1,Nlen; zlev(ilen)=Len_SPI(ilen); Enddo Allocate (X4d(Nlon,Nlat,Nlen,Ntime)) ! SPEI_ens(Nlon,Nlat,Ntime,Ncase_1,Nlen) Do icas=1,Ncase_1 Do ilen=1,Nlen Do ilon=1,Nlon Do ilat=1,Nlat Do it=1,Ntime X4d(ilon,ilat,ilen,it)=SPEI_ens(ilon,ilat,it,icas,ilen) Enddo Enddo Enddo Enddo write(Var_ln,"(A,I2.2,A,I1)") "SPEI for "//trim(Case_1(icas))//" with timescale in z=1,2,3 corresponding to 1,3,6 month(s) for Leadtime=",ilead write(Var_nm,"(A,I2.2)") trim(Case_1(icas)) !//"_",Len_SPI(ilen) if (icas==1) then Call wr_ncdf_4d (Fname, xlon,xlat,zlev,xtime, Nlon,Nlat,Nlen,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,Nlen,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) endif 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 Read_all_data ! (Model_nm) !Character (Len=*) :: Model_nm Integer :: imod,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 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))) print*,"Nlo,Nla,Nmem,Ntime=",Nlo,Nla,Nme,Nt if (imod==1.AND.icas==1.AND.ilead==1) then Allocate (T2_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel),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,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 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 !------------------------------------------------------------------------------ Subroutine Do_Ensemble Integer :: ilon,ilat,it,ilead,imem,icas, i,j,k, N Real, Allocatable :: X1d(:),Y1d(:) Real :: temp,rain ! Pr_ens(:,:,:,:,:), T2_ens(:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase_1 ! Pr_All(Nlon,Nlat,Nmem,Nlead,Ntime,Ncase,Nmodel) N = Nmodel*Nmem*(Ncase-1) Allocate (X1d(N),Y1d(N)) Allocate (T2_ens(Nlon,Nlat,Nlead,Ntime,Ncase_1),Pr_ens(Nlon,Nlat,Nlead,Ntime,Ncase_1)) T2_ens=Miss; Pr_ens=Miss Do ilead=1,Nlead Do it=1,Ntime T2_ens(:,:,ilead,it,1)=T2_All(:,:,1,ilead,it,1,1) Pr_ens(:,:,ilead,it,1)=Pr_All(:,:,1,ilead,it,1,1) Enddo Enddo Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Do ilead=1,Nlead Do it=1,Ntime i=0 Do imod=1,Nmodel Do imem=1,Nmem Do icas=3,Ncase ! "mod","cli","ann" i=i+1 X1d(i)=T2_All(ilon,ilat,imem,ilead,it,icas,imod) Y1d(i)=Pr_All(ilon,ilat,imem,ilead,it,icas,imod) !print*,i,X1d(i),Y1d(i) Enddo Enddo Enddo temp=percentile (X1d,i,0.5,Miss); rain=percentile (Y1d,i,0.5,Miss) T2_ens(ilon,ilat,ilead,it,2)=temp Pr_ens(ilon,ilat,ilead,it,2)=rain !print*,i,temp,rain Enddo Enddo endif Enddo Enddo End Subroutine Do_Ensemble !------------------------------------------------------------------------------ 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 Include "/home/tanpv/Mylibs/Sub_Percentile.f90" !-------------------------------------------------------------------- END ! Program