Use NetCDF
      ! Xac dinh so con bao hoat dong
      ! Neu con bao xuat hien nhieu hon 1 lan quan trac (6h mot) trong 1 o luoi thi chi tinh 1 con 
      ! de dam bao xac dinh dung quy dao bao
      Integer, Parameter :: Life_Mx=120, NTC_max=2000
      Real, Parameter :: Lon1=102, Lon2=140, Lat1=6, Lat2=26
      Real, Parameter :: Dlon=2, Dlat=2, Miss=-99
      Integer, Parameter :: Nlon=Int(Lon2-Lon1)/Int(Dlon), Nlat=Int(Lat2-Lat1)/Int(Dlat)
      Type TC_rec
           Integer :: ID
           Integer :: Life
           Integer :: YY(Life_Mx)
           Integer :: MM(Life_Mx)
           Integer :: DD(Life_Mx)
           Integer :: HH(Life_Mx)
           Real    :: Lon(Life_Mx)
           Real    :: Lat(Life_Mx)
           Integer :: CD(Life_Mx)
           Real    :: Pm(Life_Mx)
           Real    :: Vx(Life_Mx)
      End Type
      Integer :: TC1(Nlon,Nlat)
      Integer, Allocatable :: Track(:,:,:), tmp(:,:,:)
      Real,Allocatable ::     XTrack(:,:,:), xlon(:), xlat(:), xtime(:)

      Character (Len=100) :: Fname, St, st1, st2
      Integer :: itc, i, io, io1, ig, iid, len1, len2, k1, k2, k3, k4, k5, k6, k7
      Integer :: nyy, nmm, ndd, nhh, NTC, Nlf, Nmon
      Integer :: ilon, ilat, Yes_TC, it, iyr, iy, Yr1, Yr2, Nyear
      Real :: x,y,x1,y1,x2,y2
      Type (TC_rec) :: TC(NTC_max)

      ! Reading TC track data:
      Fname="/home/tanpv/Data/TCs/TC_Track_JMA.txt"
      Open (1, file=Fname, Status="Old")
      st1="77032706 002 2  38 1605 1004     000"
      len1=len(trim(st1))

      io=0; itc=0
      Do while (io==0)
         read (1,"(A)",IOStat=IO) St
         if (IO/=0.OR.Len(trim(ST))==0) exit
         itc=itc+1
         read (St,*) ig, iid, nlf
         TC(itc)%ID = iid
         TC(itc)%Life = nlf
         !print*,"TC_ID, Nlf=", TC(itc)%ID, TC(itc)%Life
         Do i=1, nlf
            read (1,"(A)") St
            !96071412 002 3 208 1408  998     035     00000 0000 90090 0090
            len2=len(trim(St))
            if (len2>=len1) then
               read (St,*) k1, k2, k3, k4, k5, k6, k7
            else
               read (St,*) k1, k2, k3, k4, k5, k6
               k7=-99
            endif
            if (k7/=-99) then 
               print"(A,7I10)","==>", k1, k2, k3, k4, k5, k6, k7
            else
            !   print*,"///", k1, k2, k3, k4, k5, k6
            endif
            nyy=k1/1000000; nmm=mod(k1,1000000)/10000; ndd=mod(k1,10000)/100; nhh=mod(k1,100)
            if (nyy>20) then
               nyy=1900+nyy
            else
               nyy=2000+nyy
            endif
            TC(itc)%YY(i)=nyy; TC(itc)%MM(i)=nmm; TC(itc)%DD(i)=ndd; TC(itc)%HH(i)=nhh
            TC(itc)%CD(i)=k3
            TC(itc)%lat(i)=real(k4)/10.
            TC(itc)%lon(i)=real(k5)/10.
            TC(itc)%Pm(i)=k6
            TC(itc)%Vx(i)=k7
         Enddo
         !print"(A,3I5.4,3I2.2,I6)","itc, TC(itc)%ID, nlf= ",itc, &
         !       TC(itc)%ID,TC(itc)%YY(1),TC(itc)%MM(1),TC(itc)%DD(1),TC(itc)%HH(1), nlf
         !print*,TC(itc)%CD(1), TC(itc)%lat(1), TC(itc)%lon(1), TC(itc)%Pm(1)
      Enddo
      NTC = itc    ! Total Number of TCs

      ! Determining the TC tracks:
      Yr1=TC(1)%YY(1); Yr2=TC(NTC)%YY(1); Nyear=Yr2-Yr1+1
      Allocate (Track(1:Nlon,1:Nlat,1:Nyear),tmp(1:Nlon,1:Nlat,1:Nyear)) 
      Track(:,:,:) =0         ! Set to Zero

      Do itc=1, NTC                   ! Search for each TC
         Nlf=TC(itc)%Life             ! A TC's life time (number of observations)
                 
         tmp(:,:,:) = 0               ! Total TC count in all years at each gridbox
         Do ilon=1, Nlon
            x1=Lon1+(ilon-1)*Dlon; x2=x1+Dlon
            Do ilat=1, Nlat
               y1=Lat1+(ilat-1)*Dlat; y2=y1+Dlat
               it=0
               Do iyr=Yr1, Yr2        ! Search for all period
                  it=it+1
                  Do i=1,Nlf          ! Is this TC located in given gridbox? 
                     iy=TC(itc)%YY(i) ! This TC(itc) occurred in iy year
                     x=TC(itc)%Lon(i); y=TC(itc)%Lat(i); Yes_TC=0 
                     Yes_TC=TC_Location (x,y,x1,y1,x2,y2,iy,iyr)  ! =0 or 1 if iyr=iy, otherwise =0
                     
                     if (Yes_TC>0) exit  ! This TC has passed in gridbox in iy year
                  Enddo
                  tmp(ilon,ilat,it)=tmp(ilon,ilat,it)+Yes_TC  ! To be increased if Yes_TC=1
               Enddo
            Enddo
         Enddo
         Track(:,:,:)=Track(:,:,:)+tmp(:,:,:) ! Update TC count for TC(itc)
      Enddo
      Do it=1, Nyear
         print*,"Year=",it+Yr1-1
         Do ilon=1,NLon
            print"(20I6)",(Track(ilon,ilat,it),ilat=Nlat,1,-1)
         Enddo
         !pause
      Enddo
      !Fname="TC_track_jma.nc"
      !Call Write_Track_to_NC(Fname)
      !Fname="TC_track_infor.nc"
      !Call Track_to_plot (Fname)  ! All track information as 4D array
      Fname="TC_tr_for_GrADS.txt"
      Call Track_to_TXT_Grads (Fname) 
      Fname="TC_B12_for_GrADS.txt"  ! Beaufort >=12
      Call Track_to_TXT_Grads_1 (Fname,5)
      Nmon=12
      Do imon=1,Nmon
         write (Fname,"(A,I2.2,A)") "TC_tr_for_GrADS_",imon,".txt"
         Call Track_to_TXT_Grads_Monthly (Fname,imon)
      Enddo
      !--------------------------------------------------------------------
      CONTAINS
      !--------------------------------------------------------------------
      Subroutine Track_to_TXT_Grads_Monthly (Fname,imon)
         ! Write out to text file for plotting tracks using grads
         Character (Len=*) :: Fname
         Integer :: itt, itc, ilon,ilat,it, Nlf, Nlf_x, Nx, Ny, Nz, itr

         N=0
         Do itc=1, NTC
            itr=Int(TC(itc)%MM(1))
            if (itr==imon) N=N+1
         Enddo

         Open (3,file=Fname, Status="unknown")
         write (3,*) N ! NTC
         itt=0
         Do itc=1, NTC
            itr=Int(TC(itc)%MM(1))
            if (itr==imon) then
               Nlf=TC(itc)%Life; itt=itt+1
               write (3,*) Nlf
               Do it=1, Nlf
                  write (3,"(10I6)") itt, Int(TC(itc)%YY(it)), Int(TC(itc)%MM(it)), &
                                       Int(TC(itc)%DD(it)), Int(TC(itc)%HH(it)), &
                                       Int(TC(itc)%Lon(it)*10), Int(TC(itc)%Lat(it)*10), &
                                       Int(TC(itc)%CD(it))
               Enddo
            endif
         Enddo
         Close (3)
      End Subroutine Track_to_TXT_Grads_Monthly
      !--------------------------------------------------------------------

      Subroutine Track_to_TXT_Grads_1 (Fname, Inten)
         ! Write out to text file for plotting tracks using grads
         ! for the TCs with specified intensity
         Character (Len=*) :: Fname
         Integer :: itc, ilon,ilat,it, Nlf, Nlf_x, Nx, Ny, Nz, itr, NTC1
         Integer :: Inten
         Logical :: yes

         Open (3,file=Fname, Status="unknown")
         NTC1=0
         Do itc=1, NTC
            Nlf=TC(itc)%Life
            yes=.false.
            Do it=1,Nlf
               if (Int(TC(itc)%CD(it))==Inten) then
                  yes=.True.; NTC1=NTC1+1;exit
               endif
            Enddo
         Enddo
         ! NTC1=Num of TC with intensity >=Inten
         write (3,*) NTC1
         Do itc=1, NTC
            Nlf=TC(itc)%Life
            yes=.false.
            Do it=1,Nlf
               if (Int(TC(itc)%CD(it))==Inten) then
                  yes=.True.; exit
               endif
            Enddo
            if (yes) then
               write (3,*) Nlf
               Do it=1, Nlf
                  write (3,"(10I6)") itc, Int(TC(itc)%YY(it)), Int(TC(itc)%MM(it)), &
                                          Int(TC(itc)%DD(it)), Int(TC(itc)%HH(it)), &
                                          Int(TC(itc)%Lon(it)*10), Int(TC(itc)%Lat(it)*10), &
                                          Int(TC(itc)%CD(it))
               Enddo
            endif
         Enddo
         Close (3)
      End Subroutine Track_to_TXT_Grads_1
      !--------------------------------------------------------------------

      Subroutine Track_to_TXT_Grads (Fname)
         ! Write out to text file for plotting tracks using grads
         Character (Len=*) :: Fname
         Integer :: itc, ilon,ilat,it, Nlf, Nlf_x, Nx, Ny, Nz, itr

         Open (3,file=Fname, Status="unknown")
         write (3,*) NTC
         Do itc=1, NTC
            Nlf=TC(itc)%Life
            write (3,*) Nlf
            Do it=1, Nlf
               write (3,"(10I6)") itc, Int(TC(itc)%YY(it)), Int(TC(itc)%MM(it)), &
                                       Int(TC(itc)%DD(it)), Int(TC(itc)%HH(it)), &
                                       Int(TC(itc)%Lon(it)*10), Int(TC(itc)%Lat(it)*10), &
                                       Int(TC(itc)%CD(it))
            Enddo
         Enddo
         Close (3)
      End Subroutine Track_to_TXT_Grads
      !--------------------------------------------------------------------
      Subroutine Track_to_plot (Fname)
         Character (Len=*) :: Fname
         Integer, Parameter :: time_inteval=6   ! hours
         Real, Allocatable :: lontc(:,:,:,:), lattc(:,:,:,:), tcid(:), tclf(:)
         Real, Allocatable :: yr(:,:,:,:), mn(:,:,:,:), day(:,:,:,:), hr(:,:,:,:)
         Real, Allocatable :: inten(:,:,:,:), pmin(:,:,:,:), lifetime(:,:,:,:)
         Real, Allocatable :: lonx(:), latx(:) 
         Integer :: itc, ilon,ilat,it, Nlf, Nlf_x, Nx, Ny, Nz, Nt, itr
         Character (Len=50) :: var_nm, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, t_unit

         Nlf_x=0
         Do itc=1, NTC
            Nlf=TC(itc)%Life
            if (Nlf_x<Nlf) Nlf_x=Nlf
         Enddo
         Nx=1; Ny=1; Nz=NTC; Nt=Nlf_x
         Allocate (tcid(Nz), tclf(Nt), lontc(nx,ny,nz,nt), lattc(nx,ny,nz,nt))
         Allocate (yr(nx,ny,nz,nt), mn(nx,ny,nz,nt), day(nx,ny,nz,nt), hr(nx,ny,nz,nt))
         Allocate (inten(nx,ny,nz,nt), pmin(nx,ny,nz,nt), lifetime(nx,ny,Nz,nt))
         Allocate (lonx(1), latx(1))
         
         lonx(1)=120; latx(1)=20               ! lon, lat
         Do itc=1, Nz; tcid(itc)=itc; Enddo    ! Z level
         !Do it =1, Nlf_x; tclf(it)=(it-1)*time_inteval; Enddo  ! xlat
         Do it =1, Nt; tclf(it)=it; Enddo      ! t

         lontc(:,:,:,:)=Miss; lattc(:,:,:,:)=Miss; yr(:,:,:,:)=Miss; mn(:,:,:,:)=Miss
         day(:,:,:,:)=Miss; hr(:,:,:,:)=Miss; inten(:,:,:,:)=Miss; pmin(:,:,:,:)=Miss
         lifetime(:,:,:,:)=Miss

         Do itc=1, Nz ! NTC
            Nlf=TC(itc)%Life
            lifetime(1,1,itc,:)=Nlf
            Do it=1,Nlf
               !lifetime(1,1,itc,it) = Nlf
               yr(1,1,itc,it)=TC(itc)%YY(it); mn(1,1,itc,it)=TC(itc)%MM(it)
               day(1,1,itc,it)=TC(itc)%DD(it); hr(1,1,itc,it)=TC(itc)%HH(it)
               lontc(1,1,itc,it)=TC(itc)%Lon(it); lattc(1,1,itc,it)=TC(itc)%Lat(it)
               inten(1,1,itc,it)=TC(itc)%CD(it); pmin(1,1,itc,it)=TC(itc)%Pm(it)
            Enddo
         Enddo

         lon_nm="lon"; lat_nm="lat"; z_nm="id"; t_nm="time";var_unit="1"
         t_unit="hours since 1951-01-01 00"

         var_nm="lontc"; var_ln="TC longitude"
         Call wr_ncdf_4d (Fname,lonx,latx,tcid,tclf,nx,ny,nz,nt,lontc,var_nm, &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
         var_nm="lattc"; var_ln="TC latitude"
         Call wr_ncdf_4d_apd (Fname, lattc, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="yr"; var_ln="year"
         Call wr_ncdf_4d_apd (Fname, yr, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="mon"; var_ln="Month"
         Call wr_ncdf_4d_apd (Fname, mn, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="day"; var_ln="day"
         Call wr_ncdf_4d_apd (Fname, day, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="hrs"; var_ln="6-hourly"
         Call wr_ncdf_4d_apd (Fname, hr, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="inten"; var_ln="TC intensity"
         Call wr_ncdf_4d_apd (Fname, inten, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_nm="pmin"; var_ln="Centers minimum pressure"
         Call wr_ncdf_4d_apd (Fname, pmin, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)

         var_nm="life"; var_ln="TC life time";
         Call wr_ncdf_4d_apd (Fname, lifetime, var_nm, nx,ny,nz,nt, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         !Call wr_ncdf_3d_apd (Fname,lifetime,var_nm,nx,ny,nz, var_ln, var_unit, lon_nm, lat_nm, z_nm, Miss)

         !Call wr_ncdf_1d_apd (Fname,lifetime,var_nm, Nz, var_ln, var_unit, lon_nm, Miss)
      End Subroutine Track_to_plot
      !--------------------------------------------------------------------

      Subroutine Write_Track_to_NC(Fname)
         Character (Len=*) :: Fname
         Integer :: ilon, ilat, it
         Character (Len=50) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm 

         Allocate (XTrack(Nlon,Nlat,Nyear), xlon(Nlon), xlat(Nlat), xtime(Nyear))
         Do ilon=1, Nlon; xlon(ilon)=(Lon1-Dlon/2)+Dlon*ilon; Enddo
         Do ilat=1, Nlat; xlat(ilat)=(Lat1-Dlat/2)+Dlat*ilat; Enddo
         Do it = 1, Nyear; xtime(it)=Yr1+it-1; Enddo
         XTrack(1:Nlon,1:Nlat,1:Nyear) = Track(1:Nlon,1:Nlat,1:Nyear)
         
         var_nm="tc";var_ln="TC tracks"; var_unit="-"; lon_nm="lon" 
         lat_nm="lat"; t_nm="time"; 
         t_unit=" years since 1951-01-01 00:00:00 "
         Call wr_ncdf_3d (Fname,xlon,xlat,xtime,Nlon,Nlat,Nyear,XTrack, &
              var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss)

      End Subroutine Write_Track_to_NC
      !--------------------------------------------------------------------
      Integer Function TC_Location (x,y,x1,y1,x2,y2,iyr,Yrs)
         Real :: x,y,x1,y1,x2,y2 
         Integer :: iyr,Yrs
         if (iyr==Yrs) then
            if (x>x1.and.x<=x2.and.y>y1.and.y<=y2) then
               !print"(6F10.3)",x1,x,x2,y1,y,y2
               TC_Location=1
            else
               TC_Location=0
            endif
         else
            TC_Location=0
         endif
      End Function TC_Location
      !--------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      END