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=220, NTC_max=3000
      Real, Parameter :: Lon1=102, Lon2=140, Lat1=6, Lat2=26
      Real, Parameter :: Dlon=2, Dlat=2, Miss=-99
      Real, Parameter :: Knots2ms=0.51444445   ! Knots to m/s
      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=255) :: Fname, St, st1
      Integer :: itc, i, io, ig, i0, iid, k1, k2, k3, k4, k5, k6
      Integer :: Nlfx, nyy, nmm, ndd, nhh, NTC, Nlf, len1, len2, i1
      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_JWTC.txt"
      Open (1, file=Fname, Status="Old")
      st1="1963090812      6       207     1297    120"
      len1=len(trim(st1))

!66666   0       19      9922    0       0       NA      NA
!1999091100      0       175     1414    15      0

      io=0; itc=0; Nlfx=0
      Do while (io==0)
         St=""
         read (1,"(A)",IOStat=IO) St
         if (IO/=0.OR.Len(trim(ST))==0) exit
         read (St,*) ig, i0, nlf, iid !, nlf
         if (ig/=66666) then
            print*,St; pause
         endif
         if (Nlfx<Nlf) Nlfx=Nlf
         if (Nlf>100) print*, "ID, nlf=",iid,nlf, St
         if (nlf>=1) then
            itc=itc+1
            TC(itc)%ID = iid
            TC(itc)%Life = nlf
            i1=0   ! So timeslice co CD>=1
            !print*,"ID, nlf=",iid,nlf
            Do i=1, nlf
               St=""; read (1,"(A)") St
               !len2=len(trim(St))
               read (St,*) k1, k2, k3, k4, k5 !, k6
               if (k2>=1) then
                  i1=i1+1
                  nyy=k1/1000000; nmm=mod(k1,1000000)/10000; ndd=mod(k1,10000)/100; nhh=mod(k1,100)
                  TC(itc)%YY(i1)=nyy; TC(itc)%MM(i1)=nmm; TC(itc)%DD(i1)=ndd; TC(itc)%HH(i1)=nhh
                  TC(itc)%CD(i1)=k2
                  TC(itc)%lat(i1)=real(k3)/10.
                  TC(itc)%lon(i1)=real(k4)/10.
                  TC(itc)%Vx(i1)=k5*Knots2ms
                  TC(itc)%Pm(i1)=Miss
               endif
            Enddo
            TC(itc)%Life = i1
            i1=1
            print"(A,3I5.4,3I2.2,I6)","itc, TC(itc)%ID, nlf= ",itc, &
                TC(itc)%ID,TC(itc)%YY(i1),TC(itc)%MM(i1),TC(itc)%DD(i1),TC(itc)%HH(i1), nlf
            if (TC(itc)%MM(i1)>12) then; print*,TC(itc)%MM(:); pause; endif
            !pause
            !print*,itc, TC(itc)%CD(1), TC(itc)%lat(1), TC(itc)%lon(1), TC(itc)%Pm(1)
         endif
      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="/home/tanpv/Data/TCs/TC_counts_grd_2x2_102-140_6-26_jwtc.nc"
      Call Write_Track_to_NC(Fname)
      Fname="/home/tanpv/Data/TCs/TC_information_NWP_jwtc.nc"
      Call Track_to_plot (Fname)
      !--------------------------------------------------------------------
      CONTAINS
      !--------------------------------------------------------------------
      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(:,:), Vx(:,:), lifetime(:) 
         Integer :: itc, ilon,ilat,it, Nlf, Nlf_x, Nx, Ny, Nz, itr
         Character (Len=50) :: var_nm, var_ln, var_unit, lon_nm, lat_nm

         Nlf_x=0
         Do itc=1, NTC
            Nlf=TC(itc)%Life
            if (Nlf_x<Nlf) Nlf_x=Nlf
         Enddo
         Nx=NTC; Ny=Nlf_x
         Allocate (tcid(NTC), tclf(Nlf_x), lontc(nx,ny), lattc(nx,ny))
         Allocate (yr(nx,ny), mn(nx,ny), day(nx,ny), hr(nx,ny))
         Allocate (inten(nx,ny), pmin(nx,ny), Vx(nx,ny), lifetime(Nx))

         Do itc=1, NTC; tcid(itc)=itc; Enddo                   ! xlon
         !Do it =1, Nlf_x; tclf(it)=(it-1)*time_inteval; Enddo  ! xlat
         Do it =1, Nlf_x; tclf(it)=it; Enddo  ! xlat

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

         Do itc=1, NTC
            Nlf=TC(itc)%Life
            lifetime(itc)=Nlf
            Do it=1,Nlf
               yr(itc,it)=TC(itc)%YY(it); mn(itc,it)=TC(itc)%MM(it)
               day(itc,it)=TC(itc)%DD(it); hr(itc,it)=TC(itc)%HH(it)
               lontc(itc,it)=TC(itc)%Lon(it); lattc(itc,it)=TC(itc)%Lat(it)
               inten(itc,it)=TC(itc)%CD(it); pmin(itc,it)=TC(itc)%Pm(it)
               Vx(itc,it)=TC(itc)%Vx(it)
            Enddo
         Enddo
         lon_nm="tc_id"; lat_nm="life_time"; var_unit="1"    ! TC_ID & Life time
         var_nm="lontc"; var_ln="Longitude of TC center"; var_unit="deg E"
         Call wr_ncdf_2d (Fname,tcid,tclf,NTC,Nlf_x,lontc,var_nm, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="lattc"; var_ln="Latitude of TC center"; var_unit="deg N"
         Call wr_ncdf_2d_apd (Fname,lattc,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="yr"; var_ln="year"; var_unit="year"
         Call wr_ncdf_2d_apd (Fname,yr,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="mon"; var_ln="Month"; var_unit="mon"
         Call wr_ncdf_2d_apd (Fname,mn,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="day"; var_ln="day"; var_unit="day"
         Call wr_ncdf_2d_apd (Fname,day,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="hrs"; var_ln="6-hourly"; var_unit="hour"
         Call wr_ncdf_2d_apd (Fname,hr,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="inten"; var_ln="TC intensity (Japanese classification)"; var_unit="1"
         Call wr_ncdf_2d_apd (Fname,inten,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="pmin"; var_ln="Minimum pressure at center"; var_unit="mb"
         Call wr_ncdf_2d_apd (Fname,pmin,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="vx"; var_ln="Maximum Wind Speed"; var_unit="m/s"
         Call wr_ncdf_2d_apd (Fname,vx,var_nm,NTC,Nlf_x, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_nm="life"; var_ln="TC life time (number of 6-hourly interval)"; var_unit="6hrs"
         Call wr_ncdf_1d_apd (Fname,lifetime,var_nm, NTC, 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="Gridded TC counts in domain 102-140E,6-26N (TC tracks)"
         var_unit="1"; 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