module RiceClass
   use, intrinsic :: iso_fortran_env
   use MathClass
   use SeedClass
   use LeafClass
   use RootClass
   use SoilClass
   use LightClass
   use PlantNodeClass
   implicit none

   type :: Rice_
      ! growth_habit = determinate, indeterminate, semi-indeterminate, or vine
      character*20 :: growth_habit
      character*2  :: growth_stage
      integer(int32) :: Num_Of_Node
      integer(int32) :: Num_Of_Root

      integer(int32) :: MaxLeafNum = 300
      integer(int32) :: MaxRootNum = 300
      integer(int32) :: MaxStemNum = 300

      integer(int32)  :: ms_node, br_node(300), br_from(300)
      real(real64)    :: ms_length, br_length(300)
      real(real64)    :: ms_width, br_width(300)
      real(real64)    :: ms_angle_ave, br_angle_ave(300)
      real(real64)    :: ms_angle_sig, br_angle_sig(300)

      integer(int32)  :: mr_node, brr_node(300), brr_from(300)
      real(real64)    :: mr_length, brr_length(300)
      real(real64)    :: mr_width, brr_width(300)
      real(real64)    :: mr_angle_ave, brr_angle_ave(300)
      real(real64)    :: mr_angle_sig, brr_angle_sig(300)

      real(real64)    :: peti_size_ave(300)
      real(real64)    :: peti_size_sig(300)
      real(real64)    :: peti_width_ave(300)
      real(real64)    :: peti_width_sig(300)
      real(real64)    :: peti_angle_ave(300)
      real(real64)    :: peti_angle_sig(300)

      real(real64)    :: leaf_angle_ave(300*3)
      real(real64)    :: leaf_angle_sig(300*3)
      real(real64)    :: leaf_length_ave(300*3)
      real(real64)    :: leaf_length_sig(300*3)
      real(real64)    :: leaf_width_ave(300*3)
      real(real64)    :: leaf_width_sig(300*3)
      real(real64)    :: leaf_thickness_ave(300*3)
      real(real64)    :: leaf_thickness_sig(300*3)

      character(3) :: Stage ! VE, CV, V1,V2, ..., R1, R2, ..., R8
      character(200) :: name
      integer(int32)::stage_id = 0
      real(real64) :: dt
      type(Seed_) :: Seed
      type(PlantNode_), allocatable :: NodeSystem(:)
      type(PlantRoot_), allocatable :: RootSystem(:)

      type(Stem_), allocatable :: Stem(:)
      type(Leaf_), allocatable :: Leaf(:)
      type(Root_), allocatable :: Root(:)
      type(Soil_), allocatable :: Soil

      ! 節-節点データ構造
      type(Mesh_) :: struct
      integer(int32), allocatable :: leaf2stem(:, :)
      integer(int32), allocatable :: stem2stem(:, :)
      integer(int32), allocatable :: root2stem(:, :)
      integer(int32), allocatable :: root2root(:, :)

      ! 器官オブジェクト配列
      type(FEMDomain_), allocatable :: leaf_list(:)
      type(FEMDomain_), allocatable :: stem_list(:)
      type(FEMDomain_), allocatable :: root_list(:)
      real(real64) :: time
      real(real64) :: seed_length
      real(real64) :: seed_width
      real(real64) :: seed_height
      real(real64), allocatable :: stem_angle(:, :)
      real(real64), allocatable :: root_angle(:, :)
      real(real64), allocatable :: leaf_angle(:, :)

      character(200) :: stemconfig = " "
      character(200) :: rootconfig = " "
      character(200) :: leafconfig = " "
   contains
      procedure, public :: addStem => addStemRice
      !procedure,public :: addRoot => addRootRice
      !procedure,public :: addLeaf => addLeafRice

      procedure, public :: Init => initRice
      procedure, public :: new => initRice
      procedure, public :: sowing => initRice
      procedure, public :: export => exportRice

      procedure, public :: grow => growRice

      procedure, public :: show => showRice
      procedure, public :: gmsh => gmshRice
      procedure, public :: msh => mshRice
      procedure, public :: stl => stlRice
      procedure, public :: json => jsonRice

      procedure, public :: WaterAbsorption => WaterAbsorptionRice
      procedure, public :: move => moveRice

      procedure, public :: numleaf => numleafRice
      procedure, public :: numstem => numstemRice
      procedure, public :: numroot => numrootRice

      procedure, public :: laytracing => laytracingRice
      procedure, public :: SinkSourceFlow => SinkSourceFlowRice

      !procedure,public :: AddNode => AddNodeRice
   end type

   type :: RiceCanopy_
      real(real64) :: inter_row, intra_row
      type(Rice_), allocatable :: Canopy(:, :)
   end type

contains

! ########################################
   subroutine initRice(obj, config, &
                       regacy, mass, water_content, radius, location, x, y, z, &
                       PlantRoot_diameter_per_seed_radius, max_PlantNode_num, Variety, FileName, &
                       max_leaf_num, max_stem_num, max_root_num)
      class(Rice_), intent(inout) :: obj

      real(real64), optional, intent(in) :: mass, water_content, radius, location(3), x, y, z
      real(real64), optional, intent(in) :: PlantRoot_diameter_per_seed_radius
      character(*), optional, intent(in) :: Variety, FileName, config
      logical, optional, intent(in) :: regacy
      character(200) :: fn, conf, line
      integer(int32), optional, intent(in) :: max_PlantNode_num, max_leaf_num, max_stem_num, max_root_num
      real(real64) :: MaxThickness, Maxwidth, loc(3), vec(3), rot(3), zaxis(3), meshloc(3), meshvec(3)
      integer(int32) :: i, j, k, blcount, id, rmc, n, node_id, node_id2, elemid, branch_id, num_stem_node
      integer(int32) :: num_leaf
      real(real64)::readvalreal
      integer(int32) :: readvalint
      type(IO_) :: soyconf
      type(Random_) :: random

      ! set default parameters
      ! stem
      obj%br_node(:) = 0
      obj%br_from(:) = 0
      obj%br_length(:) = 0.0d0

      obj%br_angle_ave(:) = 0.0d0
      obj%br_angle_sig(:) = 10.0d0
      obj%br_angle_ave(1) = 30.0d0
      obj%br_angle_sig(1) = 2.0d0

      obj%ms_angle_ave = 0.0d0
      obj%ms_angle_sig = 2.0d0

      ! for roots
      obj%brr_node(:) = 0
      obj%brr_from(:) = 0
      obj%brr_length(:) = 0.0d0

      obj%brr_angle_ave(:) = 0.0d0
      obj%brr_angle_sig(:) = 10.0d0
      obj%brr_angle_ave(1) = 30.0d0
      obj%brr_angle_sig(1) = 2.0d0

      obj%mr_angle_ave = 0.0d0
      obj%mr_angle_sig = 2.0d0
      ! peti
      ! is also stem

      obj%peti_size_ave(:) = 0.20d0
      obj%peti_size_sig(:) = 0.010d0

      obj%peti_width_ave(:) = 0.0050d0
      obj%peti_width_sig(:) = 0.00010d0

      obj%peti_angle_ave(:) = 30.0d0
      obj%peti_angle_sig(:) = 1.00d0

      ! leaf
      obj%leaf_length_ave(:) = 0.20d0
      obj%leaf_length_sig(:) = 0.01d0

      obj%leaf_width_ave(:) = 0.050d0
      obj%leaf_width_sig(:) = 0.010d0

      obj%leaf_thickness_ave(:) = 0.00100d0
      obj%leaf_thickness_sig(:) = 0.00050d0

      obj%leaf_angle_ave(:) = 40.0d0
      obj%leaf_angle_sig(:) = 10.0d0

      ! 子葉節、初生葉節、根の第1節まで種子の状態で存在

      ! 節を生成するためのスクリプトを開く
      if (.not. present(config) .or. index(config, ".json") == 0) then
         ! デフォルトの設定を生成
         print *, "New Rice-configuration >> soyconfig.json"
         call soyconf%open("soyconfig.json")
         write (soyconf%fh, *) '{'
         write (soyconf%fh, *) '   "type": "Rice",'
         write (soyconf%fh, *) '   "stemconfig": "stemconfig.json",'
         write (soyconf%fh, *) '   "rootconfig": "rootconfig.json",'
         write (soyconf%fh, *) '   "leafconfig": "leafconfig.json",'
         write (soyconf%fh, *) '   "stage": 0,'
         write (soyconf%fh, *) '   "length": 0.0090,'
         write (soyconf%fh, *) '   "width" : 0.0081,'
         write (soyconf%fh, *) '   "height": 0.0072,'
         write (soyconf%fh, *) '   "MaxLeafNum": 50,'
         write (soyconf%fh, *) '   "MaxRootNum":200,'
         write (soyconf%fh, *) '   "MaxStemNum": 50,'

         ! stem
         write (soyconf%fh, *) '   "br_node" : 0,'
         write (soyconf%fh, *) '   "br_from" : 0,'
         write (soyconf%fh, *) '   "br_length" : 0.00,'
         write (soyconf%fh, *) '   "br_angle_ave" : 0.00,'
         write (soyconf%fh, *) '   "br_angle_sig" : 10.00,'
         write (soyconf%fh, *) '   "br_angle_ave(1)": 360.00,'
         write (soyconf%fh, *) '   "br_angle_sig(1)": 2.00,'
         write (soyconf%fh, *) '   "ms_angle_ave": 0.00,'
         write (soyconf%fh, *) '   "ms_angle_sig": 2.00,'

         ! root
         write (soyconf%fh, *) '   "brr_node" : 0,'
         write (soyconf%fh, *) '   "brr_from" : 0,'
         write (soyconf%fh, *) '   "brr_length" : 0.00,'
         write (soyconf%fh, *) '   "brr_angle_ave" : 0.00,'
         write (soyconf%fh, *) '   "brr_angle_sig" : 10.00,'
         write (soyconf%fh, *) '   "brr_angle_ave(1)": 360.00,'
         write (soyconf%fh, *) '   "brr_angle_sig(1)": 2.00,'
         write (soyconf%fh, *) '   "mr_angle_ave": 0.00,'
         write (soyconf%fh, *) '   "mr_angle_sig": 2.00,'
         ! peti
         ! is also stem
         write (soyconf%fh, *) '   "peti_size_ave"  :  0.200,'
         write (soyconf%fh, *) '   "peti_size_sig"  :  0.0100,'
         write (soyconf%fh, *) '   "peti_width_ave"  :  0.00500,'
         write (soyconf%fh, *) '   "peti_width_sig"  :  0.000100,'
         write (soyconf%fh, *) '   "peti_angle_ave"  :  30.00,'
         write (soyconf%fh, *) '   "peti_angle_sig"  :  1.000,'
         ! leaf
         write (soyconf%fh, *) '   "leaf_length_ave"  :  0.200,'
         write (soyconf%fh, *) '   "leaf_length_sig"  :  0.010,'
         write (soyconf%fh, *) '   "leaf_width_ave"  :  0.0500,'
         write (soyconf%fh, *) '   "leaf_width_sig"  :  0.0100,'
         write (soyconf%fh, *) '   "leaf_thickness_ave"  :  0.001000,'
         write (soyconf%fh, *) '   "leaf_thickness_sig"  :  0.000500,'
         write (soyconf%fh, *) '   "leaf_angle_ave"  :  40.00,'
         write (soyconf%fh, *) '   "leaf_angle_sig"  :  10.00'
         write (soyconf%fh, *) '}'
         conf = "soyconfig.json"
         call soyconf%close()
      else
         conf = config
      end if

      call soyconf%open(conf)
      blcount = 0
      do
         read (soyconf%fh, '(a)') line
         print *, line
         if (adjustl(line) == "{") then
            blcount = 1
            cycle
         end if
         if (adjustl(line) == "}") then
            exit
         end if

         if (blcount == 1) then

            if (index(line, "Name") /= 0) then
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%name
            end if

            if (index(line, "Mainstem") /= 0) then
               do
                  read (soyconf%fh, '(a)') line
                  print *, line
                  if (index(line, "}") /= 0) then
                     exit
                  end if

                  if (index(line, "Length") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%ms_length
                  end if

                  if (index(line, "Width") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%ms_width
                  end if

                  if (index(line, "Node") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%ms_node
                  end if

               end do
            end if

            if (index(line, "Branch#") /= 0) then
               rmc = index(line, "{")
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, '"')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, '"')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, ':')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, "#")
               print *, line
               read (line(id + 1:), *) branch_id

               do
                  read (soyconf%fh, '(a)') line
                  print *, line
                  if (index(line, "}") /= 0) then
                     exit
                  end if

                  if (index(line, "Length") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%br_length(branch_id)
                  end if

                  if (index(line, "Width") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%br_Width(branch_id)
                  end if

                  if (index(line, "Node") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%br_node(branch_id)
                  end if

                  if (index(line, "From") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%br_from(branch_id)
                  end if

               end do
            end if

            ! for roots

            if (index(line, "Mainroot") /= 0) then
               do
                  read (soyconf%fh, '(a)') line
                  print *, line
                  if (index(line, "}") /= 0) then
                     exit
                  end if

                  if (index(line, "Length") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%mr_length
                  end if

                  if (index(line, "Width") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%mr_width
                  end if

                  if (index(line, "Node") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%mr_node
                  end if

               end do
            end if

            if (index(line, "Branchroot#") /= 0) then
               rmc = index(line, "{")
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, '"')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, '"')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               rmc = index(line, ':')
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, "#")
               print *, line
               read (line(id + 1:), *) branch_id

               do
                  read (soyconf%fh, '(a)') line
                  print *, line
                  if (index(line, "}") /= 0) then
                     exit
                  end if

                  if (index(line, "Length") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%brr_length(branch_id)
                  end if

                  if (index(line, "Width") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%brr_Width(branch_id)
                  end if

                  if (index(line, "Node") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%brr_node(branch_id)
                  end if

                  if (index(line, "From") /= 0) then
                     rmc = index(line, ",")
                     if (rmc /= 0) then
                        line(rmc:rmc) = " "
                     end if
                     id = index(line, ":")
                     read (line(id + 1:), *) obj%brr_from(branch_id)
                  end if

               end do
            end if

            if (index(line, "type") /= 0 .and. index(line, "Rice") == 0) then
               print *, "ERROR: This config-file is not for Rice"
               return
            end if

            if (index(line, "rootconfig") /= 0) then
               ! 茎の設定ファイル
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%rootconfig
            end if

            if (index(line, "stemconfig") /= 0) then
               ! 茎の設定ファイル
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%stemconfig
            end if

            if (index(line, "leafconfig") /= 0) then
               ! 茎の設定ファイル
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%leafconfig
            end if

            if (index(line, "stage") /= 0) then
               ! 生育ステージ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%stage_id
            end if

            if (index(line, "MaxLeafNum") /= 0) then
               ! 生育ステージ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%MaxLeafNum
            end if

            if (index(line, "MaxStemNum") /= 0) then
               ! 生育ステージ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%MaxStemNum
            end if

            if (index(line, "MaxRootNum") /= 0) then
               ! 生育ステージ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%MaxRootNum
            end if

            if (index(line, "length") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%seed_length
            end if

            if (index(line, "width") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%seed_width
            end if

            if (index(line, "height") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) obj%seed_height
            end if

            ! for version 2020.11.24

            ! stem
            if (index(line, "br_angle_ave") /= 0 .and. index(line, "br_angle_ave(") == 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%br_angle_ave(:) = readvalreal
            end if

            if (index(line, "br_angle_sig") /= 0 .and. index(line, "br_angle_sig(") == 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%br_angle_sig(:) = readvalreal
            end if

            if (index(line, "br_angle_ave(1)") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%br_angle_ave(1) = readvalreal
            end if
            if (index(line, "br_angle_sig(1)") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%br_angle_sig(1) = readvalreal
            end if

            if (index(line, "ms_angle_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%ms_angle_ave = readvalreal
            end if

            if (index(line, "ms_angle_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%ms_angle_sig = readvalreal
            end if
            ! peti
            ! is also stem

            if (index(line, "peti_size_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_size_ave(:) = readvalreal
            end if

            if (index(line, "peti_size_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_size_sig(:) = readvalreal
            end if

            if (index(line, "peti_width_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_width_ave(:) = readvalreal
            end if

            if (index(line, "peti_width_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_width_sig(:) = readvalreal
            end if

            if (index(line, "peti_angle_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_angle_ave(:) = readvalreal
            end if

            if (index(line, "peti_angle_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%peti_angle_sig(:) = readvalreal
            end if
            ! leaf

            if (index(line, "leaf_length_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_length_ave(:) = readvalreal
            end if

            if (index(line, "leaf_length_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_length_sig(:) = readvalreal
            end if

            if (index(line, "leaf_width_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_width_ave(:) = readvalreal
            end if

            if (index(line, "leaf_width_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_width_sig(:) = readvalreal
            end if

            if (index(line, "leaf_thickness_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_thickness_ave(:) = readvalreal
            end if

            if (index(line, "leaf_thickness_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_thickness_sig(:) = readvalreal
            end if

            if (index(line, "leaf_angle_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_angle_ave(:) = readvalreal
            end if

            if (index(line, "leaf_angle_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%leaf_angle_sig(:) = readvalreal
            end if

            ! added in 2020/12/15
            ! for roots

            if (index(line, "brr_angle_ave") /= 0 .and. index(line, "brr_angle_ave(") == 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%brr_angle_ave(:) = readvalreal
            end if

            if (index(line, "brr_angle_sig") /= 0 .and. index(line, "brr_angle_sig(") == 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%brr_angle_sig(:) = readvalreal
            end if

            if (index(line, "brr_angle_ave(1)") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%brr_angle_ave(1) = readvalreal
            end if
            if (index(line, "brr_angle_sig(1)") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%brr_angle_sig(1) = readvalreal
            end if

            if (index(line, "mr_angle_ave") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%mr_angle_ave = readvalreal
            end if

            if (index(line, "mr_angle_sig") /= 0) then
               ! 種子の長さ
               rmc = index(line, ",")
               ! カンマがあれば除く
               if (rmc /= 0) then
                  line(rmc:rmc) = " "
               end if
               id = index(line, ":")
               read (line(id + 1:), *) readvalreal
               obj%mr_angle_sig = readvalreal
            end if

            cycle

         end if

      end do
      call soyconf%close()

      if (index(config, ".json") == 0) then
         obj%stemconfig = " "
         obj%rootconfig = " "
         obj%leafconfig = " "
      end if

      if (obj%ms_node /= 0) then
         ! loaded from Mainstem-Branches relation file format
         ! ex.
!       {
!           "Name":"Rice",
!           "Mainstem":{
!               "Length":1.2,
!               "Node":13
!           },
!           "Branch#1":{
!               "From":1,
!               "Length":0.6,
!               "Node":7
!           },
!           "Branch#2":{
!               "From":3,
!               "Length":0.2,
!               "Node":2
!           },
!           "Branch#3":{
!               "From":4,
!               "Length":0.2,
!               "Node":2
!           }
!       }
         ! count number of nodes
         !num_node = countif(obj%ms_node,notEquai=.true.,0)
         !num_node = num_node + countif(obj%br_node,notEquai=.true.,0)

         allocate (obj%leaf(obj%MaxLeafNum))
         allocate (obj%root(obj%MaxrootNum))
         allocate (obj%stem(obj%MaxstemNum))
         allocate (obj%stem2stem(obj%MaxstemNum, obj%MaxstemNum))
         allocate (obj%leaf2stem(obj%MaxstemNum, obj%MaxLeafNum))
         allocate (obj%root2stem(obj%MaxrootNum, obj%MaxstemNum))
         allocate (obj%root2root(obj%MaxrootNum, obj%MaxrootNum))
         obj%stem2stem(:, :) = 0
         obj%leaf2stem(:, :) = 0
         obj%root2stem(:, :) = 0
         obj%root2root(:, :) = 0

         ! set mainstem
         do i = 1, obj%ms_node

            call obj%stem(i)%init()
            call obj%stem(i)%resize( &
               x=obj%ms_width, &
               y=obj%ms_width, &
               z=obj%ms_length/dble(obj%ms_node) &
               )
            call obj%stem(i)%rotate( &
               x=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)), &
               y=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)), &
               z=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)) &
               )
         end do

         do i = 1, obj%ms_node - 1
            call obj%stem(i + 1)%connect("=>", obj%stem(i))
            obj%stem2stem(i + 1, i) = 1
         end do

         ! set branches
         k = obj%ms_node
         do i = 1, size(obj%br_node)
            do j = 1, obj%br_node(i)
               k = k + 1
               call obj%stem(k)%init()
               call obj%stem(k)%resize( &
                  x=obj%ms_width, &
                  y=obj%ms_width, &
                  z=obj%ms_length/dble(obj%ms_node) &
                  )

               call obj%stem(k)%rotate( &
                  x=radian(random%gauss(mu=obj%br_angle_ave(j), sigma=obj%br_angle_sig(j))), &
                  y=0.0d0, &
                  z=radian(360.0d0*random%random()) &
                  )

               if (j == 1) then
                  call obj%stem(k)%connect("=>", obj%stem(obj%br_from(i)))
                  obj%stem2stem(k, obj%br_from(i)) = 1
               else
                  call obj%stem(k)%connect("=>", obj%stem(k - 1))
                  obj%stem2stem(k, k - 1) = 1
               end if

            end do
         end do

         ! peti and leaf
         num_stem_node = k
         num_leaf = 0
         do i = 1, num_stem_node
            ! ３複葉
            ! add peti
            num_stem_node = num_stem_node + 1
            call obj%stem(num_stem_node)%init()

            call obj%stem(num_stem_node)%resize( &
               x=random%gauss(mu=obj%peti_width_ave(i), sigma=obj%peti_width_sig(i)), &
               y=random%gauss(mu=obj%peti_width_ave(i), sigma=obj%peti_width_sig(i)), &
               z=random%gauss(mu=obj%peti_size_ave(i), sigma=obj%peti_size_sig(i)) &
               )

            call obj%stem(num_stem_node)%rotate( &
               x=radian(random%gauss(mu=obj%peti_angle_ave(i), sigma=obj%peti_angle_sig(i))), &
               y=0.0d0, &
               z=radian(360.0d0*random%random()) &
               )
            call obj%stem(num_stem_node)%connect("=>", obj%stem(i))
            obj%leaf2stem(num_stem_node, i) = 1

            ! add leaves
            do j = 1, 3
               num_leaf = num_leaf + 1
               call obj%leaf(num_leaf)%init()
               call obj%leaf(num_leaf)%resize( &
                  y=random%gauss(mu=obj%leaf_thickness_ave(i), sigma=obj%leaf_thickness_sig(i)), &
                  z=random%gauss(mu=obj%leaf_length_ave(i), sigma=obj%leaf_length_sig(i)), &
                  x=random%gauss(mu=obj%leaf_width_ave(i), sigma=obj%leaf_width_sig(i)) &
                  )
               call obj%leaf(num_leaf)%rotate( &
                  x=radian(random%gauss(mu=obj%leaf_angle_ave(i), sigma=obj%leaf_angle_sig(i))), &
                  y=0.0d0, &
                  z=radian(random%random()*360.0d0) &
                  )
               call obj%leaf(num_leaf)%connect("=>", obj%stem(num_stem_node))
               obj%leaf2stem(num_leaf, num_stem_node) = 1
            end do

         end do

         ! set mainroot
         do i = 1, obj%mr_node

            call obj%root(i)%init()
            call obj%root(i)%resize( &
               x=obj%mr_width, &
               y=obj%mr_width, &
               z=obj%mr_length/dble(obj%mr_node) &
               )
            call obj%root(i)%rotate( &
               x=radian(random%gauss(mu=obj%mr_angle_ave, sigma=obj%mr_angle_sig)), &
               y=radian(random%gauss(mu=obj%mr_angle_ave, sigma=obj%mr_angle_sig)), &
               z=radian(random%gauss(mu=obj%mr_angle_ave, sigma=obj%mr_angle_sig)) &
               )
         end do

         do i = 1, obj%mr_node - 1
            if (i == 1) then
               call obj%root(1)%connect("=>", obj%stem(1))
               obj%root2stem(1, 1) = 1
            end if
            call obj%root(i + 1)%connect("=>", obj%root(i))
            obj%root2root(i + 1, i) = 1
         end do

         ! set branches
         k = obj%mr_node
         do i = 1, size(obj%brr_node)
            do j = 1, obj%brr_node(i)
               k = k + 1
               call obj%root(k)%init()
               call obj%root(k)%resize( &
                  x=obj%mr_width, &
                  y=obj%mr_width, &
                  z=obj%mr_length/dble(obj%mr_node) &
                  )

               call obj%root(k)%rotate( &
                  x=radian(random%gauss(mu=obj%brr_angle_ave(j), sigma=obj%brr_angle_sig(j))), &
                  y=0.0d0, &
                  z=radian(360.0d0*random%random()) &
                  )

               if (j == 1) then
                  call obj%root(k)%connect("=>", obj%root(obj%brr_from(i)))
                  obj%root2root(k, obj%brr_from(i)) = 1
               else
                  call obj%root(k)%connect("=>", obj%root(k - 1))
                  obj%root2root(k, k - 1) = 1
               end if

            end do
         end do

         obj%stage = "V"//str(obj%ms_node)
         return
      else
         ! create leaf, root, stem
         allocate (obj%leaf(obj%MaxLeafNum))
         allocate (obj%root(obj%MaxrootNum))
         allocate (obj%stem(obj%MaxstemNum))
         allocate (obj%stem2stem(obj%MaxstemNum, obj%MaxstemNum))
         allocate (obj%leaf2stem(obj%MaxstemNum, obj%MaxLeafNum))
         allocate (obj%root2stem(obj%MaxrootNum, obj%MaxstemNum))
         allocate (obj%root2root(obj%MaxrootNum, obj%MaxrootNum))

         !allocate(obj%struct%NodCoord(4,3) )
         !allocate(obj%struct%ElemNod(3,2) )
         !allocate(obj%struct%ElemMat(3) )
         ! 子葉結節部=(0,0,0)
         !obj%struct%NodCoord(1,1:3) = 0.0d0
         call obj%leaf(1)%init(obj%leafconfig)
         call obj%leaf(1)%rotate(x=radian(90.0d0), y=radian(90.0d0), z=radian(10.0d0))
         call obj%leaf(2)%init(obj%leafconfig)
         call obj%leaf(2)%rotate(x=radian(90.0d0), y=radian(90.0d0), z=radian(-10.0d0))

         call obj%stem(1)%init(obj%stemconfig)
         call obj%stem(1)%rotate(x=radian(40.0d0))

         call obj%stem(2)%init(obj%stemconfig)
         call obj%stem(2)%rotate(x=radian(80.0d0))

         call obj%root(1)%init(obj%rootconfig)
         call obj%root(1)%fix(x=0.0d0, y=0.0d0, z=0.0d0)
         call obj%root(1)%rotate(x=radian(-60.0d0))

         call obj%leaf(1)%connect("=>", obj%stem(1))
         obj%leaf2stem(1, 1) = 1

         call obj%leaf(2)%connect("=>", obj%stem(1))
         obj%leaf2stem(2, 1) = 1

         call obj%stem(2)%connect("=>", obj%stem(1))
         obj%stem2stem(2, 1) = 1

         call obj%root(1)%connect("=>", obj%stem(1))
         obj%root2stem(1, 1) = 1

         obj%stage = "VE"
         ! 初生葉結節部
         !obj%struct%NodCoord(2,1) = 0.0d0
         !obj%struct%NodCoord(2,2) = 0.0d0
         !obj%struct%NodCoord(2,3) = 1.0d0/20.0d0*obj%seed_height
         ! 地際部
         !obj%struct%NodCoord(3,1) = 1.0d0/4.0d0*obj%seed_length
         !obj%struct%NodCoord(3,2) = 0.0d0
         !obj%struct%NodCoord(3,3) = -1.0d0/3.0d0*obj%seed_height
         ! 根冠
         !obj%struct%NodCoord(4,1) = 1.0d0/2.0d0*obj%seed_length
         !obj%struct%NodCoord(4,2) = 0.0d0
         !obj%struct%NodCoord(4,3) = -1.0d0/2.0d0*obj%seed_height

         ! 子葉-初生葉節
         !obj%struct%ElemNod(1,1) = 1
         !obj%struct%ElemNod(1,2) = 2
         ! 地際-子葉節
         !obj%struct%ElemNod(2,1) = 3
         !obj%struct%ElemNod(2,2) = 1
         ! 地際-根冠節
         !obj%struct%ElemNod(3,1) = 3
         !obj%struct%ElemNod(3,2) = 4

         ! 子葉-初生葉節 stem: 1
         !obj%struct%ElemMat(1) = 1
         ! 地際-子葉節 stem: 1
         !obj%struct%ElemMat(2) = 1
         ! 地際-根冠節 primary root: -1
         !obj%struct%ElemMat(3) = -1

         ! FEメッシュを生成
         ! 領域を確保
         !    n = input(default=80,option=max_leaf_num)
         !    allocate(obj%leaf_list(n) )
         !    n = input(default=80,option=max_stem_num)
         !    allocate(obj%stem_list(n) )
         !    n = input(default=80,option=max_root_num)
         !    allocate(obj%root_list(n) )
         !
         !    ! 子葉のメッシュを生成
         !    call obj%leaf_list(1)%create(meshtype="Sphere3D",x_num=10,y_num=10,z_num=10,&
         !        x_len=obj%seed_length,y_len=obj%seed_width,z_len=obj%seed_height)
         !    call obj%leaf_list(1)%move(x=0.0d0,y=-0.50d0*obj%seed_width,z=-0.50d0*obj%seed_height)
         !
         !    call obj%leaf_list(2)%create(meshtype="Sphere3D",x_num=10,y_num=10,z_num=10,&
         !        x_len=obj%seed_length,y_len=obj%seed_width,z_len=obj%seed_height)
         !    call obj%leaf_list(2)%rotate(x=radian(180.0d0) )
         !    call obj%leaf_list(2)%move(x=0.0d0,y=-0.50d0*obj%seed_width,z=-0.50d0*obj%seed_height)
         !
         !
         !
         !    ! 子葉-初生葉節のメッシュを生成
         !    rot(:) = 0.0d0
         !    call obj%stem_list(1)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
         !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
         !    ! 節基部の節点ID
         !    node_id = obj%struct%ElemNod(1,1)
         !    ! 節先端部の節点ID
         !    node_id2= obj%struct%ElemNod(1,2)
         !    ! 節基部の位置ベクトル
         !    loc(:) = obj%struct%NodCoord( node_id  ,:)
         !    ! 節先端部までの方向ベクトル
         !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)
         !
         !    ! structの構造データにメッシュデータを合わせる。
         !    print *, obj%stem_list(1)%Mesh%BottomElemID
         !    print *, obj%stem_list(1)%Mesh%TopElemID
         !
         !    elemid = obj%stem_list(1)%Mesh%BottomElemID
         !    node_id = obj%stem_list(1)%Mesh%ElemNod(elemID,1)
         !    meshloc(:) = obj%stem_list(1)%Mesh%NodCoord(node_id,:)
         !
         !    elemid = obj%stem_list(1)%Mesh%TopElemID
         !    node_id = obj%stem_list(1)%Mesh%ElemNod(elemID,1)
         !    meshvec(:) = obj%stem_list(1)%Mesh%NodCoord(node_id,:)-meshloc(:)

         !print *, "loc",loc
         !print *, "meshloc",meshloc
         !print *, "vec",vec
         !print *, "meshvec",meshvec

         !    ! 節中央を原点へ
         !    call obj%stem_list(1)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0)
         !
         !    print *, "loc",loc
         !    print *, "vec",vec
         !    print *, "rot",rot
         !    zaxis(:)=0.0d0
         !    zaxis(3)=obj%seed_length/5.0d0
         !    rot(:) = angles(zaxis,vec)
         !    call obj%stem_list(1)%move(x=loc(1),y=loc(2),z=loc(3) )
         !    call obj%stem_list(1)%rotate(x=0.0d0,y=0.0d0,z=0.0d0 )
    !!
         !
    !!
         !
         !
         !    ! 地際-子葉節のメッシュを生成
         !    rot(:) = 0.0d0
         !    call obj%stem_list(2)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
         !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
         !    ! 節基部の節点ID
         !    node_id = obj%struct%ElemNod(2,1)
         !    ! 節先端部の節点ID
         !    node_id2= obj%struct%ElemNod(2,2)
         !    ! 節基部の位置ベクトル
         !    loc(:) = obj%struct%NodCoord( node_id  ,:)
         !    ! 節先端部までの方向ベクトル
         !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)
         !    ! 節中央を原点へ
         !    call obj%stem_list(2)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0,&
         !        z=-obj%seed_length/8.0d0)
         !    zaxis(:)=0.0d0
         !    zaxis(3)=obj%seed_length/5.0d0
         !    rot(:) = angles(zaxis,vec)
         !    print *, "loc",loc
         !    print *, "vec",vec
         !    print *, "rot",rot
         !    !call obj%stem_list(2)%rotate(x=rot(1),y=rot(2),z=rot(3) )
         !    call obj%stem_list(2)%move(x=loc(1),y=loc(2),z=loc(3) )
         !
         !
         !
         !    ! 地際-根冠節のメッシュ生成
         !    rot(:) = 0.0d0
         !    call obj%root_list(1)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
         !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
         !    ! 節基部の節点ID
         !    node_id = obj%struct%ElemNod(3,1)
         !    ! 節先端部の節点ID
         !    node_id2= obj%struct%ElemNod(3,2)
         !    ! 節基部の位置ベクトル
         !    loc(:) = obj%struct%NodCoord( node_id  ,:)
         !    ! 節先端部までの方向ベクトル
         !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)
         !    ! 節基部へ移動
         !    call obj%root_list(1)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0,&
         !        z=-obj%seed_length/8.0d0)
         !    call obj%root_list(1)%move(x=loc(1),y=loc(2),z=loc(3) )
         !    zaxis(:)=0.0d0
         !    zaxis(3)=obj%seed_length/5.0d0
         !    rot(:) = angles(zaxis,vec)
         !    !call obj%root_list(1)%rotate(x=rot(1),y=rot(2),z=rot(3) )
         !    print *, "loc",loc
         !    print *, "vec",vec
         !    print *, "rot",rot
      end if

      ! ここからレガシーモード
      if (present(regacy)) then
         if (regacy .eqv. .true.) then
            obj%Stage = "VE"
            if (present(FileName)) then
               fn = FileName
            else
               fn = "untitled"
            end if

            loc(:) = 0.0d0

            if (present(x)) then
               loc(1) = x
            end if

            if (present(y)) then
               loc(2) = y
            end if

            if (present(z)) then
               loc(3) = z
            end if

            if (present(location)) then
               loc(:) = location(:)
            end if

            ! initialize RootSystem and NodeSystem
            if (.not. allocated(obj%RootSystem)) then
               allocate (obj%RootSystem(input(default=1000, option=max_PlantNode_num)))
               obj%num_of_root = 1
            end if
            if (.not. allocated(obj%NodeSystem)) then
               allocate (obj%NodeSystem(input(default=1000, option=max_PlantNode_num)))
               obj%num_of_node = 1
            end if

            ! setup seed
            if (Variety == "Tachinagaha" .or. Variety == "tachinagaha") then
               call obj%Seed%init(mass=mass, width1=9.70d0, width2=8.20d0, &
                                  width3=7.70d0, &
                                  water_content=water_content, radius=radius, location=loc)
               call obj%Seed%createMesh(FileName=fn//".stl", &
                                        ElemType="Tetrahedra")

               call obj%Seed%convertMeshType(Option="TetraToHexa")

            else
               print *, "Variety name :: is not implemented."
               stop
            end if

            ! setup primary node (plumule)
            call obj%NodeSystem(1)%init(Stage=obj%Stage, &
                                        Plantname="Rice", location=loc)

            ! setup primary node (radicle))
            MaxThickness = input(default=0.20d0, &
                                 option=PlantRoot_diameter_per_seed_radius)*obj%Seed%radius
            Maxwidth = input(default=0.20d0, &
                             option=PlantRoot_diameter_per_seed_radius)*obj%Seed%radius
            call obj%RootSystem(1)%init(Plantname="Rice", &
                                        Stage=obj%Stage, MaxThickness=MaxThickness, Maxwidth=Maxwidth, location=loc)

            obj%time = 0.0d0
            return
         end if
      end if

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

! ########################################
   subroutine growRice(obj, dt, light, air, temp)
      class(Rice_), intent(inout) :: obj
      type(Light_), optional, intent(inout) :: light
      type(air_), optional, intent(in) :: air
      real(real64), optional, intent(in) :: temp
      real(real64), intent(in) :: dt! time-interval
      real(real64) :: ac_temp ! time-interval
      integer(int32) :: i

      obj%dt = dt

      ! 光量子量を計算
      call obj%laytracing(light=light)

      ! 光合成量を計算
      do i = 1, size(obj%Leaf)
         if (obj%Leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%leaf(i)%photosynthesis(dt=dt, air=air)
         end if
      end do

      ! シンクソース輸送を計算
      call obj%SinkSourceFlow()

      ! ソースの消耗、拡散を計算
      !call obj%source2sink()

      ! 伸長を計算
      !call obj%extention()

      ! 分化を計算、構造の更新
      !call obj%development()

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

   subroutine SinkSourceFlowRice(obj)
      class(Rice_), intent(inout) :: obj
      type(DiffusionEq_) :: DiffusionEq

      !DiffusionEq%femdomain => obj%femdomain

   end subroutine

! ########################################
   subroutine WaterAbsorptionRice(obj, temp, dt)
      class(Rice_), intent(inout) :: obj
      real(real64), intent(in) :: temp, dt
      real(real64) :: a, b, c, d, AA, BB, w1max, w2max, w3max, time
      real(real64) :: x_rate, y_rate, z_rate, wx, wy, wz

      obj%time = obj%time + dt

      ! tested by tachinagaha, 2019
      a = 0.00910d0
      b = -1.76450d0
      c = 3.32E-04
      d = -0.0905180d0
      AA = a*temp + b
      !BB=c*exp(d*temp)
      BB = c*temp + d
      ! width1 becomes 1.7 times, width2 becomes 1.2, width3 becomes 1.1
      w1max = 1.70d0
      w2max = 1.20d0
      w3max = 1.10d0
      obj%seed%width1 = obj%seed%width1_origin*(w1max - AA*exp(-BB*obj%time))
      obj%seed%width2 = obj%seed%width2_origin*(w2max - AA*exp(-BB*obj%time))
      obj%seed%width3 = obj%seed%width3_origin*(w3max - AA*exp(-BB*obj%time))

      ! linear model; it should be changed in near future.
      if (obj%time > 60.0d0*6.0d0) then
         obj%seed%width2 = obj%seed%width2_origin*(w2max)
         obj%seed%width3 = obj%seed%width3_origin*(w3max)
      else
         obj%seed%width2 = obj%seed%width2_origin + obj%seed%width2_origin*(w2max - 1.0d0)*(obj%time)/(60.0d0*6.0d0)
         obj%seed%width3 = obj%seed%width3_origin + obj%seed%width3_origin*(w3max - 1.0d0)*(obj%time)/(60.0d0*6.0d0)
      end if

      wx = maxval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 1)) - minval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 1))
      wy = maxval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 2)) - minval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 2))
      wz = maxval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 3)) - minval(obj%Seed%FEMDomain%Mesh%NodCoord(:, 3))
      print *, wx, wy, wz
      x_rate = 1.0d0/wx
      y_rate = 1.0d0/wy
      z_rate = 1.0d0/wz
      call obj%Seed%FEMDomain%resize(x_rate=x_rate, y_rate=y_rate, z_rate=z_rate)
      x_rate = obj%seed%width1
      y_rate = obj%seed%width2
      z_rate = obj%seed%width3
      call obj%Seed%FEMDomain%resize(x_rate=x_rate, y_rate=y_rate, z_rate=z_rate)

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

! ########################################
   subroutine exportRice(obj, FilePath, FileName, SeedID, withSTL, withMesh)
      class(Rice_), intent(inout) :: obj
      character(*), optional, intent(in) :: FilePath
      character(*), intent(in) :: FileName
      integer(int32), optional, intent(inout) :: SeedID
      logical, optional, intent(in) :: withSTL, withMesh
      integer(int32) :: i, itr

      itr = SeedID
      ! if seed exists => output
      if (obj%Seed%num_of_seed >= 0) then
         if (present(withSTL)) then
            if (withSTL .eqv. .true.) then
               call obj%Seed%export(FileName=FileName, SeedID=itr, extention=".stl")
            end if
         end if
         if (present(withMesh)) then
            if (withMesh .eqv. .true.) then
               call obj%Seed%export(FileName=FileName, SeedID=itr, extention=".pos")
            end if
         end if

         if (present(FilePath)) then
            call obj%Seed%export(FileName=FilePath//"/seed.geo", SeedID=itr)
         else
            call obj%Seed%export(FileName=FileName, SeedID=itr)
         end if
      end if

      itr = itr + 1
      ! export NodeSystem
      do i = 1, size(obj%NodeSystem)

         if (present(FilePath)) then
            call obj%NodeSystem(i)%export(FileName=FilePath//"/Node.geo", objID=itr)
         else
            call obj%NodeSystem(i)%export(FileName=FileName//"_Node.geo", objID=itr)
         end if
         if (i == obj%num_of_node) then
            exit
         end if
      end do

      ! export RootSystem
      do i = 1, size(obj%RootSystem)

         if (present(FilePath)) then
            call obj%RootSystem(i)%export(FileName=FilePath//"/Root.geo", RootID=itr)
         else
            call obj%RootSystem(i)%export(FileName=FileName//"_Root.geo", RootID=itr)
         end if
         if (i == obj%num_of_root) then
            exit
         end if
      end do
      SeedID = itr

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

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

! ########################################
!subroutine initRice(obj,growth_habit,Max_Num_of_Node)
!    class(Rice_) :: obj
!    character(*),optional,intent(in) :: growth_habit
!    integer(int32),optional,intent(in)::Max_Num_of_Node
!    integer(int32) ::n
!
!    if(present(growth_habit) )then
!        obj%growth_habit=growth_habit
!    else
!        obj%growth_habit="determinate"
!    endif
!
!    obj%growth_stage="VE"
!
!    n=input(default=100,option=Max_Num_of_Node)
!
!    allocate(obj%NodeSystem(n))
!    obj%NumOfNode=0
!    obj%NumOfRoot=0
!
!    ! set an initial node and root
!    ! two leaves, one root.
!
!    call obj%AddNode()
!
!end subroutine
!! ########################################
!
!
!
!
!
!
!! ########################################
!subroutine AddNodeRice(obj,SizeRatio)
!    class(Rice_),intent(inout)::obj
!    real(real64),optional,intent(in)::SizeRatio
!    real(real64) :: magnif
!
!    magnif=input(default=1.0d0,option=SizeRatio)
!    obj%NumOfNode=obj%NumOfNode+1
!
!    ! add leaves
!    if(obj%NumOfNode==1 .or. obj%NumOfNode==2)then
!        allocate(obj%NodeSystem(obj%NumOfNode)%leaf(2) )
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=3.0d0*magnif,width=2.0d0*magnif)
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=3.0d0*magnif,width=2.0d0*magnif)
!    else
!        allocate(obj%NodeSystem(obj%NumOfNode)%leaf(3) )
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=4.0d0*magnif,width=2.0d0*magnif)
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=4.0d0*magnif,width=2.0d0*magnif)
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=4.0d0*magnif,width=2.0d0*magnif)
!    endif
!
!    ! add stem
!    if(obj%NumOfNode==1 .or. obj%NumOfNode==2)then
!        allocate(obj%NodeSystem(obj%NumOfNode)%Stem(1) )
!        call obj%NodeSystem(obj%NumOfNode)%leaf(1)%init(thickness=0.10d0*magnif,length=3.0d0*magnif,width=2.0d0*magnif)
!    endif
!
!    ! add Peti
!    if(obj%NumOfNode==1 .or. obj%NumOfNode==2)then
!        allocate(obj%NodeSystem(obj%NumOfNode)%Peti(1) )
!        call obj%NodeSystem(obj%NumOfNode)%Peti(1)%init(thickness=0.10d0*magnif,length=3.0d0*magnif,width=2.0d0*magnif)
!    endif
!
!end subroutine
!! ########################################
!

! ########################################
   subroutine showRice(obj, name)
      class(Rice_), intent(inout) :: obj
      character(*), intent(in)::name

      if (obj%struct%empty() .eqv. .true.) then
         print *, "Error :: showRice>> no structure is imported."
         return
      end if

      call obj%struct%export(name=name)

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

! ########################################
   function numleafRice(obj) result(ret)
      class(Rice_), intent(in) :: obj
      integer(int32) :: ret, i

      ret = 0
      do i = 1, size(obj%leaf_list)
         if (obj%leaf_list(i)%Mesh%empty() .eqv. .false.) then
            ret = ret + 1
         end if
      end do

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

! ########################################
   function numstemRice(obj) result(ret)
      class(Rice_), intent(in) :: obj
      integer(int32) :: ret, i

      ret = 0
      do i = 1, size(obj%stem_list)
         if (obj%stem_list(i)%Mesh%empty() .eqv. .false.) then
            ret = ret + 1
         end if
      end do

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

! ########################################
   function numrootRice(obj) result(ret)
      class(Rice_), intent(in) :: obj
      integer(int32) :: ret, i

      ret = 0
      do i = 1, size(obj%root_list)
         if (obj%root_list(i)%Mesh%empty() .eqv. .false.) then
            ret = ret + 1
         end if
      end do

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

! ########################################
   subroutine gmshRice(obj, name)
      class(Rice_), intent(inout) :: obj
      character(*), intent(in) :: name
      integer(int32) :: i

      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%stem(i)%gmsh(name=name//"_stem"//str(i))
         end if
      end do

      do i = 1, size(obj%root)
         if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%root(i)%gmsh(name=name//"_root"//str(i))
         end if
      end do

      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%leaf(i)%gmsh(name=name//"_leaf"//str(i))
         end if
      end do

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

! ########################################
   subroutine mshRice(obj, name)
      class(Rice_), intent(inout) :: obj
      character(*), intent(in) :: name
      integer(int32) :: i

      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%stem(i)%msh(name=name//"_stem"//str(i))
         end if
      end do

      do i = 1, size(obj%root)
         if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%root(i)%msh(name=name//"_root"//str(i))
         end if
      end do

      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%leaf(i)%msh(name=name//"_leaf"//str(i))
         end if
      end do

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

! ########################################
   subroutine jsonRice(obj, name)
      class(Rice_), intent(inout) :: obj
      character(*), intent(in) :: name
      integer(int32) :: i, countnum
      type(IO_) :: f

      call f%open(name//".json")
      call f%write("{")
      countnum = 0
      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            countnum = countnum + 1
            call f%write('"'//"stem"//str(i)//'":')
            call obj%stem(i)%femdomain%json(name=name//"_stem"//str(i), fh=f%fh, endl=.false.)
         end if
      end do
      call f%write('"num_stem":'//str(countnum)//',')

      countnum = 0
      do i = 1, size(obj%root)
         if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then
            countnum = countnum + 1
            call f%write('"'//"root"//str(i)//'":')
            call obj%root(i)%femdomain%json(name=name//"_root"//str(i), fh=f%fh, endl=.false.)
         end if
      end do
      call f%write('"num_root":'//str(countnum)//',')

      countnum = 0
      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            countnum = countnum + 1
            call f%write('"'//"leaf"//str(i)//'":')
            call obj%leaf(i)%femdomain%json(name=name//"_leaf"//str(i), fh=f%fh, endl=.false.)
         end if
      end do
      call f%write('"num_leaf":'//str(countnum)//',')
      call f%write('"return_Rice":0')
      call f%write("}")
      call f%close()
   end subroutine
! ########################################

! ########################################
   subroutine stlRice(obj, name)
      class(Rice_), intent(inout) :: obj
      character(*), intent(in) :: name
      integer(int32) :: i

      !call execute_command_line("echo ' ' > "//name//".stl")
      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%stem(i)%stl(name=name//"_stem"//str(i))
            !call execute_command_line("cat "//name//"_stem"//str(i)//"_000001.stl >> "//name//".stl")
         end if
      end do

      do i = 1, size(obj%root)
         if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%root(i)%stl(name=name//"_root"//str(i))
            !call execute_command_line("cat "//name//"_root"//str(i)//"_000001.stl >> "//name//".stl")
         end if
      end do

      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%leaf(i)%stl(name=name//"_leaf"//str(i))
            !call execute_command_line("cat "//name//"_leaf"//str(i)//"_000001.stl >> "//name//".stl")
         end if
      end do

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

! ########################################
   subroutine moveRice(obj, x, y, z)
      class(Rice_), intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z
      integer(int32) :: i

      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%stem(i)%move(x=x, y=y, z=z)
         end if
      end do

      do i = 1, size(obj%root)
         if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%root(i)%move(x=x, y=y, z=z)
         end if
      end do

      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            call obj%leaf(i)%move(x=x, y=y, z=z)
         end if
      end do

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

! ########################################
   subroutine laytracingRice(obj, light)
      class(Rice_), intent(inout) :: obj
      type(Light_), intent(in) :: light
      real(real64), allocatable :: stemcenter(:, :), stemradius(:)
      real(real64), allocatable :: leafcenter(:, :), leafradius(:)
      real(real64), allocatable :: elemnodcoord(:, :), x(:), x2(:)
      real(real64) :: max_PPFD, r, rc, r0
      real(real64), parameter :: extinction_ratio = 100.0d0 ! ratio/m
      !real(real64),parameter :: radius_ratio = 0.01d0 ! radius_of_gauss_point/element_length
      type(IO_) :: f
      integer(int32) :: i, j, n, num_particle, k, l, nodeid, m, totcount

      max_PPFD = light%maxPPFD
      ! 総当りで、総遮蔽長を割り出す
      ! 茎は光を通さない、葉は透過率あり、空間は透過率ゼロ
      ! 要素中心から頂点への平均長さを半径に持ち、要素中心を中心とする球
      ! を考え、Layとの公差判定を行う。
      num_particle = 0
      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            num_particle = num_particle + size(obj%leaf(i)%femdomain%mesh%ElemNod, 1)
         end if
      end do
      allocate (leafcenter(num_particle, 3), leafradius(num_particle))
      leafcenter(:, :) = 0.0d0
      leafradius(:) = 0.0d0

      num_particle = 0
      do i = 1, size(obj%leaf)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            num_particle = num_particle + size(obj%stem(i)%femdomain%mesh%ElemNod, 1)
         end if
      end do
      allocate (stemcenter(num_particle, 3), stemradius(num_particle))
      stemcenter(:, :) = 0.0d0
      stemradius(:) = 0.0d0

      num_particle = 0

      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            n = size(obj%leaf(i)%femdomain%mesh%Elemnod, 2)
            m = size(obj%leaf(i)%femdomain%mesh%Nodcoord, 2)
            allocate (elemnodcoord(n, m))
            allocate (x(m))
            do j = 1, size(obj%leaf(i)%femdomain%mesh%elemnod, 1)
               do k = 1, size(obj%leaf(i)%femdomain%mesh%elemnod, 2)
                  nodeid = obj%leaf(i)%femdomain%mesh%elemnod(j, k)
                  elemnodcoord(k, :) = obj%leaf(i)%femdomain%mesh%Nodcoord(nodeid, :)
               end do
               num_particle = num_particle + 1
               do k = 1, size(elemnodcoord, 1)
                  do l = 1, size(elemnodcoord, 2)
                     leafcenter(num_particle, l) = &
                        +leafcenter(num_particle, l) &
                        + 1.0d0/dble(size(elemnodcoord, 1))*elemnodcoord(k, l)
                  end do
               end do
               do k = 1, size(elemnodcoord, 1)
                  x(:) = elemnodcoord(k, :)
                  x(:) = x(:) - leafcenter(num_particle, :)
                  if (k >= 2 .and. leafradius(num_particle) > sqrt(dot_product(x, x))) then
                     leafradius(num_particle) = sqrt(dot_product(x, x))
                  elseif (k == 1) then
                     leafradius(num_particle) = sqrt(dot_product(x, x))
                  else
                     cycle
                  end if
               end do
            end do
            deallocate (elemnodcoord)
            deallocate (x)
         end if
      end do

      num_particle = 0
      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then
            n = size(obj%stem(i)%femdomain%mesh%Elemnod, 2)
            m = size(obj%stem(i)%femdomain%mesh%Nodcoord, 2)
            allocate (elemnodcoord(n, m))
            allocate (x(m))
            do j = 1, size(obj%stem(i)%femdomain%mesh%elemnod, 1)
               do k = 1, size(obj%stem(i)%femdomain%mesh%elemnod, 2)
                  nodeid = obj%stem(i)%femdomain%mesh%elemnod(j, k)
                  elemnodcoord(k, :) = obj%stem(i)%femdomain%mesh%Nodcoord(nodeid, :)
               end do
               num_particle = num_particle + 1
               do k = 1, size(elemnodcoord, 1)
                  do l = 1, size(elemnodcoord, 2)
                     stemcenter(num_particle, l) = &
                        +stemcenter(num_particle, l) &
                        + 1.0d0/dble(size(elemnodcoord, 1))*elemnodcoord(k, l)
                  end do
               end do
               do k = 1, size(elemnodcoord, 1)
                  x(:) = elemnodcoord(k, :)
                  x(:) = x(:) - stemcenter(num_particle, :)
                  !最小半径で考える
                  if (k >= 2 .and. stemradius(num_particle) > sqrt(dot_product(x, x))) then
                     stemradius(num_particle) = sqrt(dot_product(x, x))
                  elseif (k == 1) then
                     stemradius(num_particle) = sqrt(dot_product(x, x))
                  else
                     cycle
                  end if
               end do
            end do
            deallocate (elemnodcoord)
            deallocate (x)
         end if
      end do

      ! DEBUG
      call f%open("leaf.txt")
      do i = 1, size(leafcenter, 1)
         write (f%fh, *) leafcenter(i, :)
      end do
      call f%close()

      call f%open("stem.txt")
      do i = 1, size(stemcenter, 1)
         write (f%fh, *) stemcenter(i, :)
      end do
      call f%close()

      allocate (x(3), x2(3))

      num_particle = 0
      totcount = 0
      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            ! 葉あり
            obj%leaf(i)%PPFD(:) = max_PPFD
            do j = 1, size(obj%leaf(i)%PPFD)
               totcount = totcount + 1
               num_particle = num_particle + 1
               ! それぞれの要素について、遮蔽particleを探索
               ! 茎：全減衰
               ! 葉：半減衰
               ! 簡単のため上からのみ
               ! x-yのみについて見て、上方かつx-y平面距離が半径以内で覆陰判定
               x(:) = leafcenter(num_particle, :)
               r0 = leafradius(num_particle)
               ! 枝による覆陰判定

               do k = 1, size(stemcenter, 1)
                  x2(:) = stemcenter(k, :)
                  r = stemradius(k)
                  rc = (x(1) - x2(1))**(2.0d0) + (x(2) - x2(2))**(2.0d0)
                  rc = sqrt(rc)
                  if (rc <= r0 + r .and. x(3) < x2(3)) then
                     ! 茎により覆陰されてる
                     obj%leaf(i)%PPFD(j) = 0.0d0
                     exit
                  end if
               end do
               if (obj%leaf(i)%PPFD(j) == 0.0d0) then
                  cycle
               end if

               do k = 1, size(leafcenter, 1)
                  ! もし自信だったら除外
                  if (totcount == k) then
                     cycle
                  end if

                  x2(:) = leafcenter(k, :)
                  r = leafradius(k)
                  rc = (x(1) - x2(1))**(2.0d0) + (x(2) - x2(2))**(2.0d0)
                  rc = sqrt(rc)
                  if (rc <= (r0 + r)/2.0d0 .and. x(3) < x2(3)) then
                     ! 茎により覆陰されてる
                     obj%leaf(i)%PPFD(j) = &
                        obj%leaf(i)%PPFD(j)*(1.0d0 - extinction_ratio*2.0d0*r)
                     if (obj%leaf(i)%PPFD(j) <= 0.0d0) then
                        obj%leaf(i)%PPFD(j) = 0.0d0
                     end if
                  end if
               end do

            end do
         end if
      end do

      call f%open("PPFD.txt")
      do i = 1, size(obj%leaf)
         if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then
            ! 葉あり
            do j = 1, size(obj%leaf(i)%PPFD, 1)
               write (f%fh, *) obj%leaf(i)%PPFD(j), "leaf_id: ", str(i), "elem_id: ", str(j)
            end do
         end if
      end do
      call f%close()

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

   subroutine addStemRice(obj, stemid, rotx, roty, rotz, json)
      class(Rice_), intent(inout) :: obj
      integer(int32), intent(in) :: stemid
      character(*), optional, intent(in) :: json
      real(real64), optional, intent(in) :: rotx, roty, rotz
      integer(int32) :: i

      ! add a stem after stem(stemid)
      do i = 1, size(obj%stem)
         if (obj%stem(i)%femdomain%mesh%empty() .eqv. .true.) then
            if (present(json)) then
               call obj%stem(i)%init(json)
               call obj%stem(i)%rotate(x=rotx, y=roty, z=rotz)
               call obj%stem(i)%connect("=>", obj%stem(stemid))
               return
            else
               call obj%stem(i)%init()
               call obj%stem(i)%rotate(x=rotx, y=roty, z=rotz)
               call obj%stem(i)%connect("=>", obj%stem(stemid))
               return
            end if
         else
            cycle
         end if
      end do

   end subroutine

end module
