module growthmod
   implicit none
contains
   !===================================
   subroutine initialize_gw3(d_connect, d_node, peti_coord, leaf_coord, max_leaf_angle, &
                             max_leaf_length, max_leaf_width, leaf_shape_rate, peti_length, max_length)
      integer, allocatable, intent(out)::d_connect(:, :)
      real(8), allocatable, intent(out)::d_node(:, :), peti_coord(:, :, :), leaf_coord(:, :, :, :)
      real(8), intent(in)::max_leaf_angle, &
                            max_leaf_length, max_leaf_width, leaf_shape_rate, peti_length, max_length
      real(8) r1, r2, x_coord, &
         y_coord, z_coord, direction, pi, angle, x_coord_1, y_coord_1, z_coord_1, &
         peti_RUE, mid_length, width, length
      real(8), allocatable::a(:), b(:), z(:), nvec(:), a_unit(:), b_unit(:)

      integer int_num, real_value

      integer i, j, n, m, o, parent, itr
      integer :: seedsize
      integer, allocatable :: seed(:)
      real :: rnd
      integer :: c !時間を入れる
      pi = 3.14159265350d0

      int_num = 3
      real_value = 5
      allocate (d_connect(2, int_num), d_node(2, real_value))
      allocate (peti_coord(2, 4, 3), leaf_coord(2, 3, 4, 3))
      allocate (a(3), b(3), z(3), nvec(3), a_unit(3), b_unit(3))
      leaf_coord(:, :, :, :) = 0.0d0
      peti_coord(:, :, :) = 0.0d0
      !initisalize d_connect
      d_connect(1, 1) = 1
      d_connect(1, 2) = -1 !Apical No.
      d_connect(1, 3) = 0

      d_connect(2, 1) = 1
      d_connect(2, 2) = 0 !Apical No.
      d_connect(2, 3) = 0

      !initialize d_node
      d_node(1, 1) = 0.0d0
      d_node(1, 2) = 0.0d0
      d_node(1, 3) = 0.0d0
      d_node(1, 4) = 1.0d0
      d_node(1, 5) = 0.0d0

      d_node(2, 1) = 0.0d0
      d_node(2, 2) = 0.0d0
      d_node(2, 3) = max_length
      d_node(2, 4) = 1.0d0
      d_node(2, 5) = 0.0d0

      !initialize peti_coord(:,:,:)============

      !random number generation
      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      length = ((rnd - 0.50d0) + 1.0d0)*peti_length

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      direction = rnd*360.0d0

      r1 = cos(pi*(angle/180.0d0))*length
      r2 = sin(pi*(angle/180.0d0))*length

      x_coord = d_node(2, 1)
      y_coord = d_node(2, 2)
      z_coord = d_node(2, 3)

      x_coord = x_coord + r2*cos(pi*(direction/180.0d0))
      y_coord = y_coord + r2*sin(pi*(direction/180.0d0))
      z_coord = z_coord + r1

      peti_coord(1, :, :) = 0.0d0

      peti_coord(2, 1, 1) = x_coord
      peti_coord(2, 1, 2) = y_coord
      peti_coord(2, 1, 3) = z_coord

      !second-fourth
      do i = 2, 4
         !random number generation
         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 2) then
            peti_RUE = 0.1
         else
            peti_RUE = 0.05
         end if

         length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length*peti_RUE

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 2) then
            direction = rnd*360.0d0
         elseif (i == 3) then
            direction = direction + 120
         elseif (i == 4) then
            direction = direction - 240
         end if
         r1 = cos(pi*(angle/180.0d0))*length
         r2 = sin(pi*(angle/180.0d0))*length

         x_coord_1 = x_coord
         y_coord_1 = y_coord
         z_coord_1 = z_coord

         x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0))
         y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0))
         z_coord_1 = z_coord_1 + r1

         peti_coord(2, i, 1) = x_coord_1
         peti_coord(2, i, 2) = y_coord_1
         peti_coord(2, i, 3) = z_coord_1
      end do

      !leaf generate
      do i = 1, 3
         !random number generation
         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length
         width = ((rnd - 0.50d0) + 1.0d0)*max_leaf_width

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 1) then
            direction = rnd*360.0d0
         elseif (i == 2) then
            direction = direction + 120
         elseif (i == 3) then
            direction = direction - 240
         end if

         r1 = cos(pi*(angle/180.0d0))*length
         r2 = sin(pi*(angle/180.0d0))*length

         x_coord_1 = peti_coord(2, i, 1)
         y_coord_1 = peti_coord(2, i, 2)
         z_coord_1 = peti_coord(2, i, 3)

         x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0))
         y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0))
         z_coord_1 = z_coord_1 + r1

         leaf_coord(2, i, 3, 1) = x_coord_1
         leaf_coord(2, i, 3, 2) = y_coord_1
         leaf_coord(2, i, 3, 3) = z_coord_1

         leaf_coord(2, i, 1, 1) = peti_coord(2, i, 1)
         leaf_coord(2, i, 1, 2) = peti_coord(2, i, 2)
         leaf_coord(2, i, 1, 3) = peti_coord(2, i, 3)
         !
         a(1) = leaf_coord(2, i, 3, 1) - peti_coord(2, i, 1)
         a(2) = leaf_coord(2, i, 3, 2) - peti_coord(2, i, 2)
         a(3) = leaf_coord(2, i, 3, 3) - peti_coord(2, i, 3)

         z(1) = 0.0d0
         z(2) = 0.0d0
         z(3) = 1.0d0

         a_unit(:) = a(:)/dsqrt(dot_product(a, a))

         b(:) = z(:) - dot_product(a, z)*a_unit(:)

         b_unit(:) = b(:)/dsqrt(dot_product(b, b))

         nvec(1) = a_unit(2)*b_unit(3) - a_unit(3)*b_unit(2)
         nvec(2) = a_unit(3)*b_unit(1) - a_unit(1)*b_unit(3)
         nvec(3) = a_unit(1)*b_unit(2) - a_unit(2)*b_unit(1)

         nvec(:) = nvec(:)/dsqrt(dot_product(nvec, nvec))

         mid_length = 1.0d0/(1.0d0 + leaf_shape_rate)*length

         j = 2

         x_coord_1 = peti_coord(2, i, 1)
         y_coord_1 = peti_coord(2, i, 2)
         z_coord_1 = peti_coord(2, i, 3)

         x_coord_1 = x_coord_1 + mid_length*a_unit(1)
         y_coord_1 = y_coord_1 + mid_length*a_unit(2)
         z_coord_1 = z_coord_1 + mid_length*a_unit(3)

         leaf_coord(2, i, j, 1) = x_coord_1 - 0.50d0*width*nvec(1)
         leaf_coord(2, i, j, 2) = y_coord_1 - 0.50d0*width*nvec(2)
         leaf_coord(2, i, j, 3) = z_coord_1 - 0.50d0*width*nvec(3)

         j = 4
         x_coord_1 = peti_coord(2, i, 1)
         y_coord_1 = peti_coord(2, i, 2)
         z_coord_1 = peti_coord(2, i, 3)

         x_coord_1 = x_coord_1 + mid_length*a_unit(1)
         y_coord_1 = y_coord_1 + mid_length*a_unit(2)
         z_coord_1 = z_coord_1 + mid_length*a_unit(3)

         leaf_coord(2, i, j, 1) = x_coord_1 + 0.50d0*width*nvec(1)
         leaf_coord(2, i, j, 2) = y_coord_1 + 0.50d0*width*nvec(2)
         leaf_coord(2, i, j, 3) = z_coord_1 + 0.50d0*width*nvec(3)

      end do

   end subroutine initialize_gw3
   !===================================
   subroutine growth_gw3(node, d_connect, d_node, max_angle, max_length, &
                         Apical_dominance_p, max_leaf_angle, &
                         max_leaf_length, max_leaf_width, leaf_shape_rate, &
                         peti_coord, leaf_coord, peti_length)

      integer, allocatable::d_connect_es(:, :)
      integer, allocatable, intent(inout)::d_connect(:, :)

      real(8), allocatable::d_node_es(:, :), peti_coord_es(:, :, :), &
                             leaf_coord_es(:, :, :, :)
      real(8), allocatable, intent(inout)::d_node(:, :), peti_coord(:, :, :), &
                                            leaf_coord(:, :, :, :)

      real(8), intent(in)::max_angle, max_length, peti_length
      real(8), intent(in)::max_leaf_angle, &
                            max_leaf_length, max_leaf_width, leaf_shape_rate
      real(8) angle, length, direction, r1, r2, pi, x_coord, y_coord, z_coord
      integer, intent(in)::node, Apical_dominance_p

      real(8) x_coord_1, y_coord_1, z_coord_1, peti_RUE, mid_length, width
      real(8), allocatable::a(:), b(:), z(:), nvec(:), a_unit(:), b_unit(:)

      integer i, j, n, m, o, parent, itr
      integer :: seedsize
      integer, allocatable :: seed(:)
      real :: rnd
      integer :: c !時間を入れる

      pi = 3.14159265350d0

      allocate (a(3), b(3), z(3), nvec(3), a_unit(3), b_unit(3))

      n = size(d_connect, 1)
      m = size(d_connect, 2)
      o = size(d_node, 2)
      allocate (d_connect_es(n + 1, m), d_node_es(n + 1, o))
      d_connect_es(:, :) = 0
      allocate (peti_coord_es(n, 4, 3), leaf_coord_es(n, 3, 4, 3))
      !copy
      do i = 1, n
         d_connect_es(i, :) = d_connect(i, :)
         d_node_es(i, :) = d_node(i, :)
      end do

      !generate new node
      d_connect_es(n + 1, 1) = node
      if (d_connect_es(node, 2) == -1) then
         !basis
         return
      elseif (d_connect_es(node, 2) == 0) then
         !aptical node => Get aptical dominance
         d_connect_es(n + 1, 2) = 0

         d_connect_es(node, 2) = 1
         parent = d_connect_es(node, 1)
         itr = 1
         do
            itr = itr + 1
            if (d_connect_es(parent, 2) == -1) then
               !base
               exit
            else
               if (d_connect_es(parent, 2) + 1 >= itr) then
                  d_connect_es(parent, 2) = d_connect_es(parent, 2) + 1
                  parent = d_connect_es(parent, 1)
               else
                  cycle
               end if
            end if
         end do

      elseif (d_connect_es(node, 2) >= Apical_dominance_p) then
         !generate new aptical
         d_connect_es(n + 1, 2) = 0
         parent = d_connect_es(node, 1)
         itr = 1
         do
            itr = itr + 1
            if (d_connect_es(parent, 2) == -1) then
               !base
               exit
            else
               if (d_connect_es(parent, 2) + 1 == itr) then
                  d_connect_es(parent, 2) = d_connect_es(parent, 2) + 1
                  parent = d_connect_es(parent, 1)
               else
                  exit
               end if
            end if
         end do
      else
         !not generate
         return
      end if

      d_connect_es(n + 1, 3) = 0

      !random number generation
      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      angle = ((rnd - 0.50d0) + 1.0d0)*max_angle

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      length = ((rnd - 0.50d0) + 1.0d0)*max_length

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      direction = rnd*360.0d0

      r1 = abs(cos(pi*(angle/180.0d0)))*length
      r2 = abs(sin(pi*(angle/180.0d0)))*length

      x_coord = d_node(node, 1)
      y_coord = d_node(node, 2)
      z_coord = d_node(node, 3)

      x_coord = x_coord + r2*cos(pi*(direction/180.0d0))
      y_coord = y_coord + r2*sin(pi*(direction/180.0d0))
      z_coord = z_coord + r1

      d_node_es(n + 1, 1) = x_coord
      d_node_es(n + 1, 2) = y_coord
      d_node_es(n + 1, 3) = z_coord
      d_node_es(n + 1, 4) = 1.0d0
      d_node_es(n + 1, 5) = 0.0d0

      deallocate (d_connect, d_node)
      n = size(d_connect_es, 1)
      m = size(d_connect_es, 2)
      o = size(d_node_es, 2)
      allocate (d_connect(n, m), d_node(n, o))
      d_connect(:, :) = d_connect_es(:, :)
      d_node(:, :) = d_node_es(:, :)
      deallocate (d_connect_es, d_node_es)

      !====================================
      !generate peti and leaf
      !========================

      !expand peti_coord,leaf_coord
      !copy
      n = size(d_connect, 1)
      leaf_coord_es(:, :, :, :) = leaf_coord(:, :, :, :)
      peti_coord_es(:, :, :) = peti_coord(:, :, :)
      deallocate (peti_coord, leaf_coord)
      allocate (peti_coord(n, 4, 3), leaf_coord(n, 3, 4, 3))

      do i = 1, size(peti_coord, 1)
         peti_coord(i, :, :) = peti_coord_es(i, :, :)
         leaf_coord(i, :, :, :) = leaf_coord_es(i, :, :, :)
      end do
      deallocate (peti_coord_es, leaf_coord_es)

      n = size(d_connect, 1)
      !random number generation
      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      length = ((rnd - 0.50d0) + 1.0d0)*peti_length

      call system_clock(count=c) !時間を取得
      call random_seed(size=seedsize)
      allocate (seed(seedsize))
      call random_seed(get=seed)
      seed = c !時間を全部に代入
      call random_number(rnd) !rndに乱数をセット
      deallocate (seed)

      direction = rnd*360.0d0

      r1 = cos(pi*(angle/180.0d0))*length
      r2 = sin(pi*(angle/180.0d0))*length

      x_coord = d_node(n, 1)
      y_coord = d_node(n, 2)
      z_coord = d_node(n, 3)

      x_coord = x_coord + r2*cos(pi*(direction/180.0d0))
      y_coord = y_coord + r2*sin(pi*(direction/180.0d0))
      z_coord = z_coord + r1

      peti_coord(1, :, :) = 0.0d0

      peti_coord(n, 1, 1) = x_coord
      peti_coord(n, 1, 2) = y_coord
      peti_coord(n, 1, 3) = z_coord

      !second-fourth
      do i = 2, 4
         !random number generation
         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 2) then
            peti_RUE = 0.10d0
         else
            peti_RUE = 0.050d0
         end if

         length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length*peti_RUE

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 2) then
            direction = rnd*360.0d0
         elseif (i == 3) then
            direction = direction + 120
         elseif (i == 4) then
            direction = direction - 240
         end if

         r1 = abs(cos(pi*(angle/180.0d0)))*length
         r2 = abs(sin(pi*(angle/180.0d0)))*length

         x_coord_1 = x_coord
         y_coord_1 = y_coord
         z_coord_1 = z_coord

         x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0))
         y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0))
         z_coord_1 = z_coord_1 + r1

         peti_coord(n, i, 1) = x_coord_1
         peti_coord(n, i, 2) = y_coord_1
         peti_coord(n, i, 3) = z_coord_1
      end do

      !leaf generate
      do i = 1, 3
         !random number generation
         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length
         width = ((rnd - 0.50d0) + 1.0d0)*max_leaf_width

         call system_clock(count=c) !時間を取得
         call random_seed(size=seedsize)
         allocate (seed(seedsize))
         call random_seed(get=seed)
         seed = c !時間を全部に代入
         call random_number(rnd) !rndに乱数をセット
         deallocate (seed)

         if (i == 1) then
            direction = rnd*360.0d0
         elseif (i == 2) then
            direction = direction + 120
         elseif (i == 3) then
            direction = direction - 240
         end if

         r1 = abs(cos(pi*(angle/180.0d0)))*length
         r2 = abs(sin(pi*(angle/180.0d0)))*length

         x_coord_1 = peti_coord(n, i, 1)
         y_coord_1 = peti_coord(n, i, 2)
         z_coord_1 = peti_coord(n, i, 3)

         x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0))
         y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0))
         z_coord_1 = z_coord_1 + r1

         leaf_coord(n, i, 3, 1) = x_coord_1
         leaf_coord(n, i, 3, 2) = y_coord_1
         leaf_coord(n, i, 3, 3) = z_coord_1

         leaf_coord(n, i, 1, 1) = peti_coord(n, i, 1)
         leaf_coord(n, i, 1, 2) = peti_coord(n, i, 2)
         leaf_coord(n, i, 1, 3) = peti_coord(n, i, 3)
         !
         a(1) = leaf_coord(n, i, 3, 1) - peti_coord(n, i, 1)
         a(2) = leaf_coord(n, i, 3, 2) - peti_coord(n, i, 2)
         a(3) = leaf_coord(n, i, 3, 3) - peti_coord(n, i, 3)

         z(1) = 0.0d0
         z(2) = 0.0d0
         z(3) = 1.0d0

         a_unit(:) = a(:)/dsqrt(dot_product(a, a))

         b(:) = z(:) - dot_product(a, z)*a_unit(:)

         b_unit(:) = b(:)/dsqrt(dot_product(b, b))

         nvec(1) = a_unit(2)*b_unit(3) - a_unit(3)*b_unit(2)
         nvec(2) = a_unit(3)*b_unit(1) - a_unit(1)*b_unit(3)
         nvec(3) = a_unit(1)*b_unit(2) - a_unit(2)*b_unit(1)
         nvec(:) = nvec(:)/dsqrt(dot_product(nvec, nvec))
         mid_length = 1.0d0/(1.0d0 + leaf_shape_rate)*length

         j = 2

         x_coord_1 = peti_coord(n, i, 1)
         y_coord_1 = peti_coord(n, i, 2)
         z_coord_1 = peti_coord(n, i, 3)

         x_coord_1 = x_coord_1 + mid_length*a_unit(1)
         y_coord_1 = y_coord_1 + mid_length*a_unit(2)
         z_coord_1 = z_coord_1 + mid_length*a_unit(3)

         leaf_coord(n, i, j, 1) = x_coord_1 - 0.50d0*width*nvec(1)
         leaf_coord(n, i, j, 2) = y_coord_1 - 0.50d0*width*nvec(2)
         leaf_coord(n, i, j, 3) = z_coord_1 - 0.50d0*width*nvec(3)

         j = 4
         x_coord_1 = peti_coord(n, i, 1)
         y_coord_1 = peti_coord(n, i, 2)
         z_coord_1 = peti_coord(n, i, 3)

         x_coord_1 = x_coord_1 + mid_length*a_unit(1)
         y_coord_1 = y_coord_1 + mid_length*a_unit(2)
         z_coord_1 = z_coord_1 + mid_length*a_unit(3)

         leaf_coord(n, i, j, 1) = x_coord_1 + 0.50d0*width*nvec(1)
         leaf_coord(n, i, j, 2) = y_coord_1 + 0.50d0*width*nvec(2)
         leaf_coord(n, i, j, 3) = z_coord_1 + 0.50d0*width*nvec(3)

      end do

   end subroutine growth_gw3
   !===================================
   subroutine update_factor_gw3(d_connect, d_node, RUE, num_node, leaf_coord &
                                , day, daylength, limit_leaf_day, day_ID)
      real(8), intent(inout)::d_node(:, :)
      real(8), intent(in)::RUE, leaf_coord(:, :, :, :), daylength(:)
      real(8), allocatable::photosynthesis(:, :), radiation(:, :)
      real(8) factor, light_area, total_area, old_RUEmeter
      integer, intent(in)::d_connect(:, :), num_node, day_ID, day(:, :), limit_leaf_day
      integer i, j, month, node_ID, leaf_ID

      allocate (photosynthesis(num_node, 3), radiation(12, 2))
      photosynthesis(:, :) = 0.0d0
      open (60, file="radiation.txt")
      do i = 1, 12
         read (60, *) radiation(i, 1:2)
      end do
      close (60)

      !生成してから1日経ってから光合成開始
      !日照面積判定
      do i = 1, num_node
         do j = 1, size(leaf_coord, 2)
            !compute light area
            node_ID = i
            leaf_ID = j
            call detect_light_area(leaf_coord, light_area, total_area, node_ID, leaf_ID)

            photosynthesis(i, 1) = light_area
            photosynthesis(i, 2) = total_area
            !光合成量＝日射時間(hr)×面積当たり時間当たり日射量(MJ/m^2/hr)×日照面積(m^2)×日射量あたり光合成(g/MJ)
            month = day(day_ID, 2)
            if (radiation(month, 1) == 0.0d0) then
               photosynthesis(i, 3) = 0.0d0
            end if

            if (d_connect(i, 3) == -1) then
               old_RUEmeter = 0.0d0
            else
               old_RUEmeter = dble(d_connect(i, 3))/dble(limit_leaf_day)
               if (old_RUEmeter < 0) stop"debugb"
            end if

            photosynthesis(i, 3) = daylength(day_ID)*radiation(month, 2)/radiation(month, 1) &
                                   *light_area*RUE*old_RUEmeter
         end do
      end do

      do i = 1, num_node
         !compute factor
         factor = d_node(i, 5)

         factor = factor + photosynthesis(i, 3)

         d_node(i, 5) = factor
      end do

   end subroutine update_factor_gw3
   !===================================
   subroutine detect_start(start, finish, seedling_day, harvesting_day, day)
      integer, intent(out)::start, finish
      integer, intent(in)::seedling_day(:), harvesting_day(:), day(:, :)
      integer i, j, n, yy, mm, dd, diff

      do i = 1, size(day, 1)
         yy = day(i, 1)
         mm = day(i, 2)
         dd = day(i, 3)

         diff = abs(seedling_day(1) - yy) &
                + abs(seedling_day(2) - mm) &
                + abs(seedling_day(3) - dd)

         if (diff == 0) then
            start = i
            exit
         end if

      end do

      do i = 1, size(day, 1)
         yy = day(i, 1)
         mm = day(i, 2)
         dd = day(i, 3)

         diff = abs(harvesting_day(1) - yy) &
                + abs(harvesting_day(2) - mm) &
                + abs(harvesting_day(3) - dd)

         if (diff == 0) then
            finish = i
            exit
         end if

      end do

   end subroutine detect_start
   !===================================
   subroutine count_br_per_node(node_number, br_per_node, d_connect)
      integer, intent(in)::node_number, d_connect(:, :)
      integer, intent(out)::br_per_node
      integer i

      br_per_node = 0
      do i = 1, size(d_connect, 1)
         if (node_number == d_connect(i, 1)) then
            br_per_node = br_per_node + 1
         else
            cycle
         end if
      end do
      if (br_per_node >= 3) stop"debug"

   end subroutine count_br_per_node
   !===================================
   subroutine detect_light_area(leaf_coord, light_area, total_area, node_ID, leaf_ID)
      real(8), intent(in)::leaf_coord(:, :, :, :)
      real(8), intent(out)::light_area, total_area
      real(8), allocatable::mid_x(:), shape_function(:), evaluate_point(:, :), eva_x(:), rvec(:) &
                             , x1(:), x2(:), y1(:), y2(:)
      real(8) gzi_1, gzi_2, eva_r, mid_r, both_r, l1, l2, l3
      integer, intent(in)::node_ID, leaf_ID
      integer i, j, k, l, m, shade_or_not, count_shade

      allocate (mid_x(3), shape_function(4), evaluate_point(4, 2), eva_x(3), rvec(2), x1(2), x2(2))
      allocate (y1(3), y2(3))

      !2D linear interpolation,4 points
      evaluate_point(1, 1) = -0.50d0
      evaluate_point(1, 2) = -0.50d0
      evaluate_point(2, 1) = -0.50d0
      evaluate_point(2, 2) = 0.50d0
      evaluate_point(3, 1) = 0.50d0
      evaluate_point(3, 2) = 0.50d0
      evaluate_point(4, 1) = 0.50d0
      evaluate_point(4, 2) = -0.50d0

      y1(1:3) = leaf_coord(node_ID, leaf_ID, 1, 1:3) - leaf_coord(node_ID, leaf_ID, 3, 1:3)
      y2(1:3) = leaf_coord(node_ID, leaf_ID, 2, 1:3) - leaf_coord(node_ID, leaf_ID, 4, 1:3)
      l1 = dsqrt(dot_product(y1, y1))
      l2 = dsqrt(dot_product(y2, y2))
      total_area = 0.50d0*l1*l2/100.0d0/100.0d0

      count_shade = 0
      do k = 1, 4 !evaluate points
         gzi_1 = evaluate_point(k, 1)
         gzi_2 = evaluate_point(k, 2)
         shape_function(1) = 0.250d0*(1.0d0 - gzi_1)*(1.0d0 - gzi_1)
         shape_function(2) = 0.250d0*(1.0d0 - gzi_1)*(1.0d0 + gzi_1)
         shape_function(3) = 0.250d0*(1.0d0 + gzi_1)*(1.0d0 + gzi_1)
         shape_function(4) = 0.250d0*(1.0d0 + gzi_1)*(1.0d0 - gzi_1)

         eva_x(:) = 0.0d0
         do l = 1, 4
            eva_x(:) = eva_x(:) + leaf_coord(node_ID, leaf_ID, l, :) &
                       *shape_function(l)
         end do
         rvec(:) = leaf_coord(node_ID, leaf_ID, k, 1:2)
         eva_r = dsqrt(dot_product(rvec - eva_x(1:2), rvec - eva_x(1:2)))

         shade_or_not = 0

         do i = 1, size(leaf_coord, 1)
            do j = 1, size(leaf_coord, 2)
               mid_x(:) = 0.0d0
               do l = 1, 4
                  mid_x(:) = mid_x(:) + 0.250d0*leaf_coord(i, j, l, :)
               end do

               !if z coordinate of the midpoint of the leaf > evaluate point
               ! => possible node
               if (mid_x(3) < eva_x(3)) then
                  cycle
               end if

               x1(1:2) = leaf_coord(i, j, 1, 1:2) - leaf_coord(i, j, 3, 1:2)
               x2(1:2) = leaf_coord(i, j, 2, 1:2) - leaf_coord(i, j, 4, 1:2)

               mid_r = dsqrt(dot_product(x1, x1)) + dsqrt(dot_product(x2, x2))
               mid_r = 0.50d0*mid_r

               both_r = dsqrt(dot_product(mid_x(1:2) - eva_x(1:2), mid_x(1:2) - eva_x(1:2)))

               if (both_r > eva_r + mid_r) then
                  cycle
               else
                  shade_or_not = 1
                  exit
               end if
            end do

            if (shade_or_not == 1) then
               exit
            end if
         end do

         if (shade_or_not == 1) then
            cycle
         else
            count_shade = count_shade + 1
         end if
      end do

      light_area = total_area/4.0d0*dble(count_shade)

   end subroutine detect_light_area
   !===================================
   subroutine get_old(d_connect, limit_leaf_day)
      integer, intent(inout)::d_connect(:, :)
      integer, intent(in)::limit_leaf_day
      integer i

      do i = 1, size(d_connect, 1)
         if (d_connect(i, 3) == 0) then
            !first day
            d_connect(i, 3) = limit_leaf_day
         elseif (d_connect(i, 3) == 1) then
            d_connect(i, 3) = -1
         elseif (d_connect(i, 3) == -1) then
            d_connect(i, 3) = -1
         elseif (limit_leaf_day >= d_connect(i, 3) .and. d_connect(i, 3) >= 2) then
            d_connect(i, 3) = d_connect(i, 3) - 1
         else
            stop"invalid d_connect,3"
         end if
      end do

   end subroutine get_old
   !===================================
   subroutine generate_pod(d_node, d_pod, source_convection_rate, source_limitation)
      real(8), intent(inout)::d_node(:, :), d_pod(:, :)
      real(8), intent(in)::source_convection_rate, source_limitation
      real(8) total_source, extra_source
      integer i, j, pn

      do i = 1, size(d_node, 1)
         total_source = d_node(i, 5)
         pn = 0
         do j = 1, size(d_pod, 2)
            if (d_pod(i, j) == source_limitation) then
               cycle
            elseif (d_pod(i, j) < source_limitation) then
               pn = j
               exit
            else
               stop"invalid source/pod"
            end if
         end do

         if (pn == 0) then
            !sink was occupied
            cycle
         end if

         if (total_source >= source_convection_rate) then
            d_node(i, 5) = total_source - source_convection_rate
            d_pod(i, pn) = d_pod(i, pn) + source_convection_rate

         elseif (source_convection_rate > total_source .and. &
                 total_source > 0.0d0) then
            d_pod(i, pn) = d_pod(i, pn) + total_source
            d_node(i, 5) = 0.0d0
         elseif (total_source == 0.0d0) then
            cycle
         else
            stop"invalid source < 0"
         end if

         extra_source = 0.0d0
         do j = 1, size(d_pod, 2)
            if (extra_source /= 0.0d0) then
               d_pod(i, j) = d_pod(i, j) + extra_source
            end if

            if (d_pod(i, j) > source_limitation) then
               extra_source = d_pod(i, j) - source_limitation
               d_pod(i, j) = source_limitation
            end if
         end do

         if (extra_source /= 0.0d0) then
            d_node(i, 5) = d_node(i, 5) + extra_source
            extra_source = 0.0d0
         end if

      end do

   end subroutine generate_pod
   !===================================
end module growthmod
