Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2, Nreg = 7, Nmem = 7, Nlead = 6 Integer, Parameter :: NMonth = 12, Nmon = NMonth Integer, Parameter :: Ncase = 4 ! obs,mod,cli,ann Character (Len=3), Parameter :: Cases(Ncase)=(/"obs","mod","cli","ann"/) Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI !Real, Allocatable :: PET(:,:) Real, Allocatable :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:) !Real, Allocatable, Dimension (:,:) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind !Real, Allocatable, Dimension (:,:) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind !Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:) ! For calculating PDSI !Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), xlon1(:), xlat1(:), Reg(:,:) Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) ! (Nlon1*Nlat1) - All gridpoints !Real, Allocatable, Dimension (:,:) :: Prec, Temp !Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp, Pet_Inp ! 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 !Real, Allocatable :: T2m(:,:,:), Pr(:,:,:) Real, Allocatable :: Pr_All(:,:,:,:,:), T2_All(:,:,:,:,:) ! Nreg,Nmem,Nlead,Ntime,Ncase Real, Allocatable :: xreg(:), xmem(:), xlead(:), xtime(:) Real, Allocatable :: SPEI_All(:,:,:,:,:,:) ! Nreg,Nmem,Nlead,Ntime,Ncase,Nlen !------------------------------------------------------------------------------ Character (Len=255) :: Fname, R_Fname, T_Fname, BL_Dir, Out_Dr_Fname Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File Character (Len=2) :: TVar, RVar Character (Len=10) :: Model_nm Logical :: status, Empty Integer :: Yr1, Yr2, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: Xt !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir, Model_nm, Mask_File, Out_Dir, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, NLen 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)) Read (1, Dr_Timescale) !Call Read_kmean ! Return Reg(Nlon,Nlat), xlon(Nlon), xlat(Nlat) Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Nst = Nreg Call rd_ncdf_2d (Mask_File,xlon,xlat,Nlon,Nlat,Regs,"reg") Allocate (Lat(Nst)) Call Estimate_Lat_for_Subreg print*,"Nst,Ntime,NLen=",Nst,Ntime,NLen !Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime)); !Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime)); !Allocate(Pet_Inp(1:Nlon,1:Nlat,1:Ntime)) ! Boc hoi 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(Nreg,Nmem,Nlead,Ntime,Ncase,Nlen)) ! Read all data: Call Read_all_data Do ilead=1,Nlead Do imem=1,Nmem Do icas=1,Ncase ! Load data: Do ist=1,Nst Do it=1,Ntime DatTMP(it,ist)=T2_All(ist,imem,ilead,it,icas) DatBL (it,ist)=Pr_All(ist,imem,ilead,it,icas) Enddo Enddo ! Calculating SPEI: !------------------------------------------------------------------------------ SPI(:,:,:) = amssng; ! SPE(:,:,:) = amssng Call Cal_SPEI ! _Reg ! <<<<<<<< Should check again if SPE is SPI or SPEI ???? Do ist=1,Nst Do it=1,Ntime Do ilen=1, Nlen SPEI_All(ist,imem,ilead,it,icas,ilen)=SPI(ist,it,ilen) Enddo Enddo Enddo ! Calculating SPI: !------------------------------------------------------------------------------ !SPI(:,:,:) = amssng !Call Cal_SPI Enddo Enddo write (Fname,"(A,I1,A)") trim(Out_Dir)//"SPEI_Reg_"//trim(Model_nm)//"_lead_",ilead,".nc" Call Write_Dr_IDXs_to_file (Fname,ilead) Enddo CONTAINS !------------------------------------------------------------------------------ Subroutine Estimate_Lat_for_Subreg Integer :: ilon,ilat, ireg, ic,j Real :: tmp Do ireg=1,Nreg tmp=0; ic=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))==ireg) then ic=ic+1; tmp=tmp+xlat(ilat) endif Enddo Enddo tmp=tmp/real(ic) Lat(ireg)=tmp Enddo End Subroutine Estimate_Lat_for_Subreg !------------------------------------------------------------------------------ Subroutine Write_Dr_IDXs_to_file (Fname,ilead) Character (Len=*) :: Fname 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 Real :: xlen(Nlen) ! SPEI_All(Nreg,Nmem,Nlead,Ntime,Ncase,Nlen) xreg,xmem,xlead,xtime t_unit="months since 1982-01-01 "; var_unit="-"; lon_nm="reg"; lat_nm="mem"; z_nm="len"; t_nm="mon" Do ilen=1,Nlen; xlen(ilen)=ilen; Enddo Allocate (X4d(Nreg,Nmem,Nlen,Ntime)) Do icas=1,Ncase Do ireg=1,Nreg Do imem=1,Nmem Do ilen=1,Nlen Do it=1,Ntime X4d(ireg,imem,ilen,it)=SPEI_All(ireg,imem,ilead,it,icas,ilen) Enddo enddo Enddo Enddo write(Var_ln,"(A,I1)") "SPEI for "//trim(Cases(icas))//" of model "//trim(Model_nm)//" with Leadtime=",ilead if (icas==1) then Call wr_ncdf_4d (Fname, xreg,xmem,xlen,xtime, Nreg,Nmem,Nlen,Ntime,X4d,trim(Cases(icas)), & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) else Call wr_ncdf_4d_apd (Fname,X4d,trim(Cases(icas)),Nreg,Nmem,Nlen,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) endif Enddo Deallocate (X4d) !Call wr_ncdf_4d (filenm,lon,lat,z,t,nlon,nlat,nz,nt,var,var_nm, & ! t_unit, var_long_name, var_unit, lon_nm, lat_nm, z_nm, t_nm, flag) !Call wr_ncdf_4d_apd (filenm,var,var_nm, nlon,nlat,nz,nt, & ! var_long_name,var_unit, lon_nm, lat_nm, z_nm, t_nm,flag) End Subroutine Write_Dr_IDXs_to_file !------------------------------------------------------------------------------ Subroutine Read_all_data Integer :: ist,Nrgs,Nmm,Nls,Nt,icas Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:) write (Fname,"(A)") trim(Inp_Dir)//trim(TVar)//"_all_reg_"//trim(Model_nm)//".nc" Do icas=1,Ncase Call rd_ncdf_4d (Fname,xreg,xmem,xlead,xtime,Nrgs,Nmm,Nls,Nt,X4d,trim(Cases(icas))) print*,"Nreg,Nmem,Nlead,Ntime=",Nrgs,Nmm,Nls,Nt if (icas==1) then; Allocate (T2_All(Nreg,Nmem,Nlead,Ntime,Ncase)); T2_All=Miss; endif T2_All(:,:,:,:,icas)=X4d(:,:,:,:) Deallocate (xreg,xmem,xlead,xtime,X4d) Enddo write (Fname,"(A)") trim(Inp_Dir)//trim(RVar)//"_all_reg_"//trim(Model_nm)//".nc" Do icas=1,Ncase Call rd_ncdf_4d (Fname,xreg,xmem,xlead,xtime,Nrgs,Nmm,Nls,Nt,X4d,trim(Cases(icas))) if (icas==1) then; Allocate (Pr_All(Nreg,Nmem,Nlead,Ntime,Ncase)); Pr_All=Miss; endif Pr_All(:,:,:,:,icas)=X4d(:,:,:,:) if (icas