Program cal_freq
! by: hapt
! to: calculate spatial distribution of the TC frequency
Implicit none
Integer, parameter:: iyr=2013
Character*100:: dataname="_ERA"
Integer, parameter:: nstep=100
Integer:: ntc, itc, istep
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
Real, parameter:: rmis=-99
Character*4:: cyr
Character*100:: idir
Character:: ctmp
! create map to make spatial frequency distribution
Real, parameter:: slon=100.0, elon=160
Real, parameter:: slat=0, elat=35
Real, parameter:: cint=2.5
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
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))

Write(cyr,"(I4)") iyr
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(cyr)//"_TCs"//trim(dataname)//".txt",status="unknown")
     write(*,*) trim(cyr)//"_TCs"//trim(dataname)//".txt"
     Read(1,*) ctmp, ntc
     Allocate(tc_lon(ntc,nstep))
     Allocate(tc_lat(ntc,nstep))
     Allocate(tc_wind(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
     tc_wind=rmis
     Do itc=1,ntc
        Read(1,*) ctmp, ctmp
        !write(*,*) ctmp, tc_step(itc)
        Read(1,*) ctmp, tc_step(itc)
        !write(*,*) ctmp, tc_step(itc)
        Do istep=1,tc_step(itc)
           Read(1,*) tc_day(itc,istep), tc_mon(itc,istep), tc_year(itc,istep), tc_hour(itc,istep), tc_lon(itc,istep), tc_lat(itc,istep), tc_wind(itc,istep)
        End do
     End do
Close(1)

tc_freq=0
Do itc=1,ntc
   flag=0
   Do istep=1,tc_step(itc)-1
      write(*,*) tc_lon(itc,istep), tc_lat(itc,istep)
      write(*,*) tc_lon(itc,istep+1), tc_lat(itc,istep+1)
      If (istep.ne.1) then
          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
                         flag(ilon,ilat)=1
                     End if
                 End do
             End if
          End do
      End if
    
      ! ########### 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
          !write(*,*) tc_lon(itc,istep)
          Do while (tmp_lon.le.(gt_lon+tmp_cint)) 
          !Do while (tmp_lon.eq.lt_lon)
         
             itmp=itmp+1
             tmp_lon=lt_lon+(itmp-1)*tmp_cint
             tmp_lat=tmp_lon*slope+intercept
             !Write(*,"2(F6.3,2x)") slope, intercept
             Write(*,*) "abc",itmp, tmp_lon, tmp_lat

             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
          !write(*,*) "abc", tmp_lat
          Do ilon=1,nlon-1
             If (lon(ilon).le.tmp_lon.and.tmp_lon.le.lon(ilon+1)) then    
                  write(*,*) "check", lon(ilon), tmp_lon, lon(ilon+1)
                  itmp=0
                  Do while (tmp_lat.le.(gt_lat+tmp_cint))
                     itmp=itmp+1
                     tmp_lat=lt_lat + (itmp-1)*tmp_cint
                     write(*,*) tmp_lon, tmp_lat
                     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
                         write(*,*) lon(ilon), tmp_lon, lon(ilon+1)
                         If (lon(ilon).le.tmp_lon.and.tmp_lon.le.lon(ilon+1)) then
                             If (flag(ilon,ilat).eq.0) then
                                 write(*,*) "ajm",tmp_lon, tmp_lat
                                 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
End do
write(*,*) nlon, nlat     

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

Open(5,file="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="lat.list",status="unknown")
     Do ilat=1,nlat-1
        Write(6,'(F7.2)') (lat(ilat)+lat(ilat+1))/2
     End do
Close(6)

!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

End