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

Integer, parameter:: iyr=2013
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
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=40
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
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))

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(oDir)//trim(cfile)//".txt",status="unknown")
     !write(*,*) trim(cfile)//".txt"
     Read(1,*) ctmp, ntc
     !Write(*,*) 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
        !Write(*,*) ctmp
        Read(1,*) ctmp, tc_step(itc)
        Read(1,*) ctmp, tc_step(itc)
        !write(*,*) ctmp, tc_step(itc)
        Do istep=1,tc_step(itc)
           !Write(*,*) itc, istep, 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)
           !Write(*,*) tc_day(itc,istep), tc_mon(itc,istep), tc_year(itc,istep), tc_hour(itc,istep), tc_lon(itc,istep), tc_lat(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=trim(oDir)//"TC_freq_"//trim(cfile)//".txt",status="unknown")
     Do ilon=1,nlon-1
        Write(4,'(16(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,*) trim(oDir)//"TC_freq_"//trim(cfile)
Close(7)

! 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