Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2 Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI 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, Dimension (:,:) :: Rain_Mon Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:) ! For calculating PDSI Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:) Real, Allocatable :: Mask_Land(:), Lon(:), Lat(:) ! (Nlon1*Nlat1) - All gridpoints Real, Allocatable, Dimension (:,:) :: Prec, Temp Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_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 !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Mask_File, Out_Dr_Fname Character (Len=3) :: TVar, RVar Character (Len=30) :: FMT Logical :: status, Empty Integer :: YBas1, YBas2, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, it Integer :: ist, iyr, imw, ivar, isys, system, ilon, ilat !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6 Real, allocatable, Dimension (:) :: Xt !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ R_Fname, T_Fname, Mask_File, Out_Dr_Fname, & TVar, RVar Namelist /Data_Info/ YBas1, YBas2, NLen Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_GPCC_CRU.txt", Status="Old") Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Len_SPI(NLen)) Read (1, Dr_Timescale) Close (1) Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Call Read_GPCC_CRU print*,"Reading Data: Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst 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 (Lat(Nst)) !Allocate (DatNor (NVar,Nst,Nmonth)) !, Rain_Mon(Nst,Ntime)) !Allocate (AWC(Nst)) !, xlon(Nst), xlat(Nst)) !Allocate (mon_ori_pdsi(Nst,Ntime), mon_ori_wplm(Nst,Ntime), mon_ori_zind(Nst,Ntime)) !Allocate (mon_sel_pdsi(Nst,Ntime), mon_sel_wplm(Nst,Ntime), mon_sel_zind(Nst,Ntime)) !------------------------------------------------------------------------------ ! For calculating SPI, SPEI: ist=0 Do ilat=1, Nlat Do ilon=1, Nlon if (Mask(ilon,ilat)/=Miss) then ist=ist+1; DatBL (1:Ntime,ist) = Rain_Inp(ilon,ilat,1:Ntime) DatTMP(1:Ntime,ist) = Temp_Inp(ilon,ilat,1:Ntime) Lat(ist)=xlat(ilat) !print"(30F6.1)", DatBL (Ntime-10:Ntime,ist) !print"(30F6.1)", DatTMP (Ntime-10:Ntime,ist) endif Enddo Enddo ! Calculating SPI: !------------------------------------------------------------------------------ SPI(:,:,:) = amssng; SPE(:,:,:) = amssng Call Cal_SPEI ! <<<<<<<< Should check again if SPE is SPI or SPEI ???? SPE(:,:,:) = SPI(:,:,:) ! Calculating SPI: !------------------------------------------------------------------------------ SPI(:,:,:) = amssng Call Cal_SPI print*,"Writing to file ..." Call Write_All_Dr_IDX_to_NetCDF (trim(Out_Dr_Fname)) CONTAINS !------------------------------------------------------------------------------ Subroutine Read_GPCC_CRU Integer, Parameter :: Yr1=1951, Yr2=2019 Real, Parameter :: Miss_gpcc=-99999.99, Miss_cru=9.96921e+36 Real, Allocatable :: xtime(:), T2m(:,:,:), Pr(:,:,:), x1d(:), y1d(:) Integer :: Nt, Nt1, Nt2, ilon, ilat, it, it1, it2 ! xlon(:), xlat(:), Mask(:,:) Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, xtime, Nlon, Nlat, Nt1, Pr, trim(RVar)) Deallocate (xlon, xlat, xtime) Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, xtime, Nlon, Nlat, Nt2, T2m, trim(TVar)) Deallocate (xlon, xlat) Call rd_ncdf_2d (trim(Mask_File),xlon,xlat,nlon,nlat,Mask,"land") Where (Pr==Miss_gpcc) Pr=Miss Where (T2m==Miss_cru) T2m=Miss it1=(YBas1-Yr1)*NMonth+1 ! 1960-01 it2=it1+Ntime-1 ! 2015-12 print*,"Nlon,Nlat,Ntime,Nt1,Nt2,it1,it2=",Nlon,Nlat,Ntime,Nt1,Nt2,it1,it2 Allocate (Temp_Inp(Nlon,Nlat,Ntime),Rain_Inp(Nlon,Nlat,Ntime)) Allocate (x1d(Ntime),y1d(Ntime)) Temp_Inp=Miss; Rain_Inp=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then x1d(1:Ntime)=Pr(ilon,ilat,it1:it2) Call Fill_up (x1d,Ntime) Rain_Inp(ilon,ilat,1:Ntime)=x1d(1:Ntime) x1d(1:Ntime)=T2m(ilon,ilat,it1:it2) Call Fill_up (x1d,Ntime) Temp_Inp(ilon,ilat,1:Ntime)=x1d(1:Ntime) endif Enddo Enddo Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then x1d(1:Ntime)=Rain_Inp(ilon,ilat,1:Ntime) y1d(1:Ntime)=Temp_Inp(ilon,ilat,1:Ntime) if (.not.FullDara(x1d,Ntime,Miss).AND..not.FullDara(y1d,Ntime,Miss)) then Mask(ilon,ilat)=Miss Rain_Inp(ilon,ilat,1:Ntime)=Miss Temp_Inp(ilon,ilat,1:Ntime)=Miss endif endif Enddo Enddo Nst=0 Do ilat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,ilat)/=Miss) Nst=Nst+1 Enddo Enddo End Subroutine Read_GPCC_CRU !------------------------------------------------------------------------------ Subroutine Fill_up (X,Ntime) Real :: X(Ntime), TB(Nmon) Integer :: ic, it, imon, imon1, Nleng Real :: tmp if (mod(Ntime,Nmon)/=0) print*,"Something wrong !" Nleng=Int(0.8*Ntime) ! 80% do dai chuoi co SL Do imon=1,Nmon tmp=0; ic=0 Do it=1,Ntime imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon if (imon1==imon) then if (X(it)/=Miss) then ic=ic+1; tmp=tmp+X(it) endif endif Enddo if (ic>=Nleng) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif Enddo Do it=1,Ntime imon=mod(it,Nmon); if (imon==0) imon=Nmon if (TB(imon)/=Miss.AND.X(it)==Miss) X(it)=TB(imon) Enddo End Subroutine Fill_up !------------------------------------------------------------------------------ Logical Function FullDara(X,N,Miss) Integer :: N, i Real :: X(N), Miss Logical :: Full Full=.True. Do i=1,N if (X(i)==Miss) then; Full=.False.; exit; endif Enddo FullDara=Full End Function FullDara !------------------------------------------------------------------------------ Subroutine Write_All_Dr_IDX_to_NetCDF(Fname) !, Nlon, Nlat, Ntime) Integer :: ist, jlat, ilon, it, it1, ilen, ix, Nt Real, Allocatable :: Lon(:), Lat(:), Xt(:), X1(:,:,:), X2(:,:,:) Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm Character (Len=*) :: Fname Nt = Ntime !-Nmon ! Bo nam dau tien print*, Nlon, Nlat, Ntime, Nt, Nlen Allocate(X1(Nlon, Nlat, Nt)) Allocate(X2(Nlon, Nlat, Nt)) Allocate (Lon(Nlon), Lat(Nlat), Xt(Nt)) Lon(:)=xlon(:); Lat(:)=xlat(:); X1=Miss; X2=Miss print*, Lon, Lat Do it=1, Nt; Xt(it) = it; Enddo Write(t_unit,"(A,I4.4,A)") "months since ",YBas1+1,"-01-01" Do ilen=1, NLen ist=0 Do jlat=1,Nlat Do ilon=1,Nlon if (Mask(ilon,jlat)/=Miss) then ist=ist+1 Do it = 1,Ntime ! Lay tu nam thu 2 X1(ilon,jlat,it) = SPE(ist,it,ilen) X2(ilon,jlat,it) = SPI(ist,it,ilen) Enddo endif Enddo Enddo var_unit="1"; lon_nm="lon"; lat_nm="lat"; t_nm="month" if (ilen==1) then write(Var_nm,"(A,I2.2)") "spei", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPEI index in timescale of ", Len_SPI(ilen) Call wr_ncdf_3d (trim(Fname),Lon,Lat,Xt,Nlon,Nlat, Nt, X1,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),Miss) write(Var_nm,"(A,I2.2)") "spi", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPI index in timescale of ", Len_SPI(ilen) Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon,Nlat,Nt,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),Miss) else write(Var_nm,"(A,I2.2)") "spei", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPEI index in timescale of ", Len_SPI(ilen) Call wr_ncdf_3d_apd (trim(Fname), X1, trim(Var_nm),Nlon,Nlat,Nt,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),Miss) write(Var_nm,"(A,I2.2)") "spi", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPI index in timescale of ", Len_SPI(ilen) Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon,Nlat,Nt,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),Miss) endif Enddo ! Luu Nhiet do va luong mua dau vao: Var_ln="Monthly T2m-temperature" Call wr_ncdf_3d_apd (trim(Fname), Temp_Inp, "T2m", Nlon, Nlat, Ntime, trim(Var_ln), "Deg.C",trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_ln="Monthly Precipitation Count" Call wr_ncdf_3d_apd (trim(Fname), Rain_Inp, "Pr", Nlon, Nlat, Ntime, trim(Var_ln), "mm/month",trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) ! Luu Mask_Land: Var_ln="Land mask: =1: Land, =-99: Sea/water"; Call wr_ncdf_2d_apd (trim(Fname), Mask, "land", nlon, nlat, trim(Var_ln), "-", trim(lon_nm), trim(lat_nm), amssng) End Subroutine Write_All_Dr_IDX_to_NetCDF !------------------------------------------------------------------------------ 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