! 
! Copyright 2013 Shin'ya Murakami <murashin@gfd-dennou.org>
! $Id:$
! 
! 渦領域を取り出す.
! 

program SelectVortexArea
  
  use vortex
  use dc_string, only: toChar
  
  implicit none
  
  integer nx, ny, nkx, nky, s, nt, vn
  real(8) t_start, t_end, time, dt, dkx, dky, sd, dpi, varea
  real(8) r, l, gamma, energy, enst, velc, vorc, totalenst, totalenergy
  type(fft_struct2d) fs
  integer,    allocatable :: ky(:)
  real(8),    allocatable :: x(:), y(:), owc(:,:), str(:,:), vor(:,:), vorf(:,:)
  complex(8), allocatable :: fsp(:,:)
  logical,    allocatable :: lcond(:,:), lcond2(:,:)
  type(PointList), pointer :: pl => null()
  type(PatchedArea), pointer :: pa => null()
  
  call readnml_grid2d(nx, ny, dkx, dky)
  call readnml_wave2d(nkx, nky)
  allocate(x(0:nx-1), y(0:ny-1), ky(0:ny-1))
  allocate(owc(0:nx-1,0:ny-1), str(0:nx-1,0:ny-1), vor(0:nx-1,0:ny-1), & 
       &   vorf(0:nx-1,0:ny-1), lcond(0:nx-1,0:ny-1), lcond2(0:nx-1,0:ny-1))
  allocate(fsp(0:nx/2,0:ny-1))
  call readnml_readtime(t_start, t_end, dt)
  call fft_mkplan_2d(fs, nx, ny)
  call mktbl_ky(ny, ky)
  call set_axis(x, y, nx, ny, dkx, dky)
  dpi = 2.0d0*acos(-1.0d0)
  
  call HistoryCreate( &
       & file='phys-vortexarea4.nc', &
       & title='vortex area using Okubo-Weiss criterion and vorticity SD', &
       & source='f90-spm 20080728 are used', &
       & institution='Shin-ya Murakami', &
       & dims=(/'x','y','t'/), dimsizes=(/nx,ny,0/), &
       & longnames=(/'x-coordinate', 'y-coordinate', 'time        '/), &
       & units=(/'1','1','1'/), origin=real(t_start), interval=real(dt))
  call HistoryPut('x', x)
  call HistoryPut('y', y)
  call HistoryAddVariable( &
    & varname='vor', dims=(/'x','y','t'/), units='1', xtype='double', &
    & longname='vorticity at vortex area')
  call HistoryAddVariable( &
       & varname='vn', dims=(/'t'/), units='1', xtype='integer', &
       & longname='number of vortex')
  call HistoryAddVariable( &
       & varname='area', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean area of each vortex area')
  call HistoryAddVariable( &
       & varname='r', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean radius of vortices')
  call HistoryAddVariable( &
       & varname='l', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean distance between each nearest vortices')
  call HistoryAddVariable( &
       & varname='gamma', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean circulation of each nearest vortices')
  call HistoryAddVariable( &
       & varname='vorc', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean vorticity at the center of the vortices')
  call HistoryAddVariable( &
       & varname='enst', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean enstrophy')  
  call HistoryAddVariable( &
       & varname='energy', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean energy')
  call HistoryAddVariable( &
       & varname='velc', dims=(/'t'/), units='1', xtype='double', &
       & longname='mean velocity at the center of the vortices')
  
  print *, "# time, number of vortices, total vortex area, "
  print *, "# mean radius, mean distance, mean circulation, "
  print *, "# mean energy, mean enstrophy, "
  print *, "# mean vorticity at the center of the vortices, "
  print *, "# mean velocity at the center of the vortices, "
  print *, "# total energy, total enstrophy"
  
  nt = nint((t_end-t_start)/dt)+1
  do s=0,nt-1
     time = dble(s)*dt+t_start
     call HistoryGet("phys-owc2d.nc", 'owc', owc, 't='//toChar(time))
     call HistoryGet("vor.nc", 'vor', vor, 't='//toChar(time))
     
     ! Okubo-Weiss criterionによる判定
     lcond = SetIsLarger2d(-owc, 0.0d0, nx, ny) ! owcの負の部分
     
     ! 渦度の標準偏差による判定
     sd = StandardDeviation2d(vor, nx, ny)
     print *, "time, sd", time, sd
     lcond2 = SetIsAbsLarger2d(vor, sd, nx, ny) ! sdより大きい渦度
     
     ! 共通部分を取り出す
     lcond = lcond .and. lcond2
     
     ! 場を分解
     pa => vc_array2pa(lcond, nx, ny)
     call vc_pa_small_structure_filter(pa, 4)
     call pa_list2array(pa, lcond, nx, ny)
     vorf = SelectTrue2d(vor, lcond, nx, ny) ! vorfにfilterされたvorを置く.
     vn = vc_count_vortex(pa)
     varea = vc_count_gridpoints(pa)/dble(nx*ny)*area(dkx, dky)
     r = vc_mean_radii(pa, nx, ny, dkx, dky)
     gamma = vc_mean_circulation(pa, vor, nx, ny, dkx, dky)
     enst = vc_mean_enstrophy(pa, vor, nx, ny, dkx, dky)
     
     ! 渦の中心位置を求める
     pl => vc_pa2positions_of_center(pa, x, y, nx, ny, dpi/dkx, dpi/dky)
     
     ! 渦の最近傍渦までの平均距離を求める
     if(vn < 50000)then
        l = vc_mean_distance_of_nearest_vortices(pl, dpi/dkx, dpi/dky)
     else
        l = 0.0d0
     end if
     ! 渦の中心の渦度の大きさの平均値を求める
     ! call grid2wave(vor, fsp, nx, ny, fs) ! fspに渦度
     ! vorcenter = vc_function_mean_at_points(pl, fsp, nx, ny, dkx, dky, ky)
     vorc = vc_abs_function_mean_at_points2(pl, vor, x, y, nx, ny, dkx, dky)
     
     ! エネルギーを求める.
     fsp = grid2wave(vor, nx, ny, fs)
     call vor2str(fsp, fsp, nx, ny, dkx, dky, ky)
     str = wave2grid(fsp, nx, ny, fs)
     energy = vc_mean_energy(pa, vor, str, nx, ny, dkx, dky)
     
     ! 総エンストロフィーを求める
     totalenst = enstrophy_p(vor, nx, ny)*area(dkx,dky)
     
     ! 総エネルギーを求める
     fsp = grid2wave(vor, nx, ny, fs)
     call vor2str(fsp, fsp, nx, ny, dkx, dky, ky)
     call str2velx(fsp, fs%c, nx, ny, dky, ky)
     vor = wave2grid(fs%c, nx, ny, fs)  ! vorにu
     call str2vely(fsp, fs%c, nx, ny, dkx)
     fs%r = wave2grid(fs%c, nx, ny, fs) ! fs%rにv
     totalenergy = energy_p(vor, fs%r, nx, ny)*area(dkx,dky)
     
     ! 速さを求める
     vor = sqrt(vor**2+fs%r**2) ! 速さ
     !call grid2wave(vor, fsp, nx, ny, fs) ! fspに速度の大きさ
     velc = vc_function_mean_at_points2(pl, vor, x, y, nx, ny, dkx, dky)
     
     ! 出力 & 後始末
     print *, time, vn, varea, r, l, gamma, energy, enst, vorc, velc, totalenergy, totalenst
     call pa_destroy(pa)
     call pl_destroy(pl)
     call HistoryPut('t', real(time))
     call HistoryPut('vor', vorf)
     call HistoryPut('vn', vn)
     call HistoryPut('area', varea)
     call HistoryPut('gamma', gamma)
     call HistoryPut('energy', energy)
     call HistoryPut('enst', enst)
     call HistoryPut('vorc', vorc)
     call HistoryPut('velc', velc)
     call HistoryPut('r', r)
     call HistoryPut('l', l)
  end do
  call HistoryClose()
  stop
end program SelectVortexArea
