Program PDSI_SEACLID_Stn use netcdf Real, parameter :: amssng = -99. Integer, Parameter :: NVar = 2, Nst_Max = 1000 ! So Grid Cells hoac so tram Integer, Parameter :: NMonth = 12 !Integer, Parameter :: YBas1 = 1980, YBas2 = 2014, NYear = YBas2-YBas1+1, maxyrs = NYear ! For testing with origin data: Integer, Parameter :: YBas1 = 1893, YBas2 = 2003, 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 :: XH(:), Xreg(:) 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, Prec !------------------------------------------------------------------------------ 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 ...." !FName="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon.nc" ! Stn data 1966-2005 !FName="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon_Inland.nc" ! Chi lay cac tram tren dat lien FName="/home/tanpv/SEACLID_New/Data/OBS_Stn/R_T2m_1980-2014_mon_Inland_No_Rain_Center.nc" Call Read_data_netcdf (trim(FName)) ! Test with origin data: ! Call Read_data_origin 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 ! Empty = Empty_Stn (Dat_Mon, ist) ! if (Empty) then ! print*,"Khuyet SL: "//trim(station(ist)) ! else Print*, "Calculation for station: ", ist !//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,:)) if (XLat(ist)==22.817) then; print*,"AWC=", AWC(ist), mon_ori_pdsi(ist,:), Dat_Mon(:,ist,:,:); pause; endif 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 Enddo print*,"Calculating PDSI: Done !" ! Test with origin data: ! Call Write_origin_results Fname = "/home/tanpv/SEACLID_New/PDSI/Results/PDSI_Stn_No_Rain_Center.nc" Call Write_PDSI_to_NetCDF(Fname) !Call Write_HSTQ ("HSTQ_Mon.txt", TQ_M) ! Not necessary CONTAINS ! Include "Read_data_origin.f90" ! Fore testing only !------------------------------------------------------------------------------ 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, ic 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) ! if (XLon(ist)==104.967.AND.XLat(ist)==22.817) then; print*,"AWC=", AWC(ist), xlat(ist); pause; endif 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 (DatBL(ivar,ist,iyr,imw) /= amssng) then; Tmp = Tmp + DatBL(ivar,ist,iyr,imw); ic=ic+1; endif Enddo DatNor(ivar,ist,imw) = Tmp/real(ic) ! Real(NYear) Enddo write (4,"(20F10.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week) Do iyr=YBas1, YBas2 write (3,"(I4,20F10.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week) if (XLon(ist)==104.967.AND.XLat(ist)==22.817) then write (*,"(I4,52F8.3)") iyr,(DatBL(ivar,ist,iyr,imw), imw=1, Nmon_week) ! pause; endif Enddo if (XLon(ist)==104.967.AND.XLat(ist)==22.817) then; write (*,"(52F8.3)") (DatNor(ivar,ist,imw), imw=1, Nmon_week); pause ; endif 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 :: XH(:), Xreg(:) Real, allocatable, dimension (:) :: loninfo, latinfo Integer :: nxlon, nxlat, ilon, jlat Integer :: Nlon1, Nlat1, Nt, i, j, ist, it, iyr,imw Call read_netcdf_station(trim(FileName), Nst, Nt, Pr, "R", XTime, lon, lat, Xh, Xreg) Deallocate (lon, lat, Xh, Xreg) Call read_netcdf_station(trim(FileName), Nst, Nt, T2m, "T2m", XTime, lon, lat, Xh, Xreg) 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: Nst, Nt=", Nst, Nt print*,lon, lat Where (IsNaN(Pr)); Pr=amssng; Endwhere; Where (IsNaN(T2m)); T2m=amssng; Endwhere Call FillUp_Miss_Data (Pr, Nt, Nst); Call FillUp_Miss_Data (T2m, Nt, Nst) Allocate (Temp(Nst, Nt), Prec(Nst, Nt)); Temp(:,:)=T2m(:,:); Prec(:,:)=Pr(:,:) ! Luu lai de ghi ra file KQ Do ist=1, Nst XLon(ist)=lon(ist); XLat(ist)=lat(ist); it=0 Do iyr=YBas1, YBas2 Do imw=1, Nmonth it=it+1 if (.not.IsNaN(T2m(ist,it))) then; Dat_Mon(1,ist,iyr,imw) = T2m(ist,it)*1.8 + 32.; else; Dat_Mon(1,ist,iyr,imw) = amssng; endif ! Doi sang do F if (.not.IsNaN(Pr (ist,it))) then; Dat_Mon(2,ist,iyr,imw) = Pr (ist,it)/25.4 ; else; Dat_Mon(2,ist,iyr,imw) = amssng; endif ! Doi sang Inch Enddo Enddo Enddo print*,"Nst,Ntime",Nst,Nt Deallocate (loninfo,latinfo) Ntime = Nt; print*,"Nst,Ntime",Nst,Ntime End Subroutine Read_data_netcdf !------------------------------------------------------------------------------ Subroutine FillUp_Miss_Data (X, Ntime, Nst) Integer :: ist, it, imon, it1, iyr, Nst, Ntime Real :: X(Nst,Ntime), Mean(NMonth) Do ist=1, Nst Do imon=1, Nmonth Mean(imon)=0;iyr=0 Do it=imon, Ntime, Nmonth if (X(ist,it)/=amssng) then; Mean(imon)=Mean(imon)+X(ist,it); iyr=iyr+1;endif Enddo Mean(imon)=Mean(imon)/real(iyr) Do it=imon, Ntime, Nmonth if (X(ist,it)==amssng) then; print*,X(ist,it); X(ist,it)=Mean(imon); print*,X(ist,it); pause;endif Enddo Enddo Enddo End Subroutine FillUp_Miss_Data !------------------------------------------------------------------------------ Subroutine Write_PDSI_to_NetCDF(Fname) Integer :: ist, jlat, ilon, it, ilen, ix Real, Allocatable, Dimension (:,:) :: X1, X2, X3, X4, X5, X6 Real, Allocatable, Dimension (:) :: t Real, Allocatable, Dimension (:) :: XStn Character (Len=30) :: t_unit Character (Len=*) :: Fname Allocate (X1(Nst,Ntime), X2(Nst,Ntime), X3(Nst,Ntime), X4(Nst,Ntime), X5(Nst,Ntime), X6(Nst,Ntime), XStn(Nst),t(Ntime)) Do it=1, Ntime; t(it) = it-1; Enddo; Do ist=1, Nst; XStn(ist) = ist; Enddo Do ist=1, Nst Do it = 1, Ntime X1(ist,it) = mon_ori_pdsi(ist,it) X2(ist,it) = mon_ori_wplm(ist,it) X3(ist,it) = mon_ori_zind(ist,it) X4(ist,it) = mon_sel_pdsi(ist,it) X5(ist,it) = mon_sel_wplm(ist,it) X6(ist,it) = mon_sel_zind(ist,it) Enddo ! print*,X1; pause Enddo t_unit = "months since 1980-01-15" Call wr_ncdf_2d (trim(Fname), XStn, t, Nst, Ntime, X1, "ori_pdsi", "ori_pdsi", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), X2, "ori_wplm", Nst, Ntime, "ori_wplm", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), X3, "ori_zind", Nst, Ntime, "ori_zind", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), X4, "sel_pdsi", Nst, Ntime, "sel_pdsi", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), X5, "sel_wplm", Nst, Ntime, "sel_wplm", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), X6, "sel_zind", Nst, Ntime, "sel_zind", "-", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), Temp, "T2m" , Nst, Ntime, "2m-Temperature","Deg. C", "Station", "Time", amssng) Call wr_ncdf_2d_apd (trim(Fname), Prec, "Pre" , Nst, Ntime, "Precipitation", "mm/mon", "Station", "Time", amssng) Call wr_ncdf_1d_apd (trim(Fname), lon, "lon", Nst, "longitude", "deg. E", "Station", amssng) Call wr_ncdf_1d_apd (trim(Fname), lat, "lat", Nst, "latitude", "deg. N", "Station", amssng) Call wr_ncdf_1d_apd (trim(Fname), Xreg, "reg", Nst, "region", "-", "Station", 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 !-------------------------------------------------------------- Include "../../Libs/sub_wr_ncdf_Tan.F90" Include "../../Libs/Sub_RW_Station_nc.f90" Include "sub_rd_ncdf.F90" !Include "sub_wr_ncdf.F90" !Include "Sub_Plot.f90" !Include "Sen_Slope.f90" !-------------------------------------------------------------------- END ! Program