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, Dat_Dir, Res_Dir Character (Len=30) :: FMT Logical :: status, Empty Integer :: Nst, Nlon, Nlat, IO, Ntime, it, YB1, YB2 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 !------------------------------------------------------------------------------ Integer, Parameter :: NFile = 30 ! 6 models + 6 models * 2 Scenarios * 2 Future Periods Character (Len=100) :: Fname_inp(NFile), Fname_Out(NFile) Namelist /Filenames/ Fname_inp, Fname_Out Integer :: ifile !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6 Real, allocatable, Dimension (:) :: Xt !------------------------------------------------------------------------------ Open (1, file = "File_List.txt"); Read (1, Filenames); Close(1) Do ifile=1, Nfile Dat_Mon(:,:,:,:) = -99.; DatNor(:,:,:) = 0.; Rain_Mon(:,:)=-99 Res_Dir = "/home/tanpv/SEACLID_New/PDSI/my_result/" if (ifile<=6) then; YB1=1986; YB2=2005; elseif (ifile<=18) then; YB1=2046; YB2=2065; elseif (ifile<=30) then; YB1=2080; YB2=2099; endif Print*, "Reading Meteo data and station info ...." FName=Fname_inp(ifile) ! Data file to be read Call Read_data_netcdf (trim(FName)) 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 !" print*,"Writing Results to file:", Fname_Out(ifile) Fname = Fname_Out(ifile) !"/home/tanpv/SEACLID_New/PDSI/Results/PDSI_APH_BL.nc" Call Write_PDSI_to_NetCDF(Fname, Nlon, Nlat, Ntime) !Call Write_HSTQ ("HSTQ_Mon.txt", TQ_M) ! Not necessary Enddo 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, iyr0 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) iyr0=YB1-1 ! Xac dinh lai nam dau cua giai doan SL Do iyr=YBas1, YBas2 iyr0=iyr0+1 write (3,"(I4,52F8.3)") iyr0,(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=*) :: 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, iyr0 if (Allocated(lon)) Deallocate(lon); if (Allocated(lat)) Deallocate(lat) Call rd_ncdf_3d (trim(FileName), lon, lat, XTime, Nlon1, Nlat1, Nt, Pr, "pr") ! Rainfall Deallocate (lon, lat, XTime) 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)); Temp_Inp(:,:,:)=T2m(:,:,:); Rain_Inp(:,:,:)=Pr(:,:,:) FileName="/home/tanpv/SEACLID_New/MASK/VN.grid_0.25.nc" if (Allocated(vietnam)) Deallocate(vietnam) Call rd_ncdf_2d (trim(FileName),loninfo,latinfo,nxlon,nxlat, vietnam, "all") 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; iyr0=YB1-1 ! Xac dinh lai nam dau cua giai doan SL Do iyr=YBas1, YBas2 ! 1, NYear iyr0=iyr0+1 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,iyr0)/25.4 ! Doi sang Monthly Total in Inch !print*, iyr,imw,it Enddo Enddo Enddo Enddo Deallocate ( XTime, loninfo,latinfo) Nst = ist; Ntime = Nt !-240; 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) :: t Character (Len=30) :: t_unit 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-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" Call wr_ncdf_3d (trim(Fname), Lon, Lat, Xt, Nlon, Nlat, Ntime, X1, "ori_pdsi", trim(t_unit), amssng) Call wr_ncdf_3d_apd (trim(Fname), X2, "ori_wplm", Nlon, Nlat, Ntime, amssng) 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