Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2 Integer, Parameter :: NMonth = 12 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 :: Dat_Mon(:,:,:,:), DatNor(:,:,:) ! For calculating PDSI Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), xlon1(:), xlat1(:) 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 Real, Allocatable :: T2m(:,:,:), Pr(:,:,:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname, R_Fname, T_Fname, BL_Dir, PDSI_Dir, PDSI_WDir, Data_dir, MaskFname, AWC_Fname, Out_Dir Character (Len=3) :: TVar, RVar Character (Len=50) :: Model_nm Logical :: status, Empty Integer :: YBas1, YBas2, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it Integer :: imod, 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/ PDSI_Dir, PDSI_WDir, Data_dir, Model_nm, MaskFname, AWC_Fname, Out_Dir, & TVar, RVar Namelist /Data_Info/ YBas1, YBas2, NLen Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ !Open (1, file="Namelist_APH_2021.txt", Status="Old") ! Tinh voi APH 1961-2005 Open (1, file="Namelist_CRU_GPCC_2021.txt", Status="Old") ! Tinh voi CRU+GPCC 1951-2019 Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Len_SPI(NLen)) Read (1, Dr_Timescale) Call Proc_land_mask_AWC ! Return xlon, xlat, Nlon, Nlat, Nst, Mask, Mask_Land, Lon, Lat Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime)); Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime)); 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 (Dat_Mon(NVar,Nlon*Nlat,YBas1:YBas2, Nmonth)) Allocate (DatNor (NVar,Nst,Nmonth)) !, Rain_Mon(Nst,Ntime)) 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)) !R_Fname=trim(Data_dir)//trim(RVar)//"_"//trim(Model_nm)//"_historical_SRF.nc" !T_Fname=trim(Data_dir)//trim(TVar)//"_"//trim(Model_nm)//"_historical_SRF.nc" ! This subroutine read Preci and Temp from !Call Read_Pr_T2m R_Fname="/home/tanpv/Data/GPCC/pr_CRU_GPCC_historical_SRF_1951_2019.nc" T_Fname="/home/tanpv/Data/GPCC/tas_CRU_GPCC_historical_SRF_1951_2019.nc" Call Read_Pr_T2m_GPCC_CRU_Data print*,"Reading Data: Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst DatNor(:,:,:) = 0.; ! Checking !print*,"Checking..." !print "(30F6.1)", Dat_Mon(:,:,:,:);pause ! Calculating PDSI: !------------------------------------------------------------------------------ Call Reform_4_PDSI !pause isys = system ("rm -Rf "//trim(PDSI_Dir)) isys = system ("mkdir -p "//trim(PDSI_Dir)) mon_ori_pdsi(:,:)=Miss; mon_ori_wplm(:,:)=Miss; mon_ori_zind(:,:)=Miss; mon_sel_pdsi(:,:)=Miss; mon_sel_wplm(:,:)=Miss; mon_sel_zind(:,:)=Miss; Do ist=1, Nst if (Mask_Land(ist)>0.AND.AWC(ist)>=0) then !if (AWC(ist)<0) AWC(ist)=6.7 Print*, "Calculation for station: ", ist, lon(ist),lat(ist) Call Write_to_PDSI ("monthly",NMonth, ist) isys = system ("./pdsi.exe ") ! PostProcessing for station Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"PDSI.tbl", NMonth, mon_ori_pdsi(ist,:)) !Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"WPLM.tbl", NMonth, mon_ori_wplm(ist,:)) !Call Read_Table (trim(PDSI_Dir)//"/monthly/original/"//"ZIND.tbl", NMonth, mon_ori_zind(ist,:)) Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"PDSI.tbl", NMonth, mon_sel_pdsi(ist,:)) !Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"WPLM.tbl", NMonth, mon_sel_wplm(ist,:)) !Call Read_Table (trim(PDSI_Dir)//"/monthly/self_cal/"//"ZIND.tbl", NMonth, mon_sel_zind(ist,:)) !print "(24F8.3)", mon_ori_pdsi(ist,:) !print "(24F8.3)", mon_sel_pdsi(ist,:); pause endif Enddo print*,"Calculating PDSI: Done !" ! For calculating SPI, SPEI: ist=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1; DatBL (1:Ntime,ist) = Rain_Inp(ilon,ilat,1:Ntime) DatTMP(1:Ntime,ist) = Temp_Inp(ilon,ilat,1:Ntime) 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 Fname=trim(Out_Dir)//"BL_"//trim(Model_nm)//"_SEA.nc" Call Write_All_Dr_IDX_to_NetCDF (trim(Fname)) !Call Release_All CONTAINS !------------------------------------------------------------------------------ Subroutine Release_All DeAllocate (Rain,DatBL,DatTMP) DeAllocate (beta_spi, gamm, pzero) DeAllocate (SPI, SPE) DeAllocate (Dat_Mon, DatNor) !, AWC) DeAllocate (mon_ori_pdsi, mon_ori_wplm, mon_ori_zind) DeAllocate (mon_sel_pdsi, mon_sel_wplm, mon_sel_zind) print*,"Release all array..." End Subroutine Release_All !------------------------------------------------------------------------------ Subroutine Proc_land_mask_AWC Integer :: ist, ilon, ilat, nlo, nla Real, Allocatable :: xtmp(:,:),lontmp(:),lattmp(:) Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land") Call rd_ncdf_2d (trim(AWC_Fname),lontmp,lattmp,nlo, nla, xtmp,"awc") ! Remove points located outside SEA: Do ilon=1,Nlon if (xlon(ilon)>141.25) Mask(ilon,:) = Miss Enddo Do ilat=1,Nlat if (xlat(ilat)<11.25) Mask(:,ilat) = Miss Enddo ! Nst=Nlon*Nlat Allocate (Mask_Land(Nst)) Allocate (lon(Nst), lat(Nst), AWC(Nst)) ist=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat); AWC(ist)=xtmp(ilon,ilat) Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat) Enddo Enddo Deallocate (xtmp,lontmp,lattmp) End Subroutine Proc_land_mask_AWC !------------------------------------------------------------------------------ Subroutine Read_Table (Fname, Nmon_week, Outp) Integer :: Nmon_week, it, iyr, imw, year, isys Real :: wk(Nmon_week), Outp(Nmon_week*NYear) Character (Len=*) :: Fname Logical :: exists Inquire (File=trim(Fname),exist=exists) if (.Not.exists) return Open (1, File=trim(Fname), status="Old") it = 0 Do iyr=YBas1, YBas2 read (1, *, Err=100) year, (wk(imw), imw=1, Nmon_week) Do imw=1, Nmon_week it = it+1; Outp(it) = wk(imw) Enddo Enddo goto 101 100 Outp(:)=Miss 101 continue close (1) isys = system ("rm -f "//trim(Fname)) End Subroutine Read_Table !------------------------------------------------------------------------------ Subroutine Write_to_PDSI (Dat_Dir,Nmon_week, ist) Integer :: ivar, ist, Nmon_week, iyr, imw, isys, system, ic Real :: Tmp Character (Len=*) :: Dat_Dir ! Dat_Dir="monthly" Character (Len=155) :: Fname ! Write data of one station to calculate PDSI Fname = PDSI_WDir !"/home/tanpv/SEACLID_New/PDSI/my_dat" isys = system ('mkdir -p '//trim(Fname)) Do ivar=1, NVar if (ivar == 1) then Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_T", status="Unknown") Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_T_normal",status="Unknown") else Open (3, file=trim(Fname)//"/"//trim(Dat_Dir)//"_P", status="Unknown") Open (4, file=trim(Fname)//"/"//Dat_Dir(1:3)//"_P_normal",status="Unknown") endif Do imw=1, Nmon_week ! Tinh trung binh khi hau Tmp = 0; ic=0 Do iyr=YBas1, YBas2 if (Dat_Mon(ivar,ist,iyr,imw)/=Miss) then Tmp = Tmp + Dat_Mon(ivar,ist,iyr,imw) ic=ic+1 ! So nam co SL hop le endif Enddo if (ic>Nyear/3) then DatNor(ivar,ist,imw) = Tmp/Real(ic) else DatNor(ivar,ist,imw) = Miss endif Enddo write (4,"(52F15.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week) !print "(52F8.3)", (DatNor(ivar,ist,imw), imw=1, Nmon_week) Do iyr=YBas1, YBas2 Do imw=1, Nmon_week if (Dat_Mon(ivar,ist,iyr,imw)==Miss.AND.DatNor(ivar,ist,imw)/=Miss) & Dat_Mon(ivar,ist,iyr,imw)=DatNor(ivar,ist,imw) ! Thay bang TB nhieu nam Enddo write (3,"(I4,52F15.3)") iyr,(Dat_Mon(ivar,ist,iyr,imw), imw=1, Nmon_week) !print "(I4,52F8.3)", iyr,(Dat_Mon(ivar,ist,iyr,imw), imw=1, Nmon_week) Enddo Close (3); Close (4) Enddo Fname=trim(PDSI_WDir)//"/"//"parameter" Open (Unit=12, File=trim(Fname), Status="Unknown") write (12,"(2F20.5)") AWC(ist), lat(ist) !print "(2F20.5)", AWC(ist), lat(ist); pause Close (12) End Subroutine Write_to_PDSI !------------------------------------------------------------------------------ Subroutine Read_Pr_T2m Real, Allocatable :: XTime(:) Real, Allocatable :: xlon(:), xlat(:), X3D(:,:,:) Integer :: Nlon1, Nlat1, Nt0, it1, it2, it, it0,ilon,ilat Integer, Parameter :: Y1pr=1951, Y2pr=2007, Y1ts=1961, Y2ts=2005 !YBas2 - YBas1 Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, X3D, trim(RVar)) it1=(YBas1-Y1pr)*NMonth+1 ! Chi so thang dau tien cua chuoi it2=Nt-(Y2pr-YBas2)*NMonth if (it1<0.OR.it2>Nt) then print*,"Error in Input parameter. Please check subroutine Read_Pr_T2m" stop endif Nt0=it2-it1+1 !print*,"YBas1,Y1pr,YBas2,Y2pr,it1,it2,Nt0=",YBas1,Y1pr,YBas2,Y2pr,it1,it2,Nt0 Allocate (Pr(Nlon1, Nlat1,Nt0)) it0=0 Do it=it1,it2 it0=it0+1; Do ilon=1,Nlon1 Do ilat=1,Nlat1 Pr(ilon,ilat,it0)=X3D(ilon,ilat,it) Enddo Enddo !print*,"it,it0=",it,it0 Enddo Deallocate (xlon, xlat, XTime, X3D) !pause Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, X3D, trim(TVar)) it1=(YBas1-Y1ts)*NMonth+1 ! Chi so thang dau tien cua chuoi it2=Nt-(Y2ts-YBas2)*NMonth if (it1<0.OR.it2>Nt) then print*,"Error in Input parameter. Please check subroutine Read_Pr_T2m" stop endif Nt0=it2-it1+1 Allocate (T2m(Nlon1, Nlat1,Nt0)) it0=0 Do it=it1,it2 it0=it0+1; Do ilon=1,Nlon1 Do ilat=1,Nlat1 T2m(ilon,ilat,it0)=X3D(ilon,ilat,it) Enddo Enddo !print*,"it,it0=",it,it0 Enddo Deallocate (xlon, xlat, XTime, X3D) !pause print*,"Nlon1, Nlat1, Nt=",Nlon1, Nlat1, Nt Call Proc_Pr_T2m Deallocate (Pr, T2m) End Subroutine Read_Pr_T2m !------------------------------------------------------------------------------ Subroutine Read_Pr_T2m_GPCC_CRU_Data ! Read data from file in which data are merged from APH, CRU and GPCC ! Data from 1961-2005 Real, parameter :: Miss_pr=-99999.99, Miss_tas=9.96921e+36 Real, Allocatable :: XTime(:) Real, Allocatable :: xlon(:), xlat(:) Integer :: Nlon1, Nlat1, Nt0, it1, it2, it, it0,ilon,ilat Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, Pr, trim(RVar)) Deallocate (xlon, xlat, XTime) Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, T2m, trim(TVar)) Where (IsNaN(T2m)); T2m=Miss; Endwhere Where (IsNaN(Pr)); Pr =Miss; Endwhere Where(T2m==Miss_tas) T2m=Miss Where(Pr ==Miss_pr) Pr =Miss Temp_Inp(:,:,1:Ntime)=T2m(:,:,1:Ntime) Rain_Inp(:,:,1:Ntime)=Pr (:,:,1:Ntime) End Subroutine Read_Pr_T2m_GPCC_CRU_Data !------------------------------------------------------------------------------ Subroutine Proc_Pr_T2m !Integer, parameter :: NMon_Skip = 72 ! Khong can Real, parameter :: Miss_pr=-99999.99, Miss_tas=9.96921e+36 Real :: TMean(Nmonth),RMean(Nmonth),xtas(Nyear),xpr(Nyear),xt75(Nmonth),xp75(Nmonth) Integer :: ist,iyr,i,j,it,ict,icr,imon,imon1 Real, Allocatable :: Xtmp(:,:), X1d(:), Y1d(:) print*,"Beginning Proc_Pr_T2m...:", Ntime,Nt !,NMon_Skip !Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,1:Nt); ! SL tu 1976 - 2005 !Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr (:,:,1:Nt) ! Monthly data, mm/s !Where (IsNaN(T2m)); T2m=Miss; Endwhere !Where (IsNaN(Pr)); Pr =Miss; Endwhere !Where (T2m>50.0); T2m=Miss; Endwhere !Where(T2m==Miss_tas) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,1:Nt); ! SL tu 1976 - 2005 !Where(Pr ==Miss_pr) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr (:,:,1:Nt) ! Monthly data, mm/s ! Assign Miss values to all Sea-points: Do it=1, Ntime Do j=1, Nlat Do i=1, Nlon if (Mask(i,j)==Miss) then Rain_Inp(i,j,it)=Miss; Temp_Inp(i,j,it)=Miss endif Enddo Enddo Enddo ! Loai bo cac gridpoint khuyet SL: Allocate (X1d(Ntime),Y1d(Ntime)) ist=0 Do j=1,Nlat Do i=1,Nlon X1d(1:Ntime)=Temp_Inp(i,j,1:Ntime) Y1d(1:Ntime)=Rain_Inp(i,j,1:Ntime) ist=ist+1 if ((Empty_Stn_Arr1(X1d,Ntime)).OR.(Empty_Stn_Arr1(Y1d,Ntime))) then Mask(i,j) = Miss Mask_Land(ist)=Mask(i,j) ! Update Mask_Land and Mask endif Enddo Enddo Deallocate (X1d,Y1d) ! Fill up missing data: Do j=1, Nlat Do i=1, Nlon if (Mask(i,j) == 1) then Do imon=1, Nmonth TMean(imon)=0; RMean(imon)=0; ict=0; icr=0 Do it=1, Ntime imon1=mod(it,Nmonth); if (imon1==0) imon1=Nmonth if (imon1==imon) then if (Temp_Inp(i,j,it)/=Miss) then !TMean(imon)=TMean(imon)+Temp_Inp(i,j,it); ict=ict+1 ict=ict+1; xtas(ict)=Temp_Inp(i,j,it) endif if (Rain_Inp(i,j,it)>=0) then !/=Miss) then !RMean(imon)=RMean(imon)+Rain_Inp(i,j,it); icr=icr+1 icr=icr+1; xpr(icr)=Rain_Inp(i,j,it) endif endif Enddo !TMean(imon)=TMean(imon)/Real(ict) !RMean(imon)=RMean(imon)/Real(icr) if (ict>Nyear/3*2.AND.icr>Nyear/3*2) then TMean(imon)=Percentile(xtas,ict,0.5,Miss) RMean(imon)=Percentile(xpr, icr,0.5,Miss) xt75 (imon)=Percentile(xtas,ict,0.75,Miss) xp75 (imon)=Percentile(xpr, icr,0.75,Miss) else TMean(imon)=Miss; RMean(imon)=Miss; xt75 (imon)=Miss; xp75 (imon)=Miss endif Enddo ! Do not fill up missing values: !Do it=1, Ntime ! imon=mod(it,Nmonth); if (imon==0) imon=Nmonth ! if (Temp_Inp(i,j,it)==Miss) Temp_Inp(i,j,it)=TMean(imon) ! if (Rain_Inp(i,j,it)==Miss) Rain_Inp(i,j,it)=RMean(imon) ! if (Temp_Inp(i,j,it)>2.0*xt75(imon)) Temp_Inp(i,j,it)=xt75(imon) ! if (Rain_Inp(i,j,it)>3.0*xp75(imon)) Rain_Inp(i,j,it)=xp75(imon) !Enddo endif Enddo Enddo ! Convert Pr from mm/s to mm/month and K to C: ! Don't need for GPCC: Allocate (Xtmp(1:Nlon,1:Nlat)) !it=0 !Do iyr=YBas1, YBas2 ! Do imon=1, Nmonth ! it=it+1 !Xtmp(:,:)=Rain_Inp(:,:,it) ! Where (Xtmp/=Miss) Xtmp=Xtmp*86400.0*DayOfMon(imon,iyr) ! mm/s --> mm/month !Where (Xtmp/=Miss) Xtmp=Xtmp*DayOfMon(imon,iyr) ! mm/day --> mm/month !Rain_Inp(:,:,it)=Xtmp(:,:) !Xtmp(:,:)=Temp_Inp(:,:,it) !Where (Xtmp/=Miss) Xtmp=Xtmp-273.15 ! Don't need for APH data !Temp_Inp(:,:,it)=Xtmp(:,:) ! Enddo !Enddo ! Assign Miss values to all Sea-points: ! Do it=1, Ntime ! Do j=1, Nlat ! Do i=1, Nlon ! if (Mask(i,j)==Miss) then ! Rain_Inp(i,j,it)=Miss; Temp_Inp(i,j,it)=Miss ! endif ! Enddo ! Enddo ! Enddo Deallocate(Xtmp) End Subroutine Proc_Pr_T2m !------------------------------------------------------------------------------ Subroutine Reform_4_PDSI Integer :: ist, ilon,ilat,it,iyr,imon Open (3,file="tas.txt",status="unknown") Open (4,file="pr.txt",status="unknown") write (3,*) trim(Model_nm); write (4,*) trim(Model_nm) Dat_Mon(:,:,:,:) = Miss ist=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1 !write(3,"(A,I8)") "Stn No: ", ist; write(4,"(A,I8)") "Stn No: ", ist !print*,ist it=0 Do iyr=YBas1, YBas2 Do imon=1, Nmonth it=it+1 if (Temp_Inp(ilon,ilat,it)/=Miss) Dat_Mon(1,ist,iyr,imon) = Temp_Inp(ilon,ilat,it)*1.8+32.0 ! Doi sang do F if (Rain_Inp(ilon,ilat,it)/=Miss) Dat_Mon(2,ist,iyr,imon) = Rain_Inp(ilon,ilat,it)/25.4 ! mm/month --> inch/month Enddo Enddo if (Mask_Land(ist)>0) then write(3,"(A,I8)") "Stn No: ", ist; write(4,"(A,I8)") "Stn No: ", ist Do iyr=YBas1, YBas2 write(3,"(I4,12F8.1)") iyr, (Dat_Mon(1,ist,iyr,imon),imon=1,Nmonth) write(4,"(I4,12F8.1)") iyr, (Dat_Mon(2,ist,iyr,imon),imon=1,Nmonth) Enddo endif !if (ist==29662) then; print*,Mask_Land(ist),(Dat_Mon(1,ist,:,:)); pause; endif Enddo Enddo Close(3); Close(4) !pause End Subroutine Reform_4_PDSI !------------------------------------------------------------------------------ Subroutine Write_All_Dr_IDX_to_NetCDF(Fname) !, Nlon, Nlat, Ntime) Integer :: ist, jlat, ilon, it, ilen, ix Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm Character (Len=*) :: Fname if (.Not.Allocated(X1)) Allocate(X1(Nlon, Nlat, Ntime)) if (.Not.Allocated(X2)) Allocate(X2(Nlon, Nlat, Ntime)) if (.Not.Allocated(X3)) Allocate(X3(Nlon, Nlat, Ntime)) if (.Not.Allocated(X4)) Allocate(X4(Nlon, Nlat, Ntime)) if (.Not.Allocated(X5)) Allocate(X5(Nlon, Nlat, Ntime)) if (.Not.Allocated(X6)) Allocate(X6(Nlon, Nlat, Ntime)) if (.Not.Allocated(Xt)) Allocate(Xt(Ntime)) Do it=1, Ntime; Xt(it) = it; Enddo ist = 0 ! PDSI(Nst,Ntime): Do jlat=1, Nlat Do ilon=1, Nlon ist = ist +1 Do it = 1, Ntime X1(ilon,jlat,it) = mon_ori_pdsi(ist,it) X2(ilon,jlat,it) = mon_ori_wplm(ist,it) X3(ilon,jlat,it) = mon_ori_zind(ist,it) X4(ilon,jlat,it) = mon_sel_pdsi(ist,it) X5(ilon,jlat,it) = mon_sel_wplm(ist,it) X6(ilon,jlat,it) = mon_sel_zind(ist,it) Enddo Enddo Enddo Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01" Var_nm="ori_pdsi"; Var_ln="Origin PDSI"; var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="month" print*,"===>Writing ",trim(Var_ln) Call wr_ncdf_3d (trim(Fname), xlon, xlat, Xt, Nlon, Nlat, Ntime, X1,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) !Var_nm="ori_wplm"; Var_ln="Origin wplm" !Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) !Var_nm="ori_zind"; Var_ln="Origin zind" !Call wr_ncdf_3d_apd (trim(Fname), X3, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="sel_pdsi"; Var_ln="Self calibrate PDSI" print*,"===>Writing ",trim(Var_ln) Call wr_ncdf_3d_apd (trim(Fname), X4, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) !Var_nm="sel_wplm"; Var_ln="Self calibrate wplm" !Call wr_ncdf_3d_apd (trim(Fname), X5, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) !Var_nm="sel_zind"; Var_ln="Self calibrate zind" !Call wr_ncdf_3d_apd (trim(Fname), X6, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) ! SPI, SPE(Nst,Ntime,NLen) Do ilen=NLen, NLen !1, NLen ist=0 Do jlat=1, Nlat Do ilon=1, Nlon ist = ist +1 Do it = 1, Ntime X1(ilon,jlat,it) = SPE(ist,it,ilen) X2(ilon,jlat,it) = SPI(ist,it,ilen) Enddo Enddo Enddo write(Var_nm,"(A,I2.2)") "spei_", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPEI index in timescale of ", Len_SPI(ilen) print*,"===>Writing ",trim(Var_ln) Call wr_ncdf_3d_apd (trim(Fname), X1, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) write(Var_nm,"(A,I2.2)") "spi_", Len_SPI(ilen) write(Var_ln,"(A,I2.2)") "SPI index in timescale of ", Len_SPI(ilen) print*,"===>Writing ",trim(Var_ln) Call wr_ncdf_3d_apd (trim(Fname), X2, trim(Var_nm),Nlon, Nlat, Ntime, trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Enddo ! Luu Nhiet do va luong mua dau vao: Var_ln="Monthly T2m-temperature" print*,"===>Writing ",trim(Var_ln) 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" print*,"===>Writing ",trim(Var_ln) 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"; ! print*,"===>Writing ",trim(Var_ln) !Call wr_ncdf_2d_apd (trim(Fname), Mask, "land", nlon, nlat, trim(Var_ln), "-", trim(lon_nm), trim(lat_nm), amssng) Deallocate (X1, X2, X3, X4, X5, X6, Xt) print*,"************ Write to file Done !" End Subroutine Write_All_Dr_IDX_to_NetCDF !------------------------------------------------------------------------------ Logical Function Empty_Stn_Arr1 (X,N) ! Tra ve .True. neu tram ist khuyet SL Integer :: N Real :: X(N) Integer :: i,iyr,imon, k1, k2 Logical :: Empty Empty=.false. Do i=1,N if (X(i)==Miss) then Empty=.true.; exit endif Enddo Empty_Stn_Arr1 = Empty ! ! T hoac P khuyet SL End Function Empty_Stn_Arr1 !------------------------------------------------------------------------------ Logical Function Empty_Stn (DatBL, ist) ! Tra ve .True. neu tram ist khuyet SL Real, Dimension (NVar, Nst, YBas1:YBas2, NMonth) :: DatBL Integer :: N, ist,iyr,imon, k1, k2 Logical :: L1, L2 L1=.False.; N=0 Do iyr=YBas1, YBas2 Do imon=1,NMonth if (DatBL(1,ist,iyr,imon)<0.OR.DatBL(2,ist,iyr,imon)<0) then N = N+1 if (N>10) then ; L1=.True.; Exit; endif endif Enddo Enddo Empty_Stn = L1 !.OR.L2 ! T hoac P khuyet SL End Function Empty_Stn !------------------------------------------------------------------------------ Integer Function DayOfMon(Mon, Year) Integer :: Mon, Year Select Case (Mon) Case (1,3,5,7,8,10,12) DayOfMon = 31 Case (4,6,9,11) DayOfMon = 30 Case (2) if (mod(Year,4)==0.AND.mod(Year,100)/=0) then DayOfMon = 29 else DayOfMon = 28 endif End Select End Function DayOfMon !-------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_RW_Station_nc.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/Percentile.f90" !-------------------------------------------------------------------- END ! Program