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, Nmon=12
      Real, Parameter :: Lon1=102, Lon2=140, Lat1=6, Lat2=26
      !Real, Parameter :: VES_lon1=100, VES_lon2=120, VES_lat1=6, VES_lat2=25
      Real, Parameter :: VES_lon1=105, VES_lon2=120, VES_lat1=5, VES_lat2=23

      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(:,:,:), NTrack(:,:,:), tmp(:,:,:) ! Track and NTC in grids 
      Real,Allocatable ::     XTrack(:,:,:), xlon(:), xlat(:), xtime(:)
      Real,Allocatable ::     TC_VES(:,:), TC_VES_6m(:,:), NTC_VES(:,:), NTC_VES_6m(:,:)   ! (Year,Mon)
      Real,Allocatable ::     TC_VES_3m(:,:), NTC_VES_3m(:,:)
      Real,Allocatable ::     STC_VES(:,:)   ! Storm (CD>=3) occur in VES 

      Character (Len=100) :: Fname, St, st1
      Integer :: itc, i, io, ig, iid, k1, k2, k3, k4, k5, k6
      Integer :: nyy, nmm, ndd, nhh, NTC, Nlf, len1, len2
      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
         Do i=1, nlf
            read (1,"(A)") St
            len2=len(trim(St))
            !read (St,*) k1, k2, k3, k4, k5, k6
            if (len2>=len1) then
               read (St,*) k1, k2, k3, k4, k5, k6, k7
            else
               read (St,*) k1, k2, k3, k4, k5, k6
               k7=Miss
            endif
            !print*, k1, k2, k3, k4, k5, k6
            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
            if (k7 /= Miss) then
               TC(itc)%Vx(i)=real(k7)*Knots2ms
            else
               TC(itc)%Vx(i)=k7
            endif
         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))
      Allocate (TC_VES(Nyear,Nmon),TC_VES_6m(Nyear,Nmon))
      Allocate (NTC_VES(Nyear,Nmon),NTC_VES_6m(Nyear,Nmon))   ! TC count in VES
      Allocate (TC_VES_3m(Nyear,Nmon), NTC_VES_3m(Nyear,Nmon))
      Allocate (NTrack(1:Nlon,1:Nlat,1:Nyear))
      Allocate (STC_VES(Nyear,Nmon))

      Track(:,:,:)=0; NTrack(:,:,:)=0         ! Set to Zero
      TC_VES(:,:)=0; TC_VES_6m(:,:)=0; NTC_VES(:,:)=0; NTC_VES_6m(:,:)=0
      TC_VES_3m(:,:)=0; NTC_VES_3m(:,:)=0
      STC_VES(:,:)=0

      Call TC_count_grid
      Fname="WP_Annual_TC_NTC_grid_2x2.nc"
      Call TC_count_grid_to_NC (Fname)

      Fname="VES_Monthly_TC_NTC_Count_jp_new.nc"
      Call TC_count_VES (Fname)

      Do it=1, Nyear
         print "(I5, 12F6.0)",it+Yr1-1, (TC_VES(it,imon),imon=1,Nmon)
      Enddo
      print "(12F6.1)",(sum(TC_VES(1:Nyear,imon))/real(Nyear),imon=1,Nmon)
      !Fname="TC_track_jma.nc"
      !Call TC_Track_grid_to_NC (Fname)

      !Fname="TC_track_jma_infor.nc"
      !Call All_TC_Track_to_NC (Fname)
      !--------------------------------------------------------------------
      CONTAINS
      !--------------------------------------------------------------------
      Subroutine TC_count_VES (Fname)
      Character (Len=*) :: Fname
      Integer :: Nlf, ilon,ilat, it, iyr, i, iy, Yes_TCi, N
      Real :: x,y,x1,y1,x2,y2
      Integer, Parameter :: Sn6=6, Sn3=3   ! Tong TC trong mua
      Real, Allocatable :: x1d(:), y1d(:), xn1d(:), yn1d(:)
      Character (Len=50) :: var_nm, var_ln, var_unit, lon_nm, lat_nm

      ! Calculate the TC count over VES:
      x1=VES_lon1; x2=VES_lon2; y1=VES_lat1; y2=VES_lat2 ! VES area

      Do iyr=Yr1, Yr2
         Do imon=1, Nmon
            Do itc=1, NTC
               Nlf=TC(itc)%Life; Yes_TC=0
               Do i=1,Nlf     ! TC count
                  iy=TC(itc)%YY(i); im=TC(itc)%MM(i)   ! TC in iy and im
                  x=TC(itc)%Lon(i); y=TC(itc)%Lat(i)   ! TC's location
                  Yes_TC=TC_Location_ves (x,y,x1,y1,x2,y2,iy,iyr,im,imon)
                  if (Yes_TC>0) exit !then
               Enddo
               TC_VES(iyr-Yr1+1,imon)=TC_VES(iyr-Yr1+1,imon)+Yes_TC  ! Update after each TC
               ! Strong TC (CD>=3):
               Yes_TC=0
               Do i=1,Nlf     ! TC count
                  iy=TC(itc)%YY(i); im=TC(itc)%MM(i)   ! TC in iy and im
                  x=TC(itc)%Lon(i); y=TC(itc)%Lat(i)   ! TC's location
                  Yes_TC=TC_Location_ves (x,y,x1,y1,x2,y2,iy,iyr,im,imon)
                  if (Yes_TC>0.AND.TC(itc)%CD(i)>=3) exit !then
                  Yes_TC=0
               Enddo
               STC_VES(iyr-Yr1+1,imon)=STC_VES(iyr-Yr1+1,imon)+Yes_TC  ! Update after each TC
               ! TC days:
               Yes_TC=0
               Do i=1,Nlf     ! NTC count
                  iy=TC(itc)%YY(i); im=TC(itc)%MM(i)   ! TC in iy and im
                  x=TC(itc)%Lon(i); y=TC(itc)%Lat(i)   ! TC's location
                  
                  Yes_TC=TC_Location_ves (x,y,x1,y1,x2,y2,iy,iyr,im,imon)
                  NTC_VES(iyr-Yr1+1,imon)=NTC_VES(iyr-Yr1+1,imon)+Yes_TC
               Enddo
            Enddo !TC_loop
         Enddo
      Enddo

      N=Nyear*Nmon
      Allocate (x1d(N),y1d(N),xn1d(N),yn1d(N))
      ! Tinh cho tung mua:
      it=0
      Do iyr=1,Nyear
         Do imon=1,Nmon
            it=it+1; x1d(it)=TC_VES(iyr,imon); xn1d(it)=NTC_VES(iyr,imon)
         Enddo
      Enddo
      ! Mua 6 thang:
      y1d(:)=Miss; yn1d(:)=Miss
      Do it=1,N-Sn6+1   ! Tong truot 6 thang
         x=0; x1=0
         Do imon=it,it+Sn6-1
            x=x+x1d(imon); x1=x1+xn1d(imon)
         Enddo
         y1d(it)=x; yn1d(it)=x1
      Enddo
      it=0
      Do iyr=1,Nyear       ! Chuyen sang 2D-array
         Do imon=1,Nmon
            it=it+1; 
            TC_VES_6m(iyr,imon)=y1d(it); NTC_VES_6m(iyr,imon)=yn1d(it)
         Enddo
      Enddo
      ! Mua 3 thang:
      y1d(:)=Miss; yn1d(:)=Miss
      Do it=1,N-Sn3+1   ! Tong truot 6 thang
         x=0; x1=0
         Do imon=it,it+Sn3-1
            x=x+x1d(imon); x1=x1+xn1d(imon)
         Enddo
         y1d(it)=x; yn1d(it)=x1
      Enddo
      it=0
      Do iyr=1,Nyear       ! Chuyen sang 2D-array
         Do imon=1,Nmon
            it=it+1;
            TC_VES_3m(iyr,imon)=y1d(it); NTC_VES_3m(iyr,imon)=yn1d(it)
         Enddo
      Enddo

      ! Write to NC file:
      Deallocate (x1d,y1d,xn1d,yn1d)
      Allocate (x1d(Nyear),y1d(Nmon))
      Do iyr=1,Nyear; x1d(iyr)=Yr1+iyr-1; Enddo
      Do imon=1,Nmon; y1d(imon)=imon; Enddo
 
      lon_nm="year"; lat_nm="month";
      var_nm="tc"; var_ln="TC count over VES"; var_unit="num"
      Call wr_ncdf_2d (Fname,x1d,y1d,Nyear,Nmon,TC_VES,var_nm, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="stc"; var_ln="Monthly strong TC count (CD>=3)";
      Call wr_ncdf_2d_apd (Fname,STC_VES,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="tc3m"; var_ln="Total TC count in 3 mon from current mon";
      Call wr_ncdf_2d_apd (Fname,TC_VES_3m,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="tc6m"; var_ln="Total TC count in 6 mon from current mon";
      Call wr_ncdf_2d_apd (Fname,TC_VES_6m,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="ntc"; var_ln="NTC count over VES"; var_unit="6hrs"
      Call wr_ncdf_2d_apd (Fname,NTC_VES,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="ntc3m"; var_ln="Total NTC count in 3 mon from current mon"
      Call wr_ncdf_2d_apd (Fname,NTC_VES_3m,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      var_nm="ntc6m"; var_ln="Total NTC count in 6 mon from current mon"
      Call wr_ncdf_2d_apd (Fname,NTC_VES_6m,var_nm,Nyear,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)

      ! Write Monthly TC and strong TC to TXT:
      Open (3, file="TC_monthly_jma.txt",status="Unknown")
      write (3,*) "So luong bao va ATND hang thang" 
      write (3,"(A4,12I6)") "Nam",(imon,imon=1,Nmon)
      Do iyr=1,Nyear
         write (3,"(I4,12I6)") Yr1+iyr-1, (Int(TC_VES(iyr,imon)),imon=1,Nmon)
      Enddo
      write (3,*)
      write (3,*) "So luong bao (CD>=3) hang thang"
      write (3,"(A4,12I6)") "Nam",(imon,imon=1,Nmon)
      Do iyr=1,Nyear
         write (3,"(I4,12I6)") Yr1+iyr-1, (Int(STC_VES(iyr,imon)),imon=1,Nmon)
      Enddo
      Close(3)
      End Subroutine TC_count_VES
      !--------------------------------------------------------------------
      Subroutine TC_count_grid
      Integer :: Nlf, ilon,ilat, it, iyr, i, iy, Yes_TC
      Real :: x,y,x1,y1,x2,y2
      Real, Allocatable :: tmp1(:,:,:)

      Allocate (tmp1(1:Nlon,1:Nlat,1:Nyear)); tmp1(:,:,:)=0
      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
                  ! NTC: So obs TC trong o luoi
                  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
                     tmp1(ilon,ilat,it)=tmp1(ilon,ilat,it)+Yes_TC
                  Enddo
               Enddo
            Enddo
         Enddo
         Track(:,:,:) =Track (:,:,:)+tmp (:,:,:) ! Update TC count for TC(itc)
         NTrack(:,:,:)=NTrack(:,:,:)+tmp1(:,:,:)
      Enddo
      End Subroutine TC_count_grid
      !--------------------------------------------------------------------
      Subroutine All_TC_Track_to_NC (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="id"; lat_nm="lft"; 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"; 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 All_TC_Track_to_NC
      !--------------------------------------------------------------------

      Subroutine TC_count_grid_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)

         XTrack(1:Nlon,1:Nlat,1:Nyear) = NTrack(1:Nlon,1:Nlat,1:Nyear)
         var_nm="ntc"; var_ln="Number of TC days (6hrs)";
         Call wr_ncdf_3d_apd (Fname,XTrack,var_nm,Nlon,Nlat,Nyear, &
              var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)         
      End Subroutine TC_count_grid_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
      !--------------------------------------------------------------------
      Integer Function TC_Location_ves (x,y,x1,y1,x2,y2,iyr,Yrs,imon,Mon)
         Real :: x,y,x1,y1,x2,y2
         Integer :: iyr,Yrs,imon,Mon
         if (iyr==Yrs.AND.imon==Mon) 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_ves=1
            else
               TC_Location_ves=0
            endif
         else
            TC_Location_ves=0
         endif
      End Function TC_Location_ves
      !--------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      END