Program cal_freq
! by: hapt
! to: calculate spatial distribution of the TC frequency
Implicit none
!use netcdf

Character*100:: cfile
Integer, parameter:: nstep=100
Integer:: ntc, itc, istep, itime
Real, dimension(:,:),allocatable:: tc_lon, tc_lat, tc_wind
Integer, dimension(:,:),allocatable:: tc_day, tc_mon, tc_year, tc_hour
Integer, dimension(:),allocatable:: tc_step
Integer, parameter:: smon=6, emon=8  ! thang nghien cuu
Real, parameter:: rmis=-99
Character*100:: idir
Character:: ctmp
! create map to make spatial frequency distribution
Real, parameter:: slon=100.0, elon=180 ! gioi han kinh do
Real, parameter:: slat=0, elat=40 ! gioi han vi do
Real, parameter:: cint=2.5 !resolution
Integer:: nlon, nlat
Integer:: ilon, ilat, itmp
Real, dimension(:), allocatable:: lon, lat
Real, dimension(:,:), allocatable:: tc_freq
Real:: slope, intercept
Real, dimension(:,:), allocatable:: flag
Real, parameter:: tmp_cint=0.05
Real:: lt_lon, gt_lon
Real:: lt_lat, gt_lat
Real:: tmp_lon, tmp_lat
Integer,parameter:: ntime=1
Integer:: life_length
Real,dimension(:,:,:), allocatable:: out_var
Integer:: total_TC
Character*20:: var_name
Character*100:: fname, oDir
Namelist/nml/iDir,cfile,oDir,life_length

Open(22,file="info.txt",status="unknown")
      Read(22,nml)
Close(22)

var_name="freq"
nlon=((elon-slon)/cint)+1
nlat=((elat-slat)/cint)+1

Allocate(lon(nlon))
Allocate(lat(nlat)) 
Allocate(tc_freq(nlon-1,nlat-1))
Allocate(flag(nlon-1,nlat-1))

Do ilon=1,nlon
   lon(ilon)=slon+(ilon-1)*cint
   !write(*,*) lon(ilon)
End do
Do ilat=1,nlat
   lat(ilat)=slat+(ilat-1)*cint
   !write(*,*) lat(ilat)
End do

Open(1,file=trim(oDir)//trim(cfile)//".txt",status="unknown")
     Read(1,*) ctmp, ntc
     Allocate(tc_lon(ntc,nstep))
     Allocate(tc_lat(ntc,nstep))
     Allocate(tc_day(ntc,nstep))
     Allocate(tc_mon(ntc,nstep))
     Allocate(tc_year(ntc,nstep))
     Allocate(tc_hour(ntc,nstep))
     Allocate(tc_step(ntc))
     tc_lon=rmis
     tc_lat=rmis
     Do itc=1,ntc
        Read(1,*) ctmp
        Read(1,*) ctmp, tc_step(itc)
        Read(1,*) ctmp, tc_step(itc)
        Do istep=1,tc_step(itc)
           Read(1,*) tc_year(itc,istep), tc_mon(itc,istep), tc_day(itc,istep), tc_hour(itc,istep), tc_lon(itc,istep), tc_lat(itc,istep)
        End do
     End do
Close(1)

tc_freq=0
total_TC=0
Do itc=1,ntc
   flag=0
   !If (tc_mon(itc,1).ge.smon.and.tc_mon(itc,1).le.emon) then
  If (tc_step(itc).ge.life_length) then
  total_TC=total_TC+1
  If (tc_step(itc).eq.1) then
          istep=tc_step(itc)
          Do ilon=1,nlon-1
             If (lon(ilon).le.tc_lon(itc,istep).and.tc_lon(itc,istep).le.lon(ilon+1)) then
                 Do ilat=1,nlat-1
                     If(lat(ilat).le.tc_lat(itc,istep).and.tc_lat(itc,istep).le.lat(ilat+1)) then
                         tc_freq(ilon,ilat)=tc_freq(ilon,ilat)+1
                         flag(ilon,ilat)=1
                     End if
                 End do
             End if
          End do
  Else

   Do istep=1,tc_step(itc)-1
      !flag=0
    
      ! ########### xlon1 # xlon2; ylat1 # ylat2 #############
      If (tc_lat(itc,istep).ne.tc_lat(itc,istep+1).and.tc_lon(itc,istep).ne.tc_lon(itc,istep+1)) then
       
          slope=(tc_lat(itc,istep+1)-tc_lat(itc,istep))/(tc_lon(itc,istep+1)-tc_lon(itc,istep))
          intercept=-(tc_lon(itc,istep)*(tc_lat(itc,istep+1)-tc_lat(itc,istep)))/(tc_lon(itc,istep+1)-tc_lon(itc,istep))+tc_lat(itc,istep)
          If (tc_lon(itc,istep).lt.tc_lon(itc,istep+1)) then
              lt_lon=tc_lon(itc,istep)
              gt_lon=tc_lon(itc,istep+1)
          Else
              lt_lon=tc_lon(itc,istep+1)
              gt_lon=tc_lon(itc,istep)
          End if
          itmp=0
          tmp_lon=lt_lon
         
          Do while (tmp_lon.le.(gt_lon+tmp_cint)) 
    
         
             itmp=itmp+1
             tmp_lon=lt_lon+(itmp-1)*tmp_cint
             tmp_lat=tmp_lon*slope+intercept

             Do ilon=1,nlon-1
                If (lon(ilon).le.tmp_lon.and.tmp_lon.le.lon(ilon+1)) then
                    Do ilat=1,nlat-1 
                       If (lat(ilat).le.tmp_lat.and.tmp_lat.le.lat(ilat+1)) then
                           If (flag(ilon,ilat).eq.0) then
                               tc_freq(ilon,ilat)=tc_freq(ilon,ilat)+1
                               flag(ilon,ilat)=1
                           End if
                       End if
                    End do
                End if
             End do
         End do
      End if

      ! ################ xlon1 = xlon2 ####################
      If (tc_lon(itc,istep).eq.tc_lon(itc,istep+1)) then  
          tmp_lon=tc_lon(itc,istep)
          If (tc_lat(itc,istep).le.tc_lat(itc,istep+1)) then
              lt_lat=tc_lat(itc,istep)
              gt_lat=tc_lat(itc,istep+1)
          Else
              lt_lat=tc_lat(itc,istep+1)
              gt_lat=tc_lat(itc,istep)
          End if
          tmp_lat=lt_lat
          
          Do ilon=1,nlon-1
             If (lon(ilon).le.tmp_lon.and.tmp_lon.le.lon(ilon+1)) then    
                  
                  itmp=0
                  Do while (tmp_lat.le.(gt_lat+tmp_cint))
                     itmp=itmp+1
                     tmp_lat=lt_lat + (itmp-1)*tmp_cint

                     Do ilat=1,nlat-1
                         If (lat(ilat).le.tmp_lat.and.tmp_lat.le.lat(ilat+1)) then
                            If (flag(ilon,ilat).eq.0) then
                               tc_freq(ilon,ilat)=tc_freq(ilon,ilat)+1
                               flag(ilon,ilat)=1
                            End if
                         End if
                     End do
                  End do
             End if
          End do
      End if

     ! ################ ylat1 = ylat2 ####################
      If (tc_lat(itc,istep).eq.tc_lat(itc,istep+1)) then
          tmp_lat=tc_lat(itc,istep)
          If (tc_lon(itc,istep).le.tc_lon(itc,istep+1)) then
              lt_lon=tc_lon(itc,istep)
              gt_lon=tc_lon(itc,istep+1)
          Else
              lt_lon=tc_lon(itc,istep+1)
              gt_lon=tc_lon(itc,istep)
          End if
          tmp_lon=lt_lon
          Do ilat=1,nlat-1
             If (lat(ilat).le.tmp_lat.and.tmp_lat.le.lat(ilat+1)) then
                  itmp=0
                  Do while (tmp_lon.le.(gt_lon+tmp_cint))
                     itmp=itmp+1
                     tmp_lon=lt_lon + (itmp-1)*tmp_cint
                     Do ilon=1,nlon-1
                         
                         If (lon(ilon).le.tmp_lon.and.tmp_lon.le.lon(ilon+1)) then
                             If (flag(ilon,ilat).eq.0) then
                                 tc_freq(ilon,ilat)=tc_freq(ilon,ilat)+1
                                 flag(ilon,ilat)=1
                             End if
                         End if
                     End do
                  End do
              End if
           End do
        End if
   End do
   
   ! ################ for the last time step ####################
   istep=tc_step(itc)
   Do ilon=1,nlon-1
      If (lon(ilon).le.tc_lon(itc,istep).and.tc_lon(itc,istep).le.lon(ilon+1)) then
          Do ilat=1,nlat-1
             If (lat(ilat).le.tc_lat(itc,istep).and.tc_lat(itc,istep).le.lat(ilat+1)) then
              If (flag(ilon,ilat).eq.0) then
                 flag(ilon,ilat)=1
                 tc_freq(ilon,ilat)=tc_freq(ilon,ilat)+1
              End if
             End if
          End do
      End if
   End do
  
  End if ! for select ntc_step = 1
  End if ! for select TC with tc_step(itc) - ge life_length
End do
    

Open(4,file=trim(oDir)//"TC_freq_"//trim(cfile)//".txt",status="unknown")
     Do ilon=1,nlon-1
        Write(4,'(221(F5.1,2x))') (tc_freq(ilon,ilat),ilat=1,nlat-1)
     End do
Close(4)

Open(5,file=trim(oDir)//"lon.list",status="unknown")
     Do ilon=1,nlon-1
        Write(5,'(F7.2)') (lon(ilon)+lon(ilon+1))/2  
     End do
Close(5)

Open(6,file=trim(oDir)//"lat.list",status="unknown")
     Do ilat=1,nlat-1
        Write(6,'(F7.2)') (lat(ilat)+lat(ilat+1))/2
     End do
Close(6)

Open(7,file=trim(oDir)//"name.list",status="unknown")
     Write(7,'(a1000)') trim(oDir)//"TC_freq_"//trim(cfile)
Close(7)

Open(8,file="Total_TC.txt",status='unknown',access='append')
     write(8,'(I6)') total_TC
Close(8)
! Write ncdf file
!Allocate(out_var(nlon,nlat,ntime))
!itime=1
!Do ilon=1,nlon
!   Do ilat=1,nlat
!      out_var(ilon,ilat,itime)=tc_freq(ilon,ilat)
!   End do
!End do

!fname="TC_freq_"//trim(cfile)//".nc"
!call write_nc3d(fname,var_name,nlon,nlat,ntime,lon,lat,out_var)
     
!Do itc=1,ntc
!   Write(*,*) itc, tc_step(itc)
!End do

!itc=3
!Do istep=1,tc_step(itc)
!   Write(*,"2(I2,3x),i4,2x,I2,2x,3(F6.1,2x)") tc_day(itc,istep), tc_mon(itc,istep), tc_year(itc,istep), tc_hour(itc,istep), tc_lat(itc,istep), tc_lon(itc,istep), tc_wind(itc,istep)
!End do

!contains
!include '/work/users/hiepn/nampq/src/sub_Fortran/sub_write_nc3d.f90'
End