platelet initcond

! Create restart file with initial condition
program InitCond

  use ModDataTypes
  use ModDataStruct
  use ModRbc
  use ModWall
  use ModConf
  use ModData
  use ModIO
  use ModBasicMath
  use ModPostProcess
  ! use MPI

  implicit none

  integer, parameter :: nrbcMax = 128
  integer, parameter :: nwallMax = 4
  type(t_rbc), pointer :: rbc
  type(t_rbc)         :: rbcRef
  type(t_wall), pointer :: wall
  real(WP) :: radEqv, szCell(3)
  integer :: nlat0, nlat1, ii
  real(WP) :: theta, th, rc, xc(3), ztemp
  real(WP) :: xmin, xmax, ymin, ymax, zmin, zmax
  integer :: iz, i, p, l, dealias
  integer, parameter :: ranseed = 161269
  character(CHRLEN) :: fn
  ! real = double, fp
  real :: lengtube, lengspacing, phi, actlen
  real(WP) :: rand(27, 3), minDist, platRad, platDiam
  integer :: j, ierr, half2, index, layer, newiz
  real :: wbcs(2), tubeDiam, layerx(3), layery(3), halflen

  ! Initialize
  call InitMPI
  ! Allocate working arrays
  allocate (rbcs(nrbcMax))
  ! allocate(walls(nwallMax))

  tubeDiam = 22.0/2.82 ! 7.8, rad = 4

  nrbc = 24
  nlat0 = 12
  nlat1 = 6
  dealias = 3
  phi = 70/real(100)
  ! 11.34
  lengtube = 32.0/2.82 ! nrbc/real(phi) !XXLL

  lengspacing = (lengtube - ((2.62/2.82)*9))/9 ! lengtube/Real(nrbc)
  ! print *, "1"
  nwall = 1
  allocate (walls(nwall))
  wall => walls(1)

  call ReadWallMesh('Input/new_cyl_D6_L13_33.e', wall)
  actlen = 13.33
  wall%f = 0.
  do i = 1, wall%nvert
    th = ATAN2(wall%x(i, 1), wall%x(i, 2))
    wall%x(i, 1) = (tubeDiam/2.0)*COS(th)    !!!!!!!!!
    wall%x(i, 2) = (tubeDiam/2.0)*SIN(th)    !!!!!!!!!
    wall%x(i, 3) = lengtube/actlen*wall%x(i, 3)
  end do

  xmin = minval(wall%x(:, 1))
  xmax = maxval(wall%x(:, 1))

  ymin = minval(wall%x(:, 2))
  ymax = maxval(wall%x(:, 2))

  zmin = minval(wall%x(:, 3))
  zmax = maxval(wall%x(:, 3))

  ! size of the periodic box
  Lb(1) = xmax - xmin + 0.5
  Lb(2) = Lb(1)
  Lb(3) = zmax - zmin
  lengtube = Lb(3)

  print *, "Lb", Lb(:)
  !  Lb   8.3013429849047267        8.3013429849047267        11.347517730496456
  ! lengspacing = lengtube/real(nrbc)
  print *, 'lengtube 2', lengtube
  print *, 'lengspacing 2', lengspacing

  ! Reference cell i
  xc = 0.

  radEqv = 1.0
  platRad = .3

  print *, "nlat1", nlat1, "platRad", platRad

  call Rbc_Create(rbcRef, nlat1, dealias)
  call Rbc_MakePlatelet(rbcRef, platRad, xc)

  do ii = 1, 3
    ! dimensions of cell !
    szCell(ii) = maxval(rbcRef%x(:, :, ii)) - minval(rbcRef%x(:, :, ii))
  end do
  if (rootWorld) then
    print *, "szCell", szCell
  end if

  layerx = (/-1.63, 1.63, 0.0/)
  layery = (/-.943, -.943, 1.89/)

  ! place cells
  do iz = 1, nrbc
    index = (modulo(iz, 3) + 1)
    layer = (iz - 1)/3
    xc(1) = layerx(index) ! + ((rand(iz, 1) - 0.2) - 0.4)
    xc(2) = layery(index) ! + ((rand(iz, 2) - 0.2) - 0.4)
    xc(3) = (layer + .5) + (.5*layer) ! + ((rand(iz, 3) - 0.5)/3)
    if (iz == 20) then
      print *, 'platelet iz:', iz, 'index: ', index, "layer: ", layer, 'xc:', xc
      rbc => rbcs(iz)
      rbc%celltype = 3
      call Rbc_Create(rbc, nlat1, dealias)
      call Rbc_MakePlatelet(rbc, platRad, xc)
    else
      print *, 'iz:', iz, 'index: ', index, "layer: ", layer, 'xc:', xc
      rbc => rbcs(iz)
      rbc%celltype = 1
      call Rbc_Create(rbc, nlat0, dealias)
      call Rbc_MakeBiConcave(rbc, radEqv, xc)
    end if
  end do

  ! Put things in the middle of the periodic box
  call Recenter_Cells_and_Walls

  ! platDiam = GetDiameter(3)

  ! Output
  write (*, '(A,3F10.3)') 'Periodic domain size = ', Lb

  Nt0 = 0; time = 0.
  vBkg(1:2) = 0.; vBkg(3) = 8.

  ! Write intial conditions
  if (nrbc > 0) then
    write (fn, FMT=fn_FMT) 'D/', 'x', 0, '.dat'
    call WriteManyRBCs(fn, nrbc, rbcs)

    write(fn, FMT=fn_FMT) 'D/', '3x', 0, '.dat'
    call WriteManyRBCsByType(fn, nrbc, rbcs, 3)

    write (*, '(A,A)') 'Cell file: ', trim(fn)

    fn = 'D/restart.LATEST.dat'
    call WriteRestart(fn, Nt0, time)
    write (*, '(A,A)') 'Binary restart file: ', trim(fn)
  end if

  if (nwall > 0) then
    write (fn, FMT=fn_FMT) 'D/', 'wall', 0, '.dat'
    call WriteManyWalls(fn, nwall, walls)
    write (*, '(A,A)') 'Wall file: ', trim(fn)
  end if

  write (fn, FMT=fn_FMT) 'D/', 'restart', Nt0, '.dat'
  call WriteRestart(fn, Nt0, time)
  write (*, '(A,A)') 'Binary restart file: ', trim(fn)

  ! Deallocate working arrays
  deallocate (rbcs)

  ! Finalize
  call FinalizeMPI

  stop

contains
  subroutine Recenter_Cells_and_Walls
    integer :: irbc, iwall, ii, gen, repeat, j
    real :: x, y, a
    real, parameter :: PI = 3.14159265359
    type(t_rbc), pointer :: rbc
    type(t_wall), pointer :: wall
    real(WP) :: xmin(3), xmax(3), xc(3), xmax2, ymax2

    ! cells
    ! translate cells
    do irbc = 1, nrbc
      rbc => rbcs(irbc)
      rbc%x(:, :, 1) = rbc%x(:, :, 1) + 0.5*Lb(1)
      rbc%x(:, :, 2) = rbc%x(:, :, 2) + 0.5*Lb(2)
      print *, 'Zc', sum(rbc%x(:, :, 3))/real(rbc%nlat*rbc%nlon)
    end do

    ! walls
    ! find the bounding box
    xmin = 1.D10
    xmax = -1.D10
    do iwall = 1, nwall
      print *, "iwall", iwall
      wall => walls(iwall)
      do ii = 1, 3
        xmin(ii) = min(xmin(ii), minval(wall%x(:, ii)))
        xmax(ii) = max(xmax(ii), maxval(wall%x(:, ii)))
      end do
    end do
    xc = 0.5*(xmin + xmax)

    ! translate walls
    do iwall = 1, nwall
      wall => walls(iwall)
      do ii = 1, 3
        wall%x(:, ii) = wall%x(:, ii) + 0.5*Lb(ii) - xc(ii)
      end do
    end do

    xmax2 = maxval(wall%x(:, 1))
    print *, "xmax2: ", xmax2

    ymax2 = maxval(wall%x(:, 2))
    print *, "ymax2: ", ymax2

  end subroutine Recenter_Cells_and_Walls

end program InitCond