Program PDSI_SEACLID use netcdf Real, parameter :: amssng = -99. Integer, Parameter :: NVar = 2, Nst_Max = 5000 ! So Grid Cells hoac so tram Integer, Parameter :: NMonth = 12 Integer, Parameter :: YBas1 = 1986, YBas2 = 2005, NYear = YBas2-YBas1+1, maxyrs = NYear Character (Len=4) :: YT(NVar) = (/"T2m","R"/) Integer, Parameter :: NPDI = 3 ! 1=PDSI, 2= Character (Len=21) :: Station(Nst_Max) Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NMonth) :: Dat_Mon ! ivar=1: T2m ivar=2: Rain Real, Dimension (NVar, Nst_Max, NMonth) :: DatNor ! Normal Climate Condition Real :: AWC(Nst_Max), xlon(Nst_Max), xlat(Nst_Max) ! Outputs: Real, Dimension (Nst_Max,NYear*NMonth) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind Real, Dimension (Nst_Max,NYear*NMonth) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind Real, Dimension (Nst_Max,NYear*NMonth) :: Rain_Mon ! Integer, Parameter :: Nidx=7, Ntq = 7 ! Kich thuoc Ma tran Tuong quan (1..3: Mon_Ori; 4..6: Mon_Sel; 7: Rain Real, Dimension (Nst_Max,NYear*NMonth,Ntq) :: PMon ! Luu tat ca cac chuoi PDSI de tinh HSTQ Real, Dimension (Nst_Max,Ntq,Ntq) :: TQ_M ! Ma tran HSTQ Mon & Week !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, Dat_Dir, Res_Dir Character (Len=30) :: FMT Logical :: status, Empty Integer :: Nst, Nlon, Nlat, IO, Ntime, it Integer :: ist, iyr, imw, ivar, isys, system !! Lay thong tin vietnam Integer :: nxlon, nxlat, ilon, jlat Real, allocatable, dimension (:) :: loninfo, latinfo, lon, lat Real :: VN_Land(Nst_Max) ! Cac o luoi tren dat lien VN va thong tin cac vung khi hau: B1=2, B2=3,..., N3=8 Real, allocatable, dimension (:,:) :: vietnam, Xlonlat Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6 Real, allocatable, Dimension (:) :: Xt !------------------------------------------------------------------------------ Dat_Mon(:,:,:,:) = -99.; DatNor(:,:,:) = 0.; Rain_Mon(:,:)=-99 !All_IDX(:,:,:) = -99; Leng_Spel(:,:,:,:) = -99; Drought(:,:,:,:) = -99 BL_Dir = "/home/tanpv/SEACLID_New/Data/" Res_Dir = "/home/tanpv/SEACLID_New/PDSI/my_result/" Print*, "Reading Meteo data and station info ...." Call Read_data_netcdf Do ist=1, Nst it=0 Do iyr=YBas1, YBas2 Do imw=1, NMonth it=it+1 Rain_Mon(ist,it) = Dat_Mon(2,ist,iyr,imw) Enddo Enddo Enddo Call putAWC(AWC, xlon, xlat, Nst); Where (AWC<0) AWC=6.7 ! Assign to gridboxes located along coast !Call Print_test (AWC, .True.) !pause !Call Print_test (vietnam, .False.) !pause isys = system ("rm -Rf /home/tanpv/SEACLID_New/PDSI/my_result") mon_ori_pdsi(:,:)=-99; mon_ori_wplm(:,:)=-99; mon_ori_zind(:,:)=-99; mon_sel_pdsi(:,:)=-99; mon_sel_wplm(:,:)=-99; mon_sel_zind(:,:)=-99; Do ist=1, Nst if (VN_Land(ist)>=1.AND.VN_Land(ist)<=7) then ! Empty = Empty_Stn (Dat_Mon, ist) ! if (Empty) then ! print*,"Khuyet SL: "//trim(station(ist)) ! else Print*, "Calculation for station: "//trim(station(ist)) Call Write_to_PDSI ("monthly", Dat_Mon, NMonth, ist) isys = system ("./pdsi ") ! PostProcessing for station Call Read_Table (trim(Res_Dir)//"monthly/original/"//"PDSI.tbl", NMonth, mon_ori_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/original/"//"WPLM.tbl", NMonth, mon_ori_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/original/"//"ZIND.tbl", NMonth, mon_ori_zind(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"PDSI.tbl", NMonth, mon_sel_pdsi(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"WPLM.tbl", NMonth, mon_sel_wplm(ist,:)) Call Read_Table (trim(Res_Dir)//"monthly/self_cal/"//"ZIND.tbl", NMonth, mon_sel_zind(ist,:)) PMon(ist,:,1)=mon_ori_pdsi(ist,:); PMon(ist,:,2)=mon_ori_wplm(ist,:); PMon(ist,:,3)=mon_ori_zind(ist,:) PMon(ist,:,4)=mon_sel_pdsi(ist,:); PMon(ist,:,5)=mon_sel_wplm(ist,:); PMon(ist,:,6)=mon_sel_zind(ist,:) PMon(ist,:,7)=Rain_Mon(ist,:) Call HSTQ (PMon, TQ_M, ist) ! endif endif Enddo print*,"Calculating PDSI: Done !" Fname = "/home/tanpv/SEACLID_New/RES_Projection_2019/BL_OBS_VnGP_PDSI.nc" print*,"writing file NC: ", trim(Fname), " ", Nlon, Nlat, Ntime Call Write_PDSI_to_NetCDF(Fname, Nlon, Nlat, Ntime) !Call Write_HSTQ ("HSTQ_Mon.txt", TQ_M) ! Not necessary print*,"Program successfully terminated !!!" CONTAINS !------------------------------------------------------------------------------ Subroutine Write_HSTQ (Fname, X) Real, Dimension (Nst_Max,Ntq,Ntq) :: X Integer :: ist, j,k Character (Len=*) :: Fname Character (Len=2) :: Title(Ntq) = (/"op","ow","oz","sp","sw","sz","Pr"/) Open (3, file=trim(Fname), Status="Unknown") write (3, "(A)") " "//Fname Do ist=1, Nst write (3,"(A40)") trim(Station(ist)) write (3,"(8A8)") " ",(Title(j), j=1,Ntq) Do j=1,Ntq write (3,"(A4,4x,7F8.4)") Title(j),(X(ist,j,k),k=1,Ntq) Enddo Enddo Close (3) End Subroutine Write_HSTQ !------------------------------------------------------------------------------ Subroutine HSTQ (X, Y, ist) Real, Dimension (Nst_Max,NYear*NMonth,Ntq) :: X ! Luu tat ca cac chuoi PDSI de tinh HSTQ Real, Dimension (Nst_Max,Ntq,Ntq) :: Y Integer :: ist, N, it, ic, j, k Real :: TBj, TBk, Sj, Sk, Rjk, tmp N = NYear*NMonth Do j=1, Ntq Do k=1, Ntq TBj=0; TBk=0; Sj=0; Sk=0; Rjk=0; ic=0 Do it=1, N if (X(ist,it,j) /= -99. .AND. X(ist,it,k) /= -99.) then ic=ic+1 TBj=TBj+X(ist,it,j); TBk=TBk+X(ist,it,k); Sj=Sj+X(ist,it,j)*X(ist,it,j) Sk=Sk+X(ist,it,k)*X(ist,it,k); Rjk=Rjk + X(ist,it,j)*X(ist,it,k) endif Enddo tmp = ic TBj=TBj/tmp; TBk=TBk/tmp; Sj=Sj/tmp-TBj*TBj; Sk=Sk/tmp-TBk*TBk; Rjk=Rjk/tmp Y(ist,j,k) = (Rjk-TBj*TBk)/SQRT(Sj*Sk) Enddo Enddo End Subroutine HSTQ !------------------------------------------------------------------------------ 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 !print*,trim(Fname), "<<<<<<================" Do iyr=YBas1, YBas2 read (1, *) year, (wk(imw), imw=1, Nmon_week) ! print*, year, (wk(imw), imw=1, Nmon_week) Do imw=1, Nmon_week it = it+1; Outp(it) = wk(imw) Enddo Enddo isys = system ("rm -f "//trim(Fname)) End Subroutine Read_Table !------------------------------------------------------------------------------ Subroutine Write_to_PDSI (Dat_Dir, DatBL, Nmon_week, ist) Integer :: ivar, ist, iyr, imw, isys, system Real :: Tmp Character (Len=*) :: Dat_Dir Character (Len=255) :: Fname Real, Dimension (NVar, Nst_Max, YBas1:YBas2, NMonth) :: DatBL ! Write data of one station to calculate PDSI Fname = "/home/tanpv/SEACLID_New/PDSI/my_dat" isys = system ('mkdir -p '//trim(Fname)) Open (1, file=trim(Fname)//"/"//"parameter", status="Unknown") write (1,*) AWC(ist), xlat(ist) 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 Do iyr=YBas1, YBas2 Tmp = Tmp + DatBL(ivar,ist,iyr,imw) Enddo DatNor(ivar,ist,imw) = Tmp/Real(NYear) Enddo write (4,"(52F8.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week) Do iyr=YBas1, YBas2 write (3,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week) ! write (*,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week) Enddo Close (3); Close (4) Enddo Close (1) End Subroutine Write_to_PDSI !------------------------------------------------------------------------------ Subroutine Read_data_netcdf !(FileName) Character (Len=255) :: FileName Real, Allocatable :: XTime(:), T2m(:,:,:), Pr(:,:,:) !, VN(:,:) Real, allocatable, dimension (:) :: loninfo, latinfo Integer :: nxlon, nxlat, ilon, jlat Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr,imw FileName="/home/tanpv/SEACLID_New/Data/VnGP_mm_hist_allyear.nc" ! VnGP data (1986-2005) Call rd_ncdf_3d (trim(FileName), lon, lat, XTime, Nlon1, Nlat1, Nt, Pr, "pr") ! Rainfall Deallocate (lon, lat, XTime) FileName="/home/tanpv/SEACLID_New/Data/APH_mm_hist_allyear.nc" ! APH data (1986-2005) Call rd_ncdf_3d (trim(FileName), lon, lat, XTime, Nlon1, Nlat1, Nt, T2m, "tas") ! Temperature ! Luu Temperature and Precipitation for output: if (.Not.Allocated(Temp_Inp)) Allocate(Temp_Inp(1:Nlon1,1:Nlat1,1:Nt)); if (.Not.Allocated(Rain_Inp)) Allocate(Rain_Inp(1:Nlon1,1:Nlat1,1:Nt)); Where(T2m==-99.9) T2m=-99.; Temp_Inp(:,:,:)=T2m(:,:,:); Rain_Inp(:,:,:)=Pr(:,:,:) FileName="/home/tanpv/SEACLID_New/MASK/VN.grid_0.25.nc" Call rd_ncdf_2d (trim(FileName),loninfo,latinfo,nxlon,nxlat, vietnam, "all") ! Allocate (vietnam(nxlon,nxlat)); vietnam(:,:) = VN(:,:) Where (vietnam==2) vietnam=1; Where (vietnam==3) vietnam=2; Where (vietnam==4) vietnam=3; Where (vietnam==5) vietnam=4; Where (vietnam==6) vietnam=5; Where (vietnam==7) vietnam=6; Where (vietnam==8) vietnam=7 ! Gan lai cac vung KH: B1=1, B2=2, B3=3, B4=4, N1=5, N2=6, N3=7 ! Call Print_test (VN, .False.); pause print*,"In Read_data_netcdf: Nlon, Nlat, Nt=", Nlon1, Nlat1, Nt print*,lon, lat Where (Pr < 0) Pr = -99; Where (T2m < 0) T2m = -99; ist = 0 Do j=1, Nlat1 Do i=1, Nlon1 ist=ist+1; write(Station(ist),"(A,I4.4)") "Point-",ist XLon(ist)=lon(i); XLat(ist)=lat(j); VN_Land(ist)=vietnam(i,j) it=0 Do iyr=YBas1, YBas2 ! 1, NYear Do imw=1, Nmonth it=it+1 if (T2m(i,j,it)/=-99) Dat_Mon(1,ist,iyr,imw) = T2m(i,j,it)*1.8 + 32. ! Doi sang do F if (Pr (i,j,it)/=-99) Dat_Mon(2,ist,iyr,imw) = Pr (i,j,it)*DayOfMon(imw,iyr)/25.4 ! Doi sang Monthly Total in Inch Enddo Enddo Enddo Enddo Deallocate ( XTime, loninfo,latinfo) Nst = ist; Ntime = Nt; Nlon = Nlon1; Nlat = Nlat1 print*,"Nst,Ntime",Nst,Ntime End Subroutine Read_data_netcdf !------------------------------------------------------------------------------ Subroutine Write_PDSI_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 Character (Len=30) :: t_unit 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)) Do it=1, Ntime Xt(it) = it-1 Enddo ist = 0 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 t_unit = "months since 1986-01-15" print*, "writing ori_wplm ..." Call wr_ncdf_3d (trim(Fname), Lon, Lat, Xt, Nlon, Nlat, Ntime, X1, "ori_pdsi", trim(t_unit), amssng) print*, "writing ori_zind ..." Call wr_ncdf_3d_apd (trim(Fname), X2, "ori_wplm", Nlon, Nlat, Ntime, amssng) print*, "writing sel_pdsi ..." Call wr_ncdf_3d_apd (trim(Fname), X3, "ori_zind", Nlon, Nlat, Ntime, amssng) Call wr_ncdf_3d_apd (trim(Fname), X4, "sel_pdsi", Nlon, Nlat, Ntime, amssng) Call wr_ncdf_3d_apd (trim(Fname), X5, "sel_wplm", Nlon, Nlat, Ntime, amssng) Call wr_ncdf_3d_apd (trim(Fname), X6, "sel_zind", Nlon, Nlat, Ntime, amssng) ! Luu Nhiet do va luong mua dau vao: Call wr_ncdf_3d_apd (trim(Fname), Temp_Inp, "T2m", Nlon, Nlat, Ntime, amssng) Call wr_ncdf_3d_apd (trim(Fname), Rain_Inp, "Pr", Nlon, Nlat, Ntime, amssng) ! Cac vung KH Call wr_ncdf_2d_apd (trim(Fname), vietnam, "Subreg", Nlon, Nlat, amssng) End Subroutine Write_PDSI_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 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 !-------------------------------------------------------------------- Subroutine Print_test (X, One_D) ! Print test of 1-D array X(Nst) Integer :: ist, jlat, ilon Real :: X(Nst_Max) Logical :: One_D !if (One_D) then ! ist=0; if (.not.Allocated (Xlonlat)) Allocate(Xlonlat(Nlon,Nlat)) ! Do jlat=1,Nlat ! Do ilon=1,Nlon ! ist=ist+1; Xlonlat(ilon,jlat)=X(ist) ! Rain_Mon(ist,1) ! enddo ! Enddo !endif Do jlat=Nlat,1,-1 print "(40F6.1)", (Xlonlat(ilon,jlat),ilon=1,Nlon) Enddo !pause ; !if (One_D) Deallocate(Xlonlat) End Subroutine Print_test !------------------------------------------------------------------------------ Logical Function Empty_Stn (DatBL, ist) ! Tra ve .True. neu tram ist khuyet SL Real, Dimension (NVar, Nst_Max, 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 "sub_rd_ncdf.F90" Include "sub_wr_ncdf.F90" !Include "Sub_Plot.f90" !Include "Sen_Slope.f90" !-------------------------------------------------------------------- END ! Program