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(:,:) 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, 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, 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) Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Do isce=1, Nsce Do imod=1, Nmodels R_Fname=trim(Data_dir)//trim(RVar)//"_"//trim(Model_nm(imod))//"_RegCM4_"//trim(Sce_nm(isce))//".nc" T_Fname=trim(Data_dir)//trim(TVar)//"_"//trim(Model_nm(imod))//"_RegCM4_historical.nc" Call Read_SCE_netcdf 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 (DatNor (NVar,Nst,Nmonth)) Allocate (AWC(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.; ! Calculating PDSI: !------------------------------------------------------------------------------ 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; Allocate (Mask_Land(Nst)) ist=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat); Enddo Enddo Do ist=1, Nst if (Mask_Land(ist)>0) then if (AWC(ist)<0) AWC(ist)=6.7 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 SPE(:,:,:) = SPI(:,:,:) ! Calculating SPI: !------------------------------------------------------------------------------ SPI(:,:,:) = amssng Call Cal_SPI Fname=trim(Out_Dir)//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//".nc" Call Write_All_Dr_IDX_to_NetCDF (trim(Fname)) Call Release_All Enddo ! imod Enddo ! isce 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) Deallocate (Temp_Inp,Rain_Inp,Dat_Mon) Deallocate (Mask,Mask_Land,xlon,xlat,lon,lat) End Subroutine Release_All !------------------------------------------------------------------------------ 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, *) year, (wk(imw), imw=1, Nmon_week) Do imw=1, Nmon_week it = it+1; Outp(it) = wk(imw) Enddo Enddo 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 DatNor(ivar,ist,imw) = Tmp/Real(ic) Enddo write (4,"(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) Dat_Mon(ivar,ist,iyr,imw)=DatNor(ivar,ist,imw) ! Thay bang TB nhieu nam Enddo write (3,"(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) Close (12) End Subroutine Write_to_PDSI !------------------------------------------------------------------------------ Subroutine Read_SCE_netcdf ! Data timeseries since 01/1970 to 12/2005 (36 years) Real, Allocatable :: XTime(:), T2m(:,:,:), Pr(:,:,:), Xtmp(:,:) Integer :: nxlon, nxlat, ilon, jlat Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr,imw, ict,icr, imon1 Real :: TMean(Nmonth),RMean(Nmonth) Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, Pr, trim(RVar)) ! Rainfall Deallocate (xlon, xlat, XTime) Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, T2m, trim(TVar)) ! Temperature ! Luu Temperature and Precipitation for output: if (.Not.Allocated(Temp_Inp)) Allocate(Temp_Inp(1:Nlon1,1:Nlat1,1:Ntime)); if (.Not.Allocated(Rain_Inp)) Allocate(Rain_Inp(1:Nlon1,1:Nlat1,1:Ntime)); if (.Not.Allocated(Mask)) Allocate(Mask (1:Nlon1,1:Nlat1)) Mask(:,:) = 1 Where (T2m(:,:,NMon_Skip)==-99.9) Mask = Miss ! Danh dau Lan: Mask=1, Water: Mask=Miss Where (Pr (:,:,NMon_Skip)==-99.9) Mask = Miss ! (NMon_Skip: Chi can 1 timeslice Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,1:Nt); ! SL tu 2006-2099 Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr (:,:,1:Nt) ! Monthly data, mm/s ! Fill up missing data: Do j=1, Nlat1 Do i=1, Nlon1 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 endif if (Rain_Inp(i,j,it)/=Miss) then RMean(imon)=RMean(imon)+Rain_Inp(i,j,it); icr=icr+1 endif endif Enddo TMean(imon)=TMean(imon)/Real(ict) RMean(imon)=RMean(imon)/Real(icr) 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) ENddo endif Enddo Enddo ! Convert Pr from mm/s to mm/month: Allocate (Xtmp(1:Nlon1,1:Nlat1)) it=0 Do iyr=YBas1, YBas2 Do imw=1, Nmonth it=it+1; Xtmp(:,:)=Rain_Inp(:,:,it) Where (Xtmp/=Miss) Xtmp=Xtmp*86400.0*DayOfMon(imw,iyr) Rain_Inp(:,:,it)=Xtmp(:,:) Enddo Enddo Deallocate(T2m, Pr, Xtmp); Allocate (T2m(1:Nlon1,1:Nlat1,1:Ntime),Pr(1:Nlon1,1:Nlat1,1:Ntime)) T2m = Temp_Inp; Pr = Rain_Inp ! SL da xu ly, C and mm/month print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Ntime print*,lon, lat Allocate (lon(Nlon1*Nlat1), lat(Nlon1*Nlat1)) Allocate (Dat_Mon(NVar,Nlon1*Nlat1,YBas1:YBas2, Nmonth)) ! For calculating PDSI Dat_Mon(:,:,:,:) = Miss ist=0 Do j=1, Nlat1 Do i=1, Nlon1 ist=ist+1 Lon(ist)=xlon(i); Lat(ist)=xlat(j); it=0 Do iyr=YBas1, YBas2 ! T2m = Temp_Inp; Pr = Rain_Inp : SL da xu ly, C and mm/month Do imw=1, Nmonth it=it+1 if (T2m(i,j,it)/=Miss) Dat_Mon(1,ist,iyr,imw) = T2m(i,j,it)*1.8+32. ! Doi sang do F if (Pr (i,j,it)/=Miss) Dat_Mon(2,ist,iyr,imw) = Pr(i,j,it)/25.4 ! mm/month --> inch/month Enddo Enddo Enddo Enddo Nst = ist; Nlon = Nlon1; Nlat = Nlat1 Deallocate (Xtime, T2m, Pr) print*,"Nst,Ntime",Nst,Ntime End Subroutine Read_SCE_netcdf !------------------------------------------------------------------------------ Subroutine Write_All_Dr_IDX_to_NetCDF(Fname) !, Nlon, Nlat, Ntime) Integer :: ist, jlat, ilon, it, ilen, ix ! Integer :: Nlon, Nlat, Ntime !Real, Dimension (Nlon, Nlat, Ntime) :: X1, X2, X3, X4, X5, X6 !Real, Dimension (Ntime) :: Xt Real, Allocatable :: Lon(:), Lat(:) Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm Character (Len=*) :: Fname print*, Nlon, Nlat, Ntime print*, Lon, Lat 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)) Allocate (Lon(Nlon), Lat(Nlat)) Lon(:)=xlon(:); Lat(:)=xlat(:) 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" print*, "writing ori_wplm ..." Var_nm="ori_pdsi"; Var_ln="Origin PDSI"; var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="month" Call wr_ncdf_3d (trim(Fname), Lon, Lat, 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" 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) 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) 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" 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 !------------------------------------------------------------------------------ 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 !-------------------------------------------------------------------- END ! Program