Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2, Nmodels=7, Nsce=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(Nmodels) Character (Len=5), Parameter :: Sce_nm(Nsce)=(/"rcp45","rcp85"/) Logical :: status, Empty Integer :: YBas1, YBas2, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it Integer :: isce, 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_SCE.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 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)) Do isce=1, Nsce Do imod=6,6 !1, Nmodels R_Fname=trim(Data_dir)//trim(RVar)//"_"//trim(Model_nm(imod))//"_RegCM4_"//Sce_nm(isce)//".nc" T_Fname=trim(Data_dir)//trim(TVar)//"_"//trim(Model_nm(imod))//"_RegCM4_"//Sce_nm(isce)//".nc" Call Read_Pr_T2m print*,"Reading Data Done: Nlat,Nlon,Nst=",Nlat,Nlon,Nst print*,"Scenario, Model: ", trim(Sce_nm(isce)),", ",trim(Model_nm(imod)) !Fname=trim(Out_Dir)//"SCE_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//"_T2m_Prec.nc" !Call Write_Prec_Temp_to_File (Fname) DatNor(:,:,:) = 0.; !------------------------------------------------------------------------------ Call Reform_4_PDSI 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 Print*, "Calculation for station: ", 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,:)) 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)//"SCE_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//".nc" Call Write_All_Dr_IDX_to_NetCDF (trim(Fname)) !Fname=trim(Out_Dir)//"SCE_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//"_T2m_Prec.nc" !Call Write_Prec_Temp_to_File (Fname) Enddo Enddo 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") 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 ! Keep all timeslices 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 2006 - 2099 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 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 Dat_Mon(:,:,:,:) = Miss ist=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1 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 Enddo Enddo 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 ====> Only Timescale=12mon to be written 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 !------------------------------------------------------------------------------ Subroutine Write_Prec_Temp_to_File (Fname) Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm Character (Len=*) :: Fname Integer :: it if (.Not.Allocated(Xt)) Allocate(Xt(Ntime)) !if (.Not.Allocated(X1)) Allocate(X1(Nlon, Nlat, Ntime)) !if (.Not.Allocated(X2)) Allocate(X2(Nlon, Nlat, Ntime)) print*,"===>Writing Prec and Temp ...",Nlon, Nlat, Ntime !X1(:,:,:)=T2m(:,:,:); X2(:,:,:)=Pr(:,:,:) Do it=1, Ntime; Xt(it) = it; Enddo Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01" Var_nm="t2m"; Var_ln="Projected monthly 2m-temperature"; var_unit="C-deg"; 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, Temp_Inp,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="pr"; Var_ln="Projected monthly Precipitation"; var_unit="mm/month" print*,"===>Writing ",trim(Var_ln) Call wr_ncdf_3d_apd (trim(Fname), Rain_Inp, trim(Var_nm), Nlon, Nlat, Ntime, trim(Var_ln), trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Deallocate (Xt) !X1,X2,Xt) End Subroutine Write_Prec_Temp_to_File !------------------------------------------------------------------------------ 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 !------------------------------------------------------------------------------ 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