module PanicleClass
   use, intrinsic :: iso_fortran_env
   use KinematicClass
   use StemClass
   use FEMDomainClass

   implicit none

   type  :: panicle_config_
      integer(int32) :: num_seed_column = 30
      real(real64) :: panicle_seed_interval = 1.0d0/1000.0d0
      real(real64) :: panicle_seed_diameter = 5.0d0/1000.0d0
      real(real64) :: panicle_seed_length = 10.0d0/1000.0d0
      real(real64) :: panicle_panicle_diameter = 1.0d0/1000.0d0
   end type

   type :: Panicle_
      type(FEMDomain_)    ::  FEMDomain
      real(real64)        :: Length, Width, Angle
      type(Stem_), pointer ::  pStem
      integer(int32) :: division(1:3) = [5, 5, 5] ! for maize
      integer(int32) :: rice_seed_division(1:3) = [3, 3, 3] ! for rice
      integer(int32) :: wheat_seed_division(1:3) = [2, 2, 2] ! for rice

      integer(int32), allocatable  :: I_planeNodeID(:)
      integer(int32), allocatable  :: I_planeElementID(:)
      integer(int32), allocatable  :: II_planeNodeID(:)
      integer(int32), allocatable  :: II_planeElementID(:)
      integer(int32)  :: A_PointNodeID
      integer(int32)  :: B_PointNodeID
      integer(int32)  :: C_PointNodeID
      integer(int32)  :: D_PointNodeID

      integer(int32)  :: A_PointElementID
      integer(int32)  :: B_PointElementID

      real(real64) :: default_rice_seed_interval = 3.0d0/1000.0d0 ! 3 mm
      real(real64) :: default_rice_seed_branch_length = 3.0d0/1000.0d0 ! 2 mm
      real(real64) :: default_rice_seed_length = 6.0d0/1000.0d0 ! 2 mm
      real(real64) :: default_rice_seed_width = 4.0d0/1000.0d0 ! 2 mm
      real(real64) :: default_rice_seed_thickness = 2.0d0/1000.0d0 ! 2 mm
      real(real64) :: default_rice_panicle_curvature = 0.20d0

      real(real64) :: disp_x
      real(real64) :: disp_y
      real(real64) :: disp_z

      ! For deformation analysis
      real(real64), allocatable :: YoungModulus(:)! element-wise
      real(real64), allocatable :: PoissonRatio(:)! element-wise
      real(real64), allocatable :: Density(:)     ! element-wise
      real(real64), allocatable :: Stress(:, :, :)     ! Gauss point-wise
      real(real64), allocatable :: Displacement(:, :) ! node-wise, three dimensional

   contains
      procedure, public :: Init => initPanicle
      procedure, public :: move => movePanicle
      procedure, public :: rotate => rotatePanicle
      procedure, public :: getCoordinate => getCoordinatePanicle
      procedure, public :: connect => connectPanicle
      procedure, public :: vtk => vtkPanicle
      procedure, public :: stl => stlPanicle
      procedure, public :: ply => plyPanicle
      procedure, public :: remove => removePanicle
   end type

   interface operator(//)
      module procedure append_panicle_object_vector
   end interface

contains

   function to_panicle_config(num_seed_column, panicle_seed_interval, &
                              panicle_seed_diameter, panicle_seed_length, panicle_panicle_diameter) result(this)
      type(panicle_config_) :: this

      integer(int32), intent(in) :: num_seed_column
      real(real64), intent(in)  :: panicle_seed_interval
      real(real64), intent(in)  :: panicle_seed_diameter
      real(real64), intent(in)  :: panicle_seed_length
      real(real64), intent(in)  :: panicle_panicle_diameter

      this%num_seed_column = num_seed_column
      this%panicle_seed_interval = panicle_seed_interval
      this%panicle_seed_diameter = panicle_seed_diameter
      this%panicle_seed_length = panicle_seed_length
      this%panicle_panicle_diameter = panicle_panicle_diameter

   end function

! #####################################################
   recursive subroutine initPanicle(this, Length, Width, Node, shape_factor, debug, x_num, y_num, z_num, rice, &
                                    rice_seed_interval, rice_seed_branch_length, &
                                    rice_seed_length, rice_seed_width, rice_seed_thickness, &
                                    rice_panicle_curvature, rice_seed_division, &
                                    wheat, panicle_config, &
                                    Arabidopsis)
      class(Panicle_), intent(inout) :: this
      type(panicle_config_), optional, intent(in) :: panicle_config ! for wheat

      real(real64), intent(in) :: Length, width ! need for all panicle type

      integer(int32), optional, intent(in) :: Node ! for maize
      integer(int32), optional, intent(in) :: x_num, y_num, z_num ! for all
      integer(int32), optional, intent(in) :: rice_seed_division(1:3)
      real(real64), optional, intent(in) :: rice_seed_interval, rice_seed_branch_length, & ! for rice
                                            rice_seed_length, rice_seed_width, rice_seed_thickness, &
                                            rice_panicle_curvature
      real(real64), optional, intent(in) :: shape_factor ! only for Maize
      !integer(int32),optional,intent(in) :: rice_panicle_branch_num  ! for rice

      logical, optional, intent(in) :: Arabidopsis ! for Arabidopsis

      real(real64):: Angle
      type(Math_) :: math
      type(Random_) :: random
      real(real64) :: dist_val
      real(real64), allocatable :: x_axis(:), y_axis(:), z_axis(:), z_axis0(:)
      integer(int32) :: i, j
      integer(int32), allocatable :: buf(:), kill_element_list(:)
      logical, optional, intent(in) :: debug, rice, wheat

      type(FEMDomain_) :: femdomain
      type(FEMDomain_) :: seed, this_seed
      type(Panicle_), allocatable    :: branch_panicle(:)
      type(FEMDomain_), allocatable :: seeds(:), panicles(:)

      real(real64) :: center_coord(1:3)
      real(real64) :: shape_factor_val

      integer(int32), allocatable :: seed_joint(:)

      real(real64), allocatable :: nodcoord(:, :)
      integer(int32), allocatable :: elemnod(:, :)

      real(real64) :: seed_interval, seed_branch_length, panicle_curvature

      real(real64) :: ln, normal_vector(1:3), x1(1:3), x2(1:3)
      integer(int32) :: elem_idx, n1, n2, n3, n4, num_seed, node_list(1:4)
      integer(int32) :: case_id, nx, ny, nz
      real(real64) :: a(2, 3), b(2, 3), x, y, z, center(1:3), alpha, theta, &
                      seed_length, seed_width, seed_thickness, dx, dy, dz, r, z_max

      logical :: rice_mode, Arabidopsis_mode, wheat_mode

      rice_mode = input(default=.false., option=rice)
      wheat_mode = input(default=.false., option=wheat)
      Arabidopsis_mode = input(default=.false., option=Arabidopsis)

      if (rice_mode) then
         seed_interval = input(default=this%default_rice_seed_interval, &
                               option=rice_seed_interval)
         seed_branch_length = input(default=this%default_rice_seed_branch_length, &
                                    option=rice_seed_branch_length)
         seed_length = input(default=this%default_rice_seed_length, &
                             option=rice_seed_length)
         seed_width = input(default=this%default_rice_seed_width, &
                            option=rice_seed_width)
         seed_thickness = input(default=this%default_rice_seed_thickness, &
                                option=rice_seed_thickness)
         panicle_curvature = input(default=this%default_rice_panicle_curvature, &
                                   option=rice_panicle_curvature)

         if (present(rice_seed_division)) then
            this%rice_seed_division = rice_seed_division
         end if
         nx = this%rice_seed_division(1)
         ny = this%rice_seed_division(2)
         nz = this%rice_seed_division(3)

         call femdomain%create(meshtype="Cube3D", x_num=x_num, y_num=y_num, z_num=z_num)
         call femdomain%resize(x=Width, y=width, z=Length)

         call seed%create(meshtype="Cube3D", x_num=nx, y_num=ny, z_num=nz)

         call seed%resize(x=seed_width, y=seed_thickness, z=seed_length)
         center = seed%centerPosition()
         call seed%move(x=-center(1), y=-center(2), z=-seed%zmin())
         do i = 1, seed%nn()
            x = seed%mesh%nodcoord(i, 1)
            y = seed%mesh%nodcoord(i, 2)
            z = seed%mesh%nodcoord(i, 3)
            alpha = z/seed%zmax()
            theta = 1.00d0 - abs(alpha - 0.50d0)
            seed%mesh%nodcoord(i, 1) = x*theta
         end do

         seed_joint = zeros(0)
         num_seed = int(Length/seed_interval)

         nodcoord = zeros(num_seed*4, 3)
         elemnod = int(zeros(num_seed, FEMDomain%nne()))
         allocate (seeds(num_seed))
         do i = 1, num_seed
            do
               elem_idx = random%randint(From=1, To=femdomain%ne())

               n1 = FEMDomain%mesh%elemnod(elem_idx, 1)
               n2 = FEMDomain%mesh%elemnod(elem_idx, 2)
               n3 = FEMDomain%mesh%elemnod(elem_idx, 3)
               n4 = FEMDomain%mesh%elemnod(elem_idx, 4)
               x1 = FEMDomain%mesh%nodcoord(n2, :) - FEMDomain%mesh%nodcoord(n1, :)
               x2 = FEMDomain%mesh%nodcoord(n4, :) - FEMDomain%mesh%nodcoord(n1, :)

               case_id = 0
               if (minval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 1)) == FEMDomain%xmin()) then
                  if (minval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 2)) == FEMDomain%ymin()) then
                     case_id = 1
                  elseif (maxval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 2)) == FEMDomain%ymax()) then
                     case_id = 2
                  end if

               elseif (maxval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 1)) == FEMDomain%xmax()) then
                  if (minval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 2)) == FEMDomain%ymin()) then
                     case_id = 4
                  elseif (maxval(FEMDomain%mesh%nodcoord([n1, n2, n3, n4], 2)) == FEMDomain%ymax()) then
                     case_id = 3
                  end if
               end if
               if (case_id /= 0) exit
            end do
            seed_joint = seed_joint//[elem_idx]
            select case (case_id)
            case (1)
               !x-min and y_min
               node_list = [4, 1, 5, 8]
            case (2)
               !x-min and y_max
               node_list = [3, 4, 8, 7]
            case (4)
               !x-max and y_min
               node_list = [1, 2, 6, 5]
            case (3)
               !x-max and y_max
               node_list = [2, 3, 7, 6]
            end select

            n1 = FEMDomain%mesh%elemnod(elem_idx, node_list(1))
            n2 = FEMDomain%mesh%elemnod(elem_idx, node_list(2))
            n3 = FEMDomain%mesh%elemnod(elem_idx, node_list(3))
            n4 = FEMDomain%mesh%elemnod(elem_idx, node_list(4))

            x1 = FEMDomain%mesh%nodcoord(n2, :) - FEMDomain%mesh%nodcoord(n1, :)
            x2 = FEMDomain%mesh%nodcoord(n4, :) - FEMDomain%mesh%nodcoord(n1, :)
            normal_vector = cross_product(x1, x2)
            normal_vector = normal_vector/norm(normal_vector)
            nodcoord(4*i - 3, 1:3) = seed_branch_length*normal_vector &
                                     + FEMDomain%mesh%nodcoord(FEMDomain%mesh%elemnod(elem_idx, node_list(1)), :)
            nodcoord(4*i - 2, 1:3) = seed_branch_length*normal_vector &
                                     + FEMDomain%mesh%nodcoord(FEMDomain%mesh%elemnod(elem_idx, node_list(2)), :)
            nodcoord(4*i - 1, 1:3) = seed_branch_length*normal_vector &
                                     + FEMDomain%mesh%nodcoord(FEMDomain%mesh%elemnod(elem_idx, node_list(3)), :)
            nodcoord(4*i - 0, 1:3) = seed_branch_length*normal_vector &
                                     + FEMDomain%mesh%nodcoord(FEMDomain%mesh%elemnod(elem_idx, node_list(4)), :)
            elemnod(i, 1:4) = FEMDomain%mesh%elemnod(elem_idx, node_list(1:4))
            elemnod(i, 5:8) = [4*i - 3, 4*i - 2, 4*i - 1, 4*i - 0] + femdomain%nn()

            ! 果梗と子実をがっちゃんこ
            ! seedの5要素目，1234番
            seeds(i) = seed
            dx = sum(nodcoord(4*i - 3:4*i, 1))/4.0d0
            dy = sum(nodcoord(4*i - 3:4*i, 2))/4.0d0
            dz = sum(nodcoord(4*i - 3:4*i, 3))/4.0d0 + seed_branch_length
            call seeds(i)%move(x=dx, y=dy, z=dz)
            call seeds(i)%rotate(z=math%pi/2.0d0*dble(case_id))

            select case (case_id)
            case (1)
               !x-min and y_min

               nodcoord(4*i - 3, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 3), :)
               nodcoord(4*i - 2, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 4), :)
               nodcoord(4*i - 1, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 1), :)
               nodcoord(4*i - 0, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 2), :)
            case (2)
               !x-min and y_max
               ![ok]
               nodcoord(4*i - 3, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 1), :)
               nodcoord(4*i - 2, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 2), :)
               nodcoord(4*i - 1, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 3), :)
               nodcoord(4*i - 0, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 4), :)
            case (4)
               !x-max and y_min
               ![ok]
               nodcoord(4*i - 3, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 1), :)
               nodcoord(4*i - 2, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 2), :)
               nodcoord(4*i - 1, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 3), :)
               nodcoord(4*i - 0, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 4), :)
            case (3)
               !x-max and y_max
               !
               nodcoord(4*i - 3, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 3), :)
               nodcoord(4*i - 2, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 4), :)
               nodcoord(4*i - 1, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 1), :)
               nodcoord(4*i - 0, 1:3) = seeds(i)%mesh%nodcoord(seeds(i)%mesh%elemnod(int(dble(nx*ny)/2.0d0), 2), :)
            end select

         end do

         FEMDomain%mesh%nodcoord = FEMDomain%mesh%nodcoord//nodcoord
         FEMDomain%mesh%elemnod = FEMDomain%mesh%elemnod//elemnod

         do i = 1, size(seeds)
            seeds(i)%mesh%elemnod = seeds(i)%mesh%elemnod + femdomain%nn()
            femdomain%mesh%nodcoord = femdomain%mesh%nodcoord//seeds(i)%mesh%nodcoord
            femdomain%mesh%elemnod = femdomain%mesh%elemnod//seeds(i)%mesh%elemnod
         end do

         ! remove duplicate nodes
         call femdomain%Deduplicate(error=dble(1.0e-8))

         ! bend
         z_max = femdomain%zmax()
         call femdomain%move(x=panicle_curvature)
         do i = 1, femdomain%nn()
            x = femdomain%mesh%nodcoord(i, 1)
            z = femdomain%mesh%nodcoord(i, 3)
            r = x
            theta = z/panicle_curvature
            femdomain%mesh%nodcoord(i, 1) = r*cos(theta)
            femdomain%mesh%nodcoord(i, 3) = r*sin(theta)

         end do

         call femdomain%move(x=-panicle_curvature)
         theta = random%gauss(mu=0.0d0, sigma=math%pi/4.0d0)
         call femdomain%rotate(z=theta)
         dx = femdomain%mesh%nodcoord(1, 1)
         dy = femdomain%mesh%nodcoord(1, 2)
         dz = femdomain%mesh%nodcoord(1, 3)
         call femdomain%move(x=-dx, y=-dy, z=-dz)
         this%femdomain = femdomain

         this%A_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=0.0d0)
         this%B_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=Length)

         ! branch exists
         !if(present(rice_panicle_branch_num) )then
         !    allocate(branch_panicle(rice_panicle_branch_num) )
         !    ! create branches
         !    do i=1,rice_panicle_branch_num
         !
         !        call branch_panicle(i)%init(&
         !            Length=Length,Width=Width,x_num=x_num,y_num=y_num,z_num=z_num,&
         !            rice=.true.,&
         !            rice_seed_interval= rice_seed_interval ,&
         !            rice_seed_branch_length= rice_seed_branch_length ,&
         !            rice_seed_length= rice_seed_length ,&
         !            rice_seed_width= rice_seed_width ,&
         !            rice_seed_thickness= rice_seed_thickness ,&
         !            rice_panicle_curvature= rice_panicle_curvature ,&
         !            rice_seed_division= rice_seed_division  )
         !
         !    enddo
         !endif

         return

      elseif (wheat_mode) then
         if (.not. present(panicle_config)) then
            print *, "ERROR :: [Panicle%init for wheat], panicle_config should be present. "
            stop
         end if

         femdomain = to_wheat_panicle_mesh( &
                     num_seed_column=30, &
                     panicle_seed_interval=1.0d0/1000.0d0, &
                     panicle_seed_diameter=5.0d0/1000.0d0, &
                     panicle_seed_length=10.0d0/1000.0d0, &
                     panicle_panicle_diameter=1.0d0/1000.0d0 &
                     )
         this%femdomain = femdomain

         this%A_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=0.0d0)
         this%B_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=Length)

         return
      elseif (Arabidopsis_mode) then
         print *, "No Arabidopsis-mode is implemented for Panicle"
         return
      else
         !<<<< MAIZE MODE >>>>>
         shape_factor_val = input(default=0.33d0, option=shape_factor)

         Angle = 0.0d0 ! vertical panicle
         this%Length = length
         this%Width = Width
         this%Angle = Angle
         this%division(1) = input(default=this%division(1), option=x_num)
         this%division(2) = input(default=this%division(2), option=y_num)
         this%division(3) = input(default=this%division(3), option=z_num)

         x_axis = [-Length*shape_factor_val, -width/2.0d0, 0.0d0, width/2.0d0, Length*shape_factor_val]
         call refine(x_axis, this%division(1))

         !y_axis=[-width/2.0d0,0.0d0,width/2.0d0]
         !call refine(y_axis,this%division(2))
         ! debug
         y_axis = [-Length*shape_factor_val, -width/2.0d0, 0.0d0, width/2.0d0, Length*shape_factor_val]
         call refine(y_axis, this%division(2))

         z_axis = [0.0d0]
         do i = 1, Node
            z_axis = z_axis//[this%Length*shape_factor_val/dble(Node)*dble(i)]
            z_axis = z_axis//[z_axis(size(z_axis)) + this%width]
         end do

         z_axis = z_axis//[this%Length]
         z_axis0 = z_axis
         call refine(z_axis, this%division(3))

         call this%FEMDomain%create("Cube3D", &
                                    x_axis=x_axis, y_axis=y_axis, z_axis=z_axis)
         kill_element_list = zeros(this%FEMDomain%ne())
         !do i=1,this%FEMDomain%ne()
         !    center_coord = this%FEMDomain%centerPosition(ElementID=i)
         !    do j=1,size(z_axis0)-1,2
         !        if( z_axis0(j)< center_coord(3) .and. center_coord(3) < z_axis0(j+1) )then
         !            if( abs(center_coord(1)) > width/2.0d0  )then
         !                kill_element_list(i)  = 1
         !            endif
         !        endif
         !    enddo
         !enddo
         do i = 1, this%FEMDomain%ne()
            center_coord = this%FEMDomain%centerPosition(ElementID=i)
            if (abs(center_coord(1)) > width/2.0d0 .and. abs(center_coord(2)) > width/2.0d0) then
               kill_element_list(i) = 1
            end if
            do j = 1, size(z_axis0) - 1, 2

               if (z_axis0(j) < center_coord(3) .and. center_coord(3) < z_axis0(j + 1)) then
                  if (abs(center_coord(1)) > width/2.0d0 .and. abs(center_coord(2)) < width/2.0d0) then
                     kill_element_list(i) = 1
                  end if

                  if (abs(center_coord(1)) < width/2.0d0 .and. abs(center_coord(2)) > width/2.0d0) then
                     kill_element_list(i) = 1
                  end if

               end if
            end do
         end do

         call this%FEMDomain%killElement(blacklist=kill_element_list, flag=1)

         ! edit
         do i = 1, this%FEMDomain%nn()
            center_coord = this%FEMDomain%mesh%nodcoord(i, :)
            if (abs(center_coord(1)) > width/2.0d0) then
               alpha = center_coord(3)/Length
               if (alpha == 1.0d0) cycle
               theta = radian(alpha*90.0d0) ! angle:: alpha=0 => 0, alpha=1 => radian(90.0)

               ! new x
               this%FEMDomain%mesh%nodcoord(i, 1) = this%FEMDomain%mesh%nodcoord(i, 1)*cos(theta)
               ! new z
               this%FEMDomain%mesh%nodcoord(i, 3) = this%FEMDomain%mesh%nodcoord(i, 3) + abs(center_coord(1))*tan(theta)  ! x * tan(theta)

            end if

            if (abs(center_coord(2)) > width/2.0d0) then
               alpha = center_coord(3)/Length
               if (alpha == 1.0d0) cycle
               theta = radian(alpha*90.0d0) ! angle:: alpha=0 => 0, alpha=1 => radian(90.0)

               ! new x
               this%FEMDomain%mesh%nodcoord(i, 2) = this%FEMDomain%mesh%nodcoord(i, 2)*cos(theta)
               ! new z
               this%FEMDomain%mesh%nodcoord(i, 3) = this%FEMDomain%mesh%nodcoord(i, 3) + abs(center_coord(2))*tan(theta)  ! x * tan(theta)

            end if
         end do

         call this%FEMDomain%rotate(z=math%PI*2.0d0*random%random())

         !    ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
         !    this%I_planeNodeID = this%FEMdomain%mesh%getNodeList(zmax=0.0d0)
         !    this%I_planeElementID = this%FEMdomain%mesh%getElementList(zmax=0.0d0)
         !
         !    ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
         !    this%II_planeNodeID = this%FEMdomain%mesh%getNodeList(zmin=this%length)
         !    this%II_planeElementID = this%FEMdomain%mesh%getElementList(zmin=this%length)
         !
         !
         !
         !    center_coord(1) = sum(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),1) )&
         !        /size(this%I_planeNodeID)
         !    center_coord(2) = sum(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),2) )&
         !        /size(this%I_planeNodeID)
         !    center_coord(3) = sum(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),3) )&
         !        /size(this%I_planeNodeID)
         !
         !    dist_val = norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(1),:)-center_coord)
         !    this%A_PointNodeID = this%I_planeNodeID(1)
         !
         !    do i=2, size(this%I_planeNodeID)
         !        if(  norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(i),:)-center_coord) < dist_val  )then
         !            this%A_PointNodeID = this%I_planeNodeID(i)
         !            dist_val = norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(i),:)-center_coord)
         !        endif
         !    enddo
         !
         !
         !    center_coord(1) = sum(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),1) )&
         !        /size(this%I_planeNodeID)
         !    center_coord(2) = sum(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),2) )&
         !        /size(this%I_planeNodeID)
         !    center_coord(3) = sum(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),3) )&
         !        /size(this%I_planeNodeID)
         !
         !    dist_val = norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(1),:)-center_coord)
         !    this%B_PointNodeID = this%II_planeNodeID(1)
         !
         !    do i=2, size(this%II_planeNodeID)
         !        if(  norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(i),:)-center_coord) < dist_val  )then
         !            this%B_PointNodeID = this%II_planeNodeID(i)
         !            dist_val = norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(i),:)-center_coord)
         !        endif
         !    enddo
         !
         !
         !
         !    center_coord(1) = maxval(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),1) )
         !
         !    center_coord(2) = sum(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),2) )&
         !        /size(this%I_planeNodeID)
         !
         !    center_coord(3) = sum(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(:),3) )&
         !        /size(this%I_planeNodeID)
         !
         !    dist_val = norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(1),:)-center_coord)
         !    this%C_PointNodeID = this%I_planeNodeID(1)
         !
         !    do i=2, size(this%I_planeNodeID)
         !        if(  norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(i),:)-center_coord) < dist_val  )then
         !            this%C_PointNodeID = this%I_planeNodeID(i)
         !            dist_val = norm(this%FEMDomain%mesh%nodcoord(this%I_planeNodeID(i),:)-center_coord)
         !        endif
         !    enddo
         !
         !
         !
         !    center_coord(1) = sum(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),1) )&
         !        /size(this%II_planeNodeID)
         !
         !    center_coord(2) = maxval(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),2) )
         !
         !
         !    center_coord(3) = sum(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(:),3) )&
         !        /size(this%II_planeNodeID)
         !
         !    dist_val = norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(1),:)-center_coord)
         !    this%D_PointNodeID = this%II_planeNodeID(1)
         !
         !    do i=2, size(this%II_planeNodeID)
         !        if(  norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(i),:)-center_coord) < dist_val  )then
         !            this%D_PointNodeID = this%II_planeNodeID(i)
         !            dist_val = norm(this%FEMDomain%mesh%nodcoord(this%II_planeNodeID(i),:)-center_coord)
         !        endif
         !    enddo
         !
         !    buf    = this%FEMDomain%mesh%getElementList(&
         !        xmin=this%width/2.0d0 - this%width/dble(this%division(1) )/2.0d0 ,&
         !        xmax=this%width/2.0d0 + this%width/dble(this%division(1) )/2.0d0 ,&
         !        ymin=this%width/2.0d0 - this%width/dble(this%division(2) )/2.0d0 ,&
         !        ymax=this%width/2.0d0 + this%width/dble(this%division(2) )/2.0d0 ,&
         !        zmax=0.0d0)
         !    !this%A_PointElementID = buf(1)
         !    this%A_PointElementID = median(buf)
         !
         !
         !    buf    = this%FEMDomain%mesh%getElementList(&
         !        xmin=this%width/2.0d0 - this%width/dble(this%division(1) )/2.0d0 ,&
         !        xmax=this%width/2.0d0 + this%width/dble(this%division(1) )/2.0d0 ,&
         !        ymin=this%width/2.0d0 - this%width/dble(this%division(2) )/2.0d0 ,&
         !        ymax=this%width/2.0d0 + this%width/dble(this%division(2) )/2.0d0 ,&
         !        zmin=this%length)
         !
         !    !this%B_PointElementID = buf(1)
         !    this%B_PointElementID = median(buf)
         !
         !    if(debug) print *, this%A_PointNodeID
         !    if(debug) print *, this%B_PointNodeID
         !    if(debug) print *, this%A_PointElementID
         !    if(debug) print *, this%B_PointElementID
         !
         this%A_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=0.0d0)
         this%B_PointNodeID = this%FEMDomain%getNearestNodeID(x=0.0d0, y=0.0d0, z=Length)
      end if
   end subroutine
! #####################################################

! #####################################################
   subroutine vtkPanicle(this, name)
      class(Panicle_), intent(inout) :: this
      character(*), intent(in) :: name

      call this%FEMDomain%vtk(name=name)

   end subroutine
! #####################################################

! ########################################
   recursive subroutine movePanicle(this, x, y, z, reset)
      class(Panicle_), intent(inout) :: this
      real(real64), optional, intent(in) :: x, y, z
      logical, optional, intent(in) :: reset
      real(real64), allocatable :: origin1(:), origin2(:), disp(:)

      if (present(reset)) then
         if (reset .eqv. .true.) then
            call this%femdomain%move(-this%disp_x, -this%disp_y, -this%disp_z)
            this%disp_x = 0.0d0
            this%disp_y = 0.0d0
            this%disp_z = 0.0d0
         end if
      end if

      call this%femdomain%move(x, y, z)
      this%disp_x = this%disp_x + input(default=0.0d0, option=x)
      this%disp_y = this%disp_y + input(default=0.0d0, option=y)
      this%disp_z = this%disp_z + input(default=0.0d0, option=z)

   end subroutine
! ########################################

! ########################################
   function getCoordinatePanicle(this, nodetype) result(ret)
      class(Panicle_), intent(in) :: this
      character(*), intent(in) :: nodetype
      real(real64), allocatable :: ret(:)
      integer(int32) :: dimnum, n, i
      integer(int32), allocatable :: buf(:)

      dimnum = size(this%femdomain%mesh%nodcoord, 2)

      allocate (ret(dimnum))
      ret(:) = 0.0d0
      if (nodetype == "A" .or. nodetype == "a") then

         ! 20220701 this may be correct
         ret = this%femdomain%mesh%nodcoord(this%A_PointNodeID, :)
         return

         n = size(this%I_planeNodeID)
         if (n == 0) then
            print *, "ERROR >> getCoordinatePanicle >> size(this%I_planeNodeID) = 0"
         end if
         if (.not. allocated(this%I_planeNodeID)) then

            print *, "ERROR >> getCoordinatePanicle >> .not. allocated(this%I_planeNodeID) "

         end if
         do i = 1, n
            ret(:) = ret(:) + this%femdomain%mesh%nodcoord(this%I_planeNodeID(i), :)
         end do
         ret(:) = 1.0d0/dble(n)*ret(:)

         !ret = this%femdomain%mesh%nodcoord(this%A_PointNodeID,:)
      end if

      if (nodetype == "B" .or. nodetype == "b") then
         !ret = this%femdomain%mesh%nodcoord(this%B_PointNodeID,:)

         ! 20220701 this may be correct
         ret = this%femdomain%mesh%nodcoord(this%B_PointNodeID, :)
         return

         n = size(this%II_planeNodeID)
         if (n == 0) then
            print *, "ERROR >> getCoordinatePanicle >> size(this%II_planeNodeID) = 0"
         end if
         if (.not. allocated(this%I_planeNodeID)) then

            print *, "ERROR >> getCoordinatePanicle >> .not. allocated(this%II_planeNodeID) "

         end if
         do i = 1, n
            ret(:) = ret(:) + this%femdomain%mesh%nodcoord(this%II_planeNodeID(i), :)
         end do
         ret(:) = 1.0d0/dble(n)*ret(:)
      end if

   end function
! ########################################

! ########################################

   subroutine connectPanicle(obj, direct, stem)
      class(Panicle_), intent(inout) :: obj
      class(Stem_), intent(inout) :: stem

      character(2), intent(in) :: direct
      real(real64), allocatable :: x1(:), x2(:), disp(:)

      if (direct == "->" .or. direct == "=>") then
         ! move obj to connect stem (stem is not moved.)
         x1 = obj%getCoordinate("A")
         x2 = stem%getCoordinate("B")
         disp = x2 - x1
         call obj%move(x=disp(1), y=disp(2), z=disp(3))
      end if

      if (direct == "<-" .or. direct == "<=") then
         ! move obj to connect stem (stem is not moved.)
         x1 = stem%getCoordinate("A")
         x2 = obj%getCoordinate("B")
         disp = x2 - x1
         call stem%move(x=disp(1), y=disp(2), z=disp(3))
      end if
   end subroutine

! ########################################

! ########################################
   recursive subroutine rotatePanicle(this, x, y, z)
      class(Panicle_), intent(inout) :: this
      real(real64), optional, intent(in) :: x, y, z

      call this%FEMDomain%rotate(x=x, y=y, z=z)
   end subroutine
! ########################################

! ##############################################

   subroutine stlPanicle(obj, name)
      class(Panicle_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%stl(Name=name)
   end subroutine
! ########################################

   subroutine plyPanicle(obj, name)
      class(Panicle_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%ply(Name=name)
   end subroutine
! ########################################
   subroutine removePanicle(this)
      class(Panicle_), intent(inout) :: this

      call this%FEMDomain%remove()
      this%Length = 0.0d0
      this%Width = 0.0d0
      this%Angle = 0.0d0
      if (associated(this%pStem)) nullify (this%pStem)
      this%division(1:3) = [5, 5, 5]

      if (allocated(this%I_planeNodeID)) deallocate (this%I_planeNodeID)! (:)
      if (allocated(this%I_planeElementID)) deallocate (this%I_planeElementID)! (:)
      if (allocated(this%II_planeNodeID)) deallocate (this%II_planeNodeID)! (:)
      if (allocated(this%II_planeElementID)) deallocate (this%II_planeElementID)! (:)
      this%A_PointNodeID = 0
      this%B_PointNodeID = 0
      this%C_PointNodeID = 0
      this%D_PointNodeID = 0

      this%A_PointElementID = 0
      this%B_PointElementID = 0

      this%default_rice_seed_interval = 3.0d0/1000.0d0 ! 3 mm
      this%default_rice_seed_branch_length = 3.0d0/1000.0d0 ! 2 mm
      this%default_rice_seed_length = 6.0d0/1000.0d0 ! 2 mm
      this%default_rice_seed_width = 4.0d0/1000.0d0 ! 2 mm
      this%default_rice_seed_thickness = 2.0d0/1000.0d0 ! 2 mm
      this%default_rice_panicle_curvature = 0.20d0

      this%disp_x = 0.0d0
      this%disp_y = 0.0d0
      this%disp_z = 0.0d0

      ! For deformation analysis
      if (allocated(this%YoungModulus)) deallocate (this%YoungModulus)! (:)! element-wise
      if (allocated(this%PoissonRatio)) deallocate (this%PoissonRatio)! (:)! element-wise
      if (allocated(this%Density)) deallocate (this%Density)! (:)     ! element-wise
      if (allocated(this%Stress)) deallocate (this%Stress)! (:,:,:)     ! Gauss point-wise
      if (allocated(this%Displacement)) deallocate (this%Displacement)! (:,:) ! node-wise, three dimensional

   end subroutine

! ################################################################

   function to_wheat_panicle_mesh(num_seed_column, panicle_seed_interval, panicle_seed_diameter, &
                                  panicle_seed_length, panicle_panicle_diameter, seed_angle, &
                                  culm_length, culm_diameter, culm_division, heights_vs_diameters) result(output_d)
      type(FEMDomain_) :: output_d
      type(Mesh_) :: wheat, culm
      ! args
      integer(int32), intent(in) :: num_seed_column
      real(real64), intent(in) :: panicle_seed_interval
      real(real64), intent(in) :: panicle_seed_diameter
      real(real64), intent(in) :: panicle_seed_length
      real(real64), intent(in) :: panicle_panicle_diameter

      real(real64), optional, intent(in) :: seed_angle
      ! for culm
      real(real64), optional, intent(in) :: culm_diameter
      real(real64), optional, intent(in) :: heights_vs_diameters(:, :)
      real(real64), optional, intent(in) :: culm_length
      integer(int32), optional, intent(in) :: culm_division

      integer(int32) :: i, j, itr, tot_nod_num
      real(real64)   :: center(1:3), dz, dx, dy, z0
      real(real64)   :: dtheta, rotmat_z(3, 3), theta
      real(real64), allocatable :: hvd(:, :)

      if (present(seed_angle)) then
         theta = radian(seed_angle)
      else
         theta = radian(45.0d0)
      end if

      wheat%nodcoord = zeros(num_seed_column*28, 3)
      wheat%elemnod = int(zeros(num_seed_column*6, 8))
      center(:) = 0.0d0

      itr = 0
      do i = 1, num_seed_column
         z0 = center(3)
         ! stem
         wheat%nodcoord(itr + 1, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3)]
         wheat%nodcoord(itr + 2, 1:3) = &
            [panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3)]
         wheat%nodcoord(itr + 3, 1:3) = &
            [panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3)]
         wheat%nodcoord(itr + 4, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3)]

         dz = panicle_seed_interval
         wheat%nodcoord(itr + 5, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3) + dz]
         wheat%nodcoord(itr + 6, 1:3) = &
            [panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3) + dz]
         wheat%nodcoord(itr + 7, 1:3) = &
            [panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3) + dz]
         wheat%nodcoord(itr + 8, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3) + dz]

         ! right seed
         itr = itr + 8

         dz = panicle_seed_diameter*0.50d0
         dy = panicle_seed_diameter*0.50d0

         wheat%nodcoord(itr + 1, 1:3) = &
            [panicle_seed_length*0.50d0*cos(theta), &
             -panicle_panicle_diameter*0.50d0 - dy, &
             panicle_seed_length*0.50d0*sin(theta) - panicle_panicle_diameter*0.50d0 - dz + z0]
         wheat%nodcoord(itr + 2, 1:3) = &
            [panicle_seed_length*0.50d0*cos(theta), &
             -panicle_panicle_diameter*0.50d0 - dy, &
             panicle_seed_length*0.50d0*sin(theta) + panicle_panicle_diameter*0.50d0 + dz + z0]
         wheat%nodcoord(itr + 3, 1:3) = &
            [panicle_seed_length*0.50d0*cos(theta), &
             panicle_panicle_diameter*0.50d0 + dy, &
             panicle_seed_length*0.50d0*sin(theta) + panicle_panicle_diameter*0.50d0 + dz + z0]
         wheat%nodcoord(itr + 4, 1:3) = &
            [panicle_seed_length*0.50d0*cos(theta), &
             panicle_panicle_diameter*0.50d0 + dy, &
             panicle_seed_length*0.50d0*sin(theta) - panicle_panicle_diameter*0.50d0 - dz + z0]

         wheat%nodcoord(itr + 5, 1:3) = &
            [panicle_seed_length*cos(theta), &
             -panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) - panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 6, 1:3) = &
            [panicle_seed_length*cos(theta), &
             -panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) + panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 7, 1:3) = &
            [panicle_seed_length*cos(theta), &
             panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) + panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 8, 1:3) = &
            [panicle_seed_length*cos(theta), &
             panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) - panicle_panicle_diameter*0.50d0 + z0]

         ! stem
         itr = itr + 8
         dz = panicle_seed_interval
         wheat%nodcoord(itr + 1, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3) + 2*dz]
         wheat%nodcoord(itr + 2, 1:3) = &
            [panicle_panicle_diameter*0.50d0, -panicle_panicle_diameter*0.50d0, center(3) + 2*dz]
         wheat%nodcoord(itr + 3, 1:3) = &
            [panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3) + 2*dz]
         wheat%nodcoord(itr + 4, 1:3) = &
            [-panicle_panicle_diameter*0.50d0, panicle_panicle_diameter*0.50d0, center(3) + 2*dz]

         ! left seed
         itr = itr + 4
         z0 = center(3) + panicle_seed_interval

         wheat%nodcoord(itr + 1, 1:3) = &
            [-panicle_seed_length*0.50d0*cos(theta), &
             -panicle_panicle_diameter*0.50d0 - dy, &
             panicle_seed_length*0.50d0*sin(theta) - panicle_panicle_diameter*0.50d0 - dz + z0]
         wheat%nodcoord(itr + 2, 1:3) = &
            [-panicle_seed_length*0.50d0*cos(theta), &
             -panicle_panicle_diameter*0.50d0 - dy, &
             panicle_seed_length*0.50d0*sin(theta) + panicle_panicle_diameter*0.50d0 + dz + z0]
         wheat%nodcoord(itr + 3, 1:3) = &
            [-panicle_seed_length*0.50d0*cos(theta), &
             panicle_panicle_diameter*0.50d0 + dy, &
             panicle_seed_length*0.50d0*sin(theta) + panicle_panicle_diameter*0.50d0 + dz + z0]
         wheat%nodcoord(itr + 4, 1:3) = &
            [-panicle_seed_length*0.50d0*cos(theta), &
             panicle_panicle_diameter*0.50d0 + dy, &
             panicle_seed_length*0.50d0*sin(theta) - panicle_panicle_diameter*0.50d0 - dz + z0]

         wheat%nodcoord(itr + 5, 1:3) = &
            [-panicle_seed_length*cos(theta), &
             -panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) - panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 6, 1:3) = &
            [-panicle_seed_length*cos(theta), &
             -panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) + panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 7, 1:3) = &
            [-panicle_seed_length*cos(theta), &
             panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) + panicle_panicle_diameter*0.50d0 + z0]
         wheat%nodcoord(itr + 8, 1:3) = &
            [-panicle_seed_length*cos(theta), &
             panicle_panicle_diameter*0.50d0, &
             panicle_seed_length*sin(theta) - panicle_panicle_diameter*0.50d0 + z0]
         center(3) = center(3) + 2*dz
         itr = itr + 8

         dtheta = radian(90.0d0)
         rotmat_z(1, 1) = cos(dtheta*(i - 1)); rotmat_z(1, 2) = -sin(dtheta*(i - 1)); rotmat_z(1, 3) = 0.0d0; 
         rotmat_z(2, 1) = sin(dtheta*(i - 1)); rotmat_z(2, 2) = cos(dtheta*(i - 1)); rotmat_z(2, 3) = 0.0d0; 
         rotmat_z(3, 1) = 0.0d0; rotmat_z(3, 2) = 0.0d0; rotmat_z(3, 3) = 1.0d0; 
         do j = itr - 27, itr
            wheat%nodcoord(j, 1:3) = matmul(rotmat_z, wheat%nodcoord(j, 1:3))
         end do
         if (i > 2) then
            wheat%nodcoord(itr - 27, 1:3) = wheat%nodcoord(itr - 27 - 28, 1:3)
            wheat%nodcoord(itr - 26, 1:3) = wheat%nodcoord(itr - 26 - 28, 1:3)
            wheat%nodcoord(itr - 25, 1:3) = wheat%nodcoord(itr - 25 - 28, 1:3)
            wheat%nodcoord(itr - 24, 1:3) = wheat%nodcoord(itr - 24 - 28, 1:3)
         end if
      end do

! connectivity
      itr = 0
      tot_nod_num = 0
      do i = 1, num_seed_column
         ! stem
         wheat%elemnod(itr + 1, :) = [ &
                                     tot_nod_num + 1, tot_nod_num + 2, tot_nod_num + 3, tot_nod_num + 4, &
                                     tot_nod_num + 5, tot_nod_num + 6, tot_nod_num + 7, tot_nod_num + 8]
         ! seed
         wheat%elemnod(itr + 2, :) = [ &
                                     tot_nod_num + 2, tot_nod_num + 8 + 1, tot_nod_num + 8 + 4, tot_nod_num + 3, &
                                     tot_nod_num + 6, tot_nod_num + 8 + 2, tot_nod_num + 8 + 3, tot_nod_num + 7]
         ! seed
         wheat%elemnod(itr + 3, :) = [ &
                                     tot_nod_num + 8 + 1, tot_nod_num + 8 + 4 + 1, tot_nod_num + 8 + 4 + 4, tot_nod_num + 8 + 4, &
                                     tot_nod_num + 8 + 2, tot_nod_num + 8 + 4 + 2, tot_nod_num + 8 + 4 + 3, tot_nod_num + 8 + 3]

         ! stem

         wheat%elemnod(itr + 4, :) = [ &
                                     tot_nod_num + 5, tot_nod_num + 6, tot_nod_num + 7, tot_nod_num + 8, &
                                     tot_nod_num + 4*4 + 1, tot_nod_num + 4*4 + 2, tot_nod_num + 4*4 + 3, tot_nod_num + 4*4 + 4]
         ! seed
         wheat%elemnod(itr + 5, :) = [ &
                                     tot_nod_num + 4*5 + 1, tot_nod_num + 5, tot_nod_num + 8, tot_nod_num + 4*5 + 4, &
                                     tot_nod_num + 4*5 + 2, tot_nod_num + 4*4 + 1, tot_nod_num + 4*4 + 4, tot_nod_num + 4*5 + 3]

         wheat%elemnod(itr + 6, :) = [ &
                                     tot_nod_num + 4*6 + 1, tot_nod_num + 4*5 + 1, tot_nod_num + 4*5 + 4, tot_nod_num + 4*6 + 4, &
                                     tot_nod_num + 4*6 + 2, tot_nod_num + 4*5 + 2, tot_nod_num + 4*5 + 3, tot_nod_num + 4*6 + 3]
         itr = itr + 6
         tot_nod_num = tot_nod_num + 28
      end do
      output_d%mesh = wheat
      call output_d%Deduplicate(error=dble(1.0e-8))

      if (present(culm_length)) then
         if (present(heights_vs_diameters)) then
            ! create culm
            wheat = output_d%mesh
            culm%nodcoord = zeros(culm_division*4, 3)
            culm%elemnod = int(zeros(culm_division, 8))
            itr = 0
            hvd = heights_vs_diameters.n. [culm_length, panicle_panicle_diameter]

            do i = 1, culm_division
               theta = dble(i)/dble(culm_division)

               !dx = culm_diameter*0.50d0*theta &
               !    + panicle_panicle_diameter*0.50d0*(1.0d0-theta)
               !dy = culm_diameter*0.50d0*theta &
               !    + panicle_panicle_diameter*0.50d0*(1.0d0-theta)
               dx = 0.50d0*to_culm_diameter(heights_vs_diameters=hvd, height=(1.0d0 - theta)*culm_length)
               dy = dx

               z0 = -theta*culm_length
               culm%nodcoord(itr + 1, 1:3) = [-dx, -dy, z0]
               culm%nodcoord(itr + 2, 1:3) = [dx, -dy, z0]
               culm%nodcoord(itr + 3, 1:3) = [dx, dy, z0]
               culm%nodcoord(itr + 4, 1:3) = [-dx, dy, z0]
               itr = itr + 4
            end do

            itr = output_d%nn()
            culm%elemnod(1, 1:8) = &
               [itr + 1, itr + 2, itr + 3, itr + 4, 1, 2, 3, 4]

            do i = 2, culm_division

               culm%elemnod(i, 1:8) = &
                  [itr + 1, itr + 2, itr + 3, itr + 4, itr + 1 - 4, itr + 2 - 4, itr + 3 - 4, itr + 4 - 4]
               itr = itr + 4
            end do

            wheat%nodcoord = wheat%nodcoord.v.culm%nodcoord
            wheat%elemnod = wheat%elemnod.v.culm%elemnod

            output_d%mesh = wheat
            call output_d%move(z=-output_d%zmin())
         elseif (present(culm_diameter)) then
            ! create culm
            wheat = output_d%mesh
            culm%nodcoord = zeros(culm_division*4, 3)
            culm%elemnod = int(zeros(culm_division, 8))
            itr = 0

            do i = 1, culm_division
               theta = dble(i)/dble(culm_division)
               dx = culm_diameter*0.50d0*theta &
                    + panicle_panicle_diameter*0.50d0*(1.0d0 - theta)
               dy = culm_diameter*0.50d0*theta &
                    + panicle_panicle_diameter*0.50d0*(1.0d0 - theta)
               z0 = -theta*culm_length
               culm%nodcoord(itr + 1, 1:3) = [-dx, -dy, z0]
               culm%nodcoord(itr + 2, 1:3) = [dx, -dy, z0]
               culm%nodcoord(itr + 3, 1:3) = [dx, dy, z0]
               culm%nodcoord(itr + 4, 1:3) = [-dx, dy, z0]
               itr = itr + 4
            end do

            itr = output_d%nn()
            culm%elemnod(1, 1:8) = &
               [itr + 1, itr + 2, itr + 3, itr + 4, 1, 2, 3, 4]

            do i = 2, culm_division

               culm%elemnod(i, 1:8) = &
                  [itr + 1, itr + 2, itr + 3, itr + 4, itr + 1 - 4, itr + 2 - 4, itr + 3 - 4, itr + 4 - 4]
               itr = itr + 4
            end do

            wheat%nodcoord = wheat%nodcoord.v.culm%nodcoord
            wheat%elemnod = wheat%elemnod.v.culm%elemnod

            output_d%mesh = wheat
            call output_d%move(z=-output_d%zmin())
         end if
      end if

   end function

   function to_culm_diameter(heights_vs_diameters, height) result(ret)
      real(real64), intent(in) :: heights_vs_diameters(:, :) ! 1:height of observation points , 2:diameters
      real(real64), intent(in) :: height
      real(real64) :: ret
      real(real64), allocatable :: heights(:), diameters(:)
      real(real64) :: theta
      integer(int32) :: i

      heights = heights_vs_diameters(:, 1)
      diameters = heights_vs_diameters(:, 2)

      call heapsort(n=size(heights), array=heights, val=diameters)

      if (height < heights(1)) then
         ret = diameters(1)
         return
      end if

      if (height > heights(size(heights))) then
         ret = diameters(size(diameters))
         return
      end if

      do i = 1, size(heights) - 1

         if (heights(i) <= height .and. height <= heights(i + 1)) then
            theta = (height - heights(i))/(heights(i + 1) - heights(i))
            ret = theta*diameters(i + 1) + (1.0d0 - theta)*diameters(i)
            return
         else
            cycle
         end if
      end do

   end function

! ############################################################
   function append_panicle_object_vector(arg1, arg2) result(ret)
      type(Panicle_), allocatable, intent(in) :: arg1(:), arg2(:)
      type(Panicle_), allocatable :: ret(:)

      if (.not. allocated(arg1)) then
         if (.not. allocated(arg2)) then
            return
         else
            ret = arg2
         end if
      else
         if (.not. allocated(arg2)) then
            ret = arg1
            return
         else
            allocate (ret(size(arg1) + size(arg2)))
            ret(1:size(arg1)) = arg1(:)
            ret(size(arg1) + 1:) = arg2(:)
         end if
      end if

   end function
! ############################################################

end module
