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 /File_Dir_Names/ Data_dir, Model_nm, MaskFname, Out_Dir, & TVar, RVar Namelist /Data_Info/ YBas1, YBas2, NLen Namelist /Dr_Timescale/ Len_SPI !------------------------------------------------------------------------------ Open (1, file="Namelist_APH_No_PDSI_2023.txt", Status="Old") 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 !Allocate (AWC(Nst)) !Call putAWC(AWC, lon, lat, Nst) !; Where (AWC<0) AWC=6.7 ! Assign to gridboxes located along coast 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 (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)) R_Fname=trim(Data_dir)//trim(RVar)//"_"//trim(Model_nm)//"_Monthly_1976-2005.nc" T_Fname=trim(Data_dir)//trim(TVar)//"_"//trim(Model_nm)//"_Monthly_1976-2005.nc" Call Read_Pr_T2m 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 (Dat_Mon(NVar,Nlon*Nlat,YBas1:YBas2, Nmonth)) !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)) DatNor(:,:,:) = 0.; ! Checking !print*,"Checking..." !print "(30F6.1)", Dat_Mon(:,:,:,:);pause ! Calculating PDSI: !------------------------------------------------------------------------------ ! Call Reform_4_PDSI !pause !Call putAWC(AWC, lon, lat, Nst) !; Where (AWC<0) AWC=6.7 ! Assign to gridboxes located along coast ! 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)//"_No_PDSI.nc" print*,"Writing ..." Call Write_All_Dr_IDX_to_NetCDF (trim(Fname)) CONTAINS !------------------------------------------------------------------------------ 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") 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(:) Integer :: Nlon1, Nlat1 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)) print*,"Nlon1, Nlat1, Nt=",Nlon1, Nlat1, Nt Call Proc_Pr_T2m Deallocate (Pr, T2m) End Subroutine Read_Pr_T2m !------------------------------------------------------------------------------ Subroutine Proc_Pr_T2m Integer, parameter :: NMon_Skip = 0 ! 72 Real :: TMean(Nmonth),RMean(Nmonth),xtas(Nyear),xpr(Nyear),xt75(Nmonth),xp75(Nmonth) Integer :: iyr,i,j,it,ict,icr,imon,imon1 Real, Allocatable :: Xtmp(:,:) print*,"Beginning Proc_Pr_T2m...:", Ntime,Nt,NMon_Skip Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,NMon_Skip+1:Nt); ! SL tu 1976 - 2005 Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr (:,:,NMon_Skip+1:Nt) ! Monthly data, mm/s ! 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 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: ! 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) ! 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 ! 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. ! 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 ! Enddo ! Enddo ! Close(3); Close(4) !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=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) if (ilen==1) then 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) else 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) endif 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, Xt) ! X3, X4, X5, X6, Xt) print*,"************ Write to file Done !" End Subroutine Write_All_Dr_IDX_to_NetCDF !------------------------------------------------------------------------------ !SUBROUTINE putAWC(dat,lon,lat,nst) ! integer :: nst, ist ! real :: dat(nst), lon(nst), lat(nst) ! integer, parameter :: ncols = 1099, nrows = 549 ! integer :: icol, irow, itmp !character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc" ! real, parameter :: xcor = -180.0, ycor = -89.9056530001 ! real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99 ! real :: indat(ncols,nrows), x(ncols), y(nrows) ! <-~~~=|==----------------------------------------------------------- 99.0 ! endif ! print*,ist, dat(ist) ! enddo !return !END SUBROUTINE putAWC !-------------------------------------------------------------------- !integer FUNCTION nrtpos(x,y,n,flag) ! integer :: n, i, itmp ! real :: flag, x(n), y, minx ! minx = 999999 ! itmp = 0 ! do i = 1, n ! if (x(i) /= flag .and. y /= flag) then ! if (abs(x(i) - y) < minx) then ! minx = abs(x(i) - y) ! itmp = i ! endif ! endif ! enddo ! nrtpos = itmp ! return !END FUNCTION nrtpos !-------------------------------------------------------------------- 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 !------------------------------------------------------------------------------ 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