initial_grid.f90

Go to the documentation of this file.

subroutine initial_grid(fr)
!--------------------------------------------------------------------
! on exit xx and yy are filled with the values of x, xr, xs, xrs and
! y, yr, ys, yrs of every node :
! xx(1,node) = x, xx(2,node) = xr, xx(3,node) = xs, xx(4,node) = xrs
! the shape of the boundary is given by fr(m)
!--------------------------------------------------------------------
  
  use itm_types
  use mod_parameter
  use mod_dat
  use mod_gauss
  use mod_mesh
  use mod_meshacc
  use mod_output
  use plot_data
  
  implicit none
  
  real(r8), dimension(mfm + 2), intent(in) :: fr
  
  real(r8), dimension(nr) :: sg_equi, density
  real(r8), dimension(np) :: theta, dtc
  real(r8) :: dt, dr
  real(r8) :: thtj, radius
  real(r8) :: rm, drm, drmt, drmtr
  integer(itm_i4) :: i, j, m
  integer(itm_i4) :: node
  
  dt = (1._r8 + ias) * pi / dble(np - 1)
  dr = 1._r8 / dble(nr - 1)
  
!--------------------------- change theta grid to constant arclength
  if (iarc == 1) then
    call arclength(fr, mfm, theta, dtc, xgauss, wgauss)
  else
    do j = 1, np
      theta(j) = dt * (j - 1)
      dtc(j) = 1._r8
    end do
  end if

!TODO: mesh accumulation in psi or s for all imesh possible
!TODO: allow for mesh accumulation in s or in psi in general
  if (imesh == 2 .or. imesh == 3) then
    call mesh_accumulation(nr, n_acc_points, s_acc(1 : n_acc_points), &
     sig(1 : n_acc_points), weights(1 : n_acc_points), equidistant, &
     sg, dsg, ddsg, .false.)
  else if (imesh == 3) then
    call mesh_accumulation(nr, n_acc_points, s_acc(1 : n_acc_points), &
     sig(1 : n_acc_points), weights(1 : n_acc_points), equidistant, &
     sg, dsg, ddsg, .true.)
  else
    do i = 1, nr
      sg(i) = dble(i - 1) / dble(nr - 1)
      dsg(i) = 1._r8
      ddsg(i) = 0._r8
    end do
  end if

  do i = 1, nr
    sg_equi(i) = dble(i)
    density(i) = dble(nr - 1) / dsg(i)
  end do
  if (xmgrace_output) then
    call plot_data_1d('line', sg, sg_equi, nr, 'xmgr_radial_points')
    call plot_data_1d('line', sg, density, nr, 'xmgr_s_density')
  end if
  
  do i = 1, nr
    do j = 1, np
      node = np * (i - 1) + j
      thtj = theta(j)
      radius = sg(nr - i + 1)
      xx(1,node) = radius * fr(1) * cos(thtj) / 2._r8
      xx(2,node) = fr(1) * cos(thtj) / 2._r8
      xx(3,node) = -radius * fr(1) * sin(thtj) / 2._r8
      xx(4,node) = -fr(1) * sin(thtj) / 2._r8
      yy(1,node) = radius * fr(1) * sin(thtj) / 2._r8
      yy(2,node) = fr(1) * sin(thtj) / 2._r8
      yy(3,node) = radius * fr(1) * cos(thtj) / 2._r8
      yy(4,node) = fr(1) * cos(thtj) / 2._r8
!---------------------------- keep ellipticity on axis -----------
      do m = 2, mharm
        if (m == 2) then
          rm = radius * (fr(2 * m - 1) * cos((m - 1) * thtj) &
           + fr(2 * m) * sin((m - 1) * thtj))
          drm = (fr(2 * m - 1) * cos((m - 1) * thtj) + fr(2 * m) &
           * sin((m - 1) * thtj))
          drmt = radius * (-fr(2 * m - 1) * (m - 1) * sin((m - 1) &
           * thtj) + fr(2 * m) * (m - 1) * cos((m - 1) * thtj))
          drmtr = (-fr(2 * m - 1) * (m - 1) * sin((m - 1) * thtj) &
           + fr(2 * m) * (m - 1) * cos((m - 1) * thtj))
        else
          rm = radius**(m - 1) * (fr(2 * m - 1) * cos((m - 1) &
           * thtj) + fr(2 * m) * sin((m - 1) * thtj))
          drm = (m - 1) * radius**(m - 2) * (fr(2 * m - 1) &
           * cos((m - 1) * thtj) + fr(2 * m) * sin((m - 1) * thtj))
          drmt = radius**(m - 1) * (-fr(2 * m - 1) * (m - 1) &
           * sin((m - 1) * thtj) + fr(2 * m) * (m - 1) &
           * cos((m - 1) * thtj))
          drmtr = (m - 1) * radius**(m - 2) * (-fr(2 * m - 1) &
           * (m - 1) * sin((m - 1) * thtj) + fr(2 * m) * (m - 1) &
           * cos((m - 1) * thtj))
        end if
        xx(1, node) = xx(1, node) + rm * cos(thtj)
        yy(1, node) = yy(1, node) + rm * sin(thtj)
        xx(2, node) = xx(2, node) + drm * cos(thtj)
        yy(2, node) = yy(2, node) + drm * sin(thtj)
        xx(3, node) = xx(3, node) - rm * sin(thtj) + drmt * cos(thtj)
        yy(3, node) = yy(3, node) + rm * cos(thtj) + drmt * sin(thtj)
        xx(4, node) = xx(4, node) - drm * sin(thtj) + drmtr * cos(thtj)
        yy(4, node) = yy(4, node) + drm * cos(thtj) + drmtr * sin(thtj)
      end do
      xx(2, node) = -xx(2, node) * dr / 2._r8 * dsg(nr - i + 1)
      xx(3, node) = xx(3, node) * dt / 2._r8 * dtc(j)
      xx(4, node) = -xx(4, node) * dr / 2._r8 * dt / 2._r8 * dtc(j) &
       * dsg(nr - i + 1)
      yy(2, node) = -yy(2, node) * dr / 2._r8 * dsg(nr - i + 1)
      yy(3, node) = yy(3, node) * dt / 2._r8 * dtc(j)
      yy(4, node) = -yy(4, node) * dr / 2._r8 * dt / 2._r8 * dtc(j) &
       * dsg(nr - i + 1)
      psi(4 * (node - 1) + 1) = radius**2
      psi(4 * (node - 1) + 2) = -2._r8* radius * dr / 2._r8 &
       * dsg(nr - i + 1)
      psi(4 * (node - 1) + 3) = 0._r8
      psi(4 * (node - 1) + 4) = 0._r8
    end do
  end do
  
  return
end subroutine initial_grid