module RiceClass
    use LeafClass
    use StemClass
    use RootClass
    use PanicleClass
    implicit none


    type :: Rice_internode_info_
        real(real64),allocatable :: FinalInterNodeLength(:)
        real(real64),allocatable :: FinalPetioleLength(:)
        real(real64),allocatable :: FinalLeafLength(:)
        real(real64),allocatable :: FinalLeafWidth(:)
    end type
    
    type :: Rice_NodeID_Branch_
        integer(int32),allocatable :: ID(:)
    contains 
        !procedure, public :: sync => syncRice_NodeID_Branch
    end type
    
    type :: Rice_
        integer(int32) :: version 
        integer(int32) :: num_shoot=0

        ! for version >= 2.0, this instace has all states.
        ! >>>>>
        type(Rice_),allocatable :: rice_shoots(:)
        ! <<<<<

        
        integer(int32) :: TYPE_STEM    = 1
        integer(int32) :: TYPE_LEAF    = 2
        integer(int32) :: TYPE_ROOT    = 3
        integer(int32) :: TYPE_PANICLE = 4

        ! Only for pre-processing >>>>
        real(real64)   :: mainstem_length
        real(real64)   :: mainstem_width
        integer(int32) :: mainstem_node

        real(real64)   :: mainroot_length
        real(real64)   :: mainroot_width
        integer(int32) :: mainroot_node

        integer(int32) :: num_branch_root
        integer(int32) :: num_branch_root_node

        real(real64) :: ms_angle_ave = 0.0d0
        real(real64) :: ms_angle_sig = 0.0d0

        integer(int32),allocatable :: Leaf_From(:)
        real(real64),allocatable :: leaf_curvature(:)

        real(real64),allocatable :: leaf_thickness_ave(:)
        real(real64),allocatable :: leaf_thickness_sig(:)

        real(real64),allocatable :: leaf_angle_ave_x(:)
        real(real64),allocatable :: leaf_angle_sig_x(:)
        real(real64),allocatable :: leaf_angle_ave_z(:)
        real(real64),allocatable :: leaf_angle_sig_z(:)

        real(real64),allocatable :: leaf_length_ave(:)
        real(real64),allocatable :: leaf_length_sig(:)
        real(real64),allocatable :: leaf_width_ave(:)
        real(real64),allocatable :: leaf_width_sig(:)
        ! Only for pre-processing <<<<

        
        !character(:),allocatable :: LeafSurfaceData

        ! objects generated from json file
        ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
        integer(int32) :: num_leaf = 0
        integer(int32) :: num_stem = 0
        integer(int32) :: num_panicle = 0
        integer(int32) :: num_root= 0
        type(Leaf_),allocatable :: Leaf(:)
        type(Stem_),allocatable :: Stem(:)
        type(Panicle_),allocatable :: Panicle(:)
        type(Root_),allocatable :: Root(:)
        integer(int32),allocatable :: leaf2stem(:,:)
        integer(int32),allocatable :: stem2stem(:,:)
        integer(int32),allocatable :: panicle2stem(:,:)
        integer(int32),allocatable :: root2stem(:,:)
        integer(int32),allocatable :: root2root(:,:)
        ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

        ! integer(int32),allocatable :: NodeID_MainStem(:)
        ! type(Rice_NodeID_Branch_),allocatable :: NodeID_Branch(:)
        
        logical ::  inLoop = .false.
        real(real64) :: hours = 0.0d0
        
        ! growth simulation
        real(real64) :: FullyExpanded_stem_threshold = 0.10d0
        integer(int32) :: MaxBranchNum = 20
        type(Rice_internode_info_),allocatable :: InterNodeInfo(:)
        real(real64) :: default_Leaf_growth_ratio = 1.0d0/3.0d0
        real(real64) :: default_Stem_growth_ratio = 1.0d0/3.0d0
        integer(int32),allocatable :: MainStem_num_branch(:)
        real(real64) :: apical_dominance_distance = 1.0d0
        

        ! シミュレータ
        type(ContactMechanics_) :: contact
        real(real64) :: Gravity_acceralation = 9.810d0
        real(real64) :: PenaltyParameter = 100000.0d0
        logical :: GaussPointProjection = .false.

        ! setting
        integer(int32) :: stem_division(1:3) = [10, 10, 10]
        integer(int32) :: leaf_division(1:3) = [10, 10, 10]
        integer(int32) :: panicle_division(1:3) = [2, 2, 300]
        

    contains
        procedure,public :: init => createRice
        procedure,public :: create => createRice
        procedure,public :: createShoot => createRiceShoot
        procedure,public :: msh => mshRice
        procedure,public :: vtk => vtkRice
        procedure,public :: stl => stlRice
        procedure,public :: json => jsonRice
        procedure,public :: update => updateRice
        
        ! Editor
        procedure,public :: remove => removeRice
        procedure,public :: rotate => rotateRice
        procedure,public :: move => moveRice

        procedure,public :: getElementList => getElementListRice
        ! Info
        !procedure,public :: properties => propertiesRice
        ! number of subdomain
        procedure,public :: ns => nsRice
        ! number of element
        procedure,public :: ne => neRice
        ! number of points
        procedure,public :: nn => nnRice
        ! MemSize
        procedure,public :: checkMemoryRequirement => checkMemoryRequirementRice
        ! number of stem, leaf ...etc.
        procedure,public :: numleaf => numleafRice
        procedure,public :: numstem => numstemRice
        procedure,public :: numroot => numrootRice
        procedure,public :: numPanicle => numPanicleRice
        ! get pointers
        procedure,public :: getFEMDomainPointers => getFEMDomainPointersRice
        ! check data
        procedure,public :: checkYoungModulus => checkYoungModulusRice
        procedure,public :: checkPoissonRatio => checkPoissonRatioRice
        procedure,public :: checkDensity => checkDensityRice
        ! get data
        procedure,public :: getYoungModulus => getYoungModulusRice 
        procedure,public :: getPoissonRatio => getPoissonRatioRice
        procedure,public :: getDensity => getDensityRice
        
        procedure,public :: getYoungModulusField => getYoungModulusFieldRice
        procedure,public :: getPoissonRatioField => getPoissonRatioFieldRice
        procedure,public :: getDensityField => getDensityFieldRice
        procedure,public :: getStressField => getStressFieldRice
        
        ! alternative setters
        procedure,public :: setYoungModulus => setYoungModulusRice
        procedure,public :: setPoissonRatio => setPoissonRatioRice
        procedure,public :: setDensity => setDensityRice

        ! simulator
        procedure,public :: getEigenMode => getEigenModeRice
        procedure,public :: getDisplacement => getDisplacementRice
        procedure,public :: deform => deformRice

        ! export eigen modes and frequency
        procedure,public :: export_eig => export_eigRice
        
    end type

    interface operator(+)
        module procedure add_rice_to_rice
    end interface

    interface assignment(=)
        module procedure substitute_rice_into_rice
    end interface
contains

! #############################################################
recursive subroutine createRice(this,config,debug,output_detail_config)
    class(Rice_),intent(inout) :: this
    type(Rice_),allocatable :: rice_shoots(:)
    type(Rice_) :: this_shoot
    character(*),intent(in) :: config
    character(*),optional,intent(in) :: output_detail_config
    character(:),allocatable :: line
    logical,optional,intent(in) :: debug
    logical :: debug_log
    type(IO_) :: Riceconfig,input_dc,output_dc
    type(Random_) :: random
    type(Math_) :: math
    integer(int32)::i,n,j,k,num_leaf,num_stem_node,num_branch_branch,cpid,shoot_idx
    real(real64) :: x_A(1:3),rx,ry,angle,plot_angle_ave,plot_angle_sig


    if(present(output_detail_config) )then
        call output_dc%open(output_detail_config,"a")
        call output_dc%write("{")
    endif

    debug_log = input(default=.false.,option=debug)
    cpid = 0
    
    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
    
    this%version = fint(Riceconfig%parse_json(config,to_list("Version")))
    if(this%version == 2)then
        

        this%num_shoot = fint(Riceconfig%parse_json(config,to_list("num_shoot")))
        
        allocate(rice_shoots(this%num_shoot))
        do shoot_idx=1,this%num_shoot
            
            
            angle = random%random()*2.0d0*math%pi
            rx = random%random()*freal(Riceconfig%parse_json(config,to_list("plot_radius_x") ) )
            ry = random%random()*freal(Riceconfig%parse_json(config,to_list("plot_radius_y") ) )
            plot_angle_ave = radian(freal(Riceconfig%parse_json(config,to_list("plot_angle_ave") ) ))
            plot_angle_sig = radian(freal(Riceconfig%parse_json(config,to_list("plot_angle_sig") ) ))
    

            if(present(output_detail_config) )then
                call output_dc%open(output_detail_config,"a")
                call output_dc%write('"shoot_idx_'+str(shoot_idx)+'rx'+'":'+str(rx)+',')
                call output_dc%write('"shoot_idx_'+str(shoot_idx)+'ry'+'":'+str(ry)+',')
                call output_dc%write('"shoot_idx_'+str(shoot_idx)+'angle'+'":'+str(angle)+',')
            endif
            !実装途中!!@20240624

             
            call rice_shoots(shoot_idx)%createShoot(config=config,ShootIdx=shoot_idx,debug=debug)
            call rice_shoots(shoot_idx)%rotate(x=random%gauss(mu=plot_angle_ave,sigma=plot_angle_sig))
            
            call rice_shoots(shoot_idx)%move(x=rx*cos(angle),y=ry*sin(angle))
            call rice_shoots(shoot_idx)%rotate(z=random%random()*2.0d0*math%pi)




            if(shoot_idx==1)then
                this_shoot = rice_shoots(1)
            else
                this_shoot = this_shoot + rice_shoots(shoot_idx)
            endif
        enddo 
        this = this_shoot

        this%rice_shoots = rice_shoots
        ! integrate rice_shoots to a rice object
        !call rice%add(rice_shoots)
        return
    endif


    this%mainstem_length = freal(Riceconfig%parse(config,key1="Mainstem",key2="Length"))
    this%mainstem_width = freal(Riceconfig%parse(config,key1="Mainstem",key2="Width"))
    this%mainstem_node = fint(Riceconfig%parse(config,key1="Mainstem",key2="Node"))
    this%ms_angle_ave = freal(Riceconfig%parse(config,key1="Mainstem",key2="ms_angle_ave"))
    this%ms_angle_sig = freal(Riceconfig%parse(config,key1="Mainstem",key2="ms_angle_sig"))


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
    
    ! get number of leaf
    this%num_leaf=1
    do
        if(this%num_leaf ==  this%mainstem_node)then
            this%num_leaf = this%num_leaf -1
            exit
        endif

        line = Riceconfig%parse(config,key1="Leaf_"//str(this%num_leaf)//"_",key2="From")
        
        if(len(trim(line))==0)then
            this%num_leaf = this%num_leaf -1
            exit
        else
            this%num_leaf = this%num_leaf  + 1
            cycle
        endif
        
    enddo

    

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    allocate(this%leaf_curvature(this%num_leaf))
    allocate(this%leaf_thickness_ave(this%num_leaf))
    allocate(this%leaf_thickness_sig(this%num_leaf))
    allocate(this%leaf_angle_ave_x(this%num_leaf))
    allocate(this%leaf_angle_sig_x(this%num_leaf))
    allocate(this%leaf_angle_ave_z(this%num_leaf))
    allocate(this%leaf_angle_sig_z(this%num_leaf))
    allocate(this%leaf_length_ave(this%num_leaf))
    allocate(this%leaf_length_sig(this%num_leaf))
    allocate(this%leaf_width_ave(this%num_leaf))
    allocate(this%leaf_width_sig(this%num_leaf))
    allocate(this%leaf_From(this%num_leaf))
    !allocate(this%leaf_Length(this%num_leaf))
    !allocate(this%leaf_Width(this%num_leaf))


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    do i=1,this%num_leaf
        this%leaf_From(i)= fint(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="From"))
        !this%leaf_Length(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="Length"))
        !this%leaf_Width(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="Width"))
        this%leaf_curvature(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_curvature"))
        this%leaf_thickness_ave(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_thickness_ave"))
        this%leaf_thickness_sig(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_thickness_sig"))
        this%leaf_angle_ave_x(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_ave_x"))
        this%leaf_angle_sig_x(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_sig_x"))
        this%leaf_angle_ave_z(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_ave_z"))
        this%leaf_angle_sig_z(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_sig_z"))
        this%leaf_length_ave(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_length_ave"))
        this%leaf_length_sig(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_length_sig"))
        this%leaf_width_ave(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_width_ave"))
        this%leaf_width_sig(i)= freal(Riceconfig%parse(&
            config,key1="Leaf_"//str(i)//"_",key2="leaf_width_sig"))
    enddo


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
!    this%mainroot_length = freal(Riceconfig%parse(config,key1="Mainroot",key2="Length"))
!    this%mainroot_width = freal(Riceconfig%parse(config,key1="Mainroot",key2="Width"))
!    this%mainroot_node = fint(Riceconfig%parse(config,key1="Mainroot",key2="Node"))



    ! get number of branch && number of node
!    this%num_branch_root=1
!    this%num_branch_root_node=0
!    do 
!        line = Riceconfig%parse(config,key1="Branchroot_"//str(this%num_branch_root)//"_",key2="Node" )
!        if(len(trim(line))==0)then
!            this%num_branch_root = this%num_branch_root -1
!            exit
!        else
!            this%num_branch_root = this%num_branch_root  + 1
!            this%num_branch_root_node = this%num_branch_root_node + fint(line)
!            cycle
!        endif
!    enddo


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
    this%num_stem = this%mainstem_node 
    !this%num_root =this%num_branch_root_node + this%mainroot_node

    !allocate(this%leaf_list(this%num_leaf))
    !allocate(this%stem_list(this%num_stem))
    !allocate(this%root_list(this%num_root))

    allocate(this%leaf(this%num_leaf))
    allocate(this%stem(this%num_stem))
    !allocate(this%root(this%num_root))



    this%leaf2stem = zeros( this%num_leaf , this%num_stem) 
    this%stem2stem = zeros( this%num_stem, this%num_stem)
    this%panicle2stem = zeros( this%num_panicle, this%num_stem) 
    !this%root2stem = zeros( this%num_root , this%num_stem) 
    !this%root2root = zeros( this%num_root , this%num_root ) 

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    ! set mainstem
    do i=1,this%mainstem_node

        call this%stem(i)%init(&
            x_num = this%stem_division(1),&
            y_num = this%stem_division(2),&
            z_num = this%stem_division(3) &
            )
             
        call this%stem(i)%resize(&
            x = this%mainstem_width, &
            y = this%mainstem_width, &
            z = this%mainstem_length/dble(this%mainstem_node) &
            )
        call this%stem(i)%rotate(&
            x = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig)),  &
            y = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig)),  &
            z = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig))   &
            )
    enddo

    do i=1,this%mainstem_node-1
        call this%stem(i+1)%connect("=>",this%stem(i))
        this%stem2stem(i+1,i) = 1
    enddo
    

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    !set leaf
    num_leaf = 0
    do i=1,this%num_leaf
        ! 1葉/1節
        ! add leaves
        
        num_leaf=num_leaf+1
        
        call this%leaf(num_leaf)%init(species=PF_RICE,&
            x_num = this%leaf_division(1),&
            y_num = this%leaf_division(2),&
            z_num = this%leaf_division(3) &
            )
        call this%leaf(num_leaf)%femdomain%resize(&
                y = random%gauss(mu=this%leaf_thickness_ave(i),sigma=this%leaf_thickness_sig(i))  , &
                z = random%gauss(mu=this%leaf_length_ave(i)   ,sigma=this%leaf_length_sig(i)) , &
                x = random%gauss(mu=this%leaf_width_ave(i)    ,sigma=this%leaf_width_sig(i)) &
            )
        call this%leaf(num_leaf)%curve(curvature=this%leaf_curvature(i) )
        
        call this%leaf(num_leaf)%femdomain%rotate(&
                x = radian(random%gauss(mu=this%leaf_angle_ave_x(i),sigma=this%leaf_angle_sig_x(i))), &
                y = 0.0d0, &
                z = radian(random%gauss(mu=this%leaf_angle_ave_z(i),sigma=this%leaf_angle_sig_z(i)) ) &
            )
        call this%leaf(num_leaf)%connect("=>",this%stem( this%Leaf_From(i) ))
            this%leaf2stem(num_leaf, this%Leaf_From(i) ) = 1
    enddo

    ! set panicle

    ! set panicles
    ! get number of panicles
    this%num_panicle=1
    do
        if(this%num_panicle ==  this%mainstem_node)then
            this%num_panicle = this%num_panicle -1
            exit
        endif

        line = Riceconfig%parse(config,key1="Panicle_"//str(this%num_panicle)//"_",key2="From")
        
        if(len(trim(line))==0)then
            this%num_panicle = this%num_panicle -1
            exit
        else
            this%num_panicle = this%num_panicle  + 1
            cycle
        endif
        
    enddo

    allocate(this%panicle(this%num_panicle) )
    this%panicle2stem = zeros(this%num_panicle, size(this%stem2stem,1 ) )
    do i=1,this%num_panicle
        call this%panicle(i)%init(&
            x_num = this%panicle_division(1),&
            y_num = this%panicle_division(2),&
            z_num = this%panicle_division(3),&
            rice=.true.,&
            Length=freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="Length")),&
            Width= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="Width")),&
            rice_seed_interval= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_seed_interval")),&
            rice_seed_branch_length= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_seed_branch_length")),&
            rice_seed_length= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_seed_length")),&
            rice_seed_width= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_seed_width")),&
            rice_seed_thickness= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_seed_thickness")),&
            rice_panicle_curvature= freal(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="rice_panicle_curvature")) &
            )
        n = fint(Riceconfig%parse(config,key1="Panicle_"//str(i)//"_",key2="From"))
        this%panicle2stem(i,n) = 1
        
        ! こいつを実装する．
        call this%panicle(i)%connect("=>",this%stem(n) )

    enddo


    
    call this%update()
    


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

! #############################################################
recursive subroutine createRiceShoot(this,config,ShootIdx,debug)
    class(Rice_),intent(inout) :: this
    character(*),intent(in) :: config
    integer(int32),intent(in) :: ShootIdx

    character(:),allocatable :: line
    logical,optional,intent(in) :: debug
    logical :: debug_log
    type(IO_) :: Riceconfig
    type(Random_) :: random
    integer(int32)::i,n,j,k,num_leaf,num_stem_node,num_branch_branch,cpid,shoot_idx
    real(real64) :: x_A(1:3)



    debug_log = input(default=.false.,option=debug)
    cpid = 0

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

!    this%mainstem_length = freal(Riceconfig%parse(config,key1="Mainstem",key2="Length"))
!    this%mainstem_width = freal(Riceconfig%parse(config,key1="Mainstem",key2="Width"))
!    this%mainstem_node = fint(Riceconfig%parse(config,key1="Mainstem",key2="Node"))
!    this%ms_angle_ave = freal(Riceconfig%parse(config,key1="Mainstem",key2="ms_angle_ave"))
!    this%ms_angle_sig = freal(Riceconfig%parse(config,key1="Mainstem",key2="ms_angle_sig"))
!
!    print *, this%mainstem_length &
!    ,this%mainstem_width &
!    ,this%mainstem_node &
!    ,this%ms_angle_ave &
!    ,this%ms_angle_sig

    this%mainstem_length = freal(Riceconfig%parse_json(&
        config,to_list("Shoot_"+str(ShootIdx)+"_","Mainstem","Length")))
    this%mainstem_width = freal(Riceconfig%parse_json(&
        config,to_list("Shoot_"+str(ShootIdx)+"_","Mainstem","Width")))
    this%mainstem_node = fint(Riceconfig%parse_json(&
        config,to_list("Shoot_"+str(ShootIdx)+"_","Mainstem","Node")))
    this%ms_angle_ave = freal(Riceconfig%parse_json(&
        config,to_list("Shoot_"+str(ShootIdx)+"_","Mainstem","ms_angle_ave")))
    this%ms_angle_sig = freal(Riceconfig%parse_json(&
        config,to_list("Shoot_"+str(ShootIdx)+"_","Mainstem","ms_angle_sig")))
    
    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
    
    ! get number of leaf
    this%num_leaf=0
    do
        if(this%num_leaf ==  this%mainstem_node)then
            exit
        endif

        ! line = Riceconfig%parse(config,key1="Leaf_"//str(this%num_leaf)//"_",key2="From")
        line = Riceconfig%parse_json(config,to_list(&
            "Shoot_"+str(ShootIdx)+"_","Leaf_"//str(this%num_leaf+1)//"_","From"))
        
        if ("not found" .in. line)then
            exit
        else
            this%num_leaf = this%num_leaf  + 1
        endif

        !if(len(trim(line))==0)then
        !    this%num_leaf = this%num_leaf -1
        !    exit
        !else
        !    this%num_leaf = this%num_leaf  + 1
        !    cycle
        !endif
        
    enddo


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    allocate(this%leaf_curvature(this%num_leaf))
    allocate(this%leaf_thickness_ave(this%num_leaf))
    allocate(this%leaf_thickness_sig(this%num_leaf))
    allocate(this%leaf_angle_ave_x(this%num_leaf))
    allocate(this%leaf_angle_sig_x(this%num_leaf))
    allocate(this%leaf_angle_ave_z(this%num_leaf))
    allocate(this%leaf_angle_sig_z(this%num_leaf))
    allocate(this%leaf_length_ave(this%num_leaf))
    allocate(this%leaf_length_sig(this%num_leaf))
    allocate(this%leaf_width_ave(this%num_leaf))
    allocate(this%leaf_width_sig(this%num_leaf))
    allocate(this%leaf_From(this%num_leaf))
    !allocate(this%leaf_Length(this%num_leaf))
    !allocate(this%leaf_Width(this%num_leaf))


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    do i=1,this%num_leaf
        !this%leaf_From(i)= fint(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="From"))
        !this%leaf_curvature(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_curvature"))
        !this%leaf_thickness_ave(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_thickness_ave"))
        !this%leaf_thickness_sig(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_thickness_sig"))
        !this%leaf_angle_ave_x(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_ave_x"))
        !this%leaf_angle_sig_x(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_sig_x"))
        !this%leaf_angle_ave_z(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_ave_z"))
        !this%leaf_angle_sig_z(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_angle_sig_z"))
        !this%leaf_length_ave(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_length_ave"))
        !this%leaf_length_sig(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_length_sig"))
        !this%leaf_width_ave(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_width_ave"))
        !this%leaf_width_sig(i)= freal(Riceconfig%parse(&
        !    config,key1="Leaf_"//str(i)//"_",key2="leaf_width_sig"))
        
        this%leaf_From(i)= fint(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","From")))
        this%leaf_curvature(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_curvature")))
        this%leaf_thickness_ave(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_thickness_ave")))
        this%leaf_thickness_sig(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_thickness_sig")))
        this%leaf_angle_ave_x(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_angle_ave_x")))
        this%leaf_angle_sig_x(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_angle_sig_x")))
        this%leaf_angle_ave_z(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_angle_ave_z")))
        this%leaf_angle_sig_z(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_angle_sig_z")))
        this%leaf_length_ave(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_length_ave")))
        this%leaf_length_sig(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_length_sig")))
        this%leaf_width_ave(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_width_ave")))
        this%leaf_width_sig(i)= freal(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Leaf_"//str(i)//"_","leaf_width_sig")))
        
    enddo
    
    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
!    this%mainroot_length = freal(Riceconfig%parse(config,key1="Mainroot",key2="Length"))
!    this%mainroot_width = freal(Riceconfig%parse(config,key1="Mainroot",key2="Width"))
!    this%mainroot_node = fint(Riceconfig%parse(config,key1="Mainroot",key2="Node"))



    ! get number of branch && number of node
!    this%num_branch_root=1
!    this%num_branch_root_node=0
!    do 
!        line = Riceconfig%parse(config,key1="Branchroot_"//str(this%num_branch_root)//"_",key2="Node" )
!        if(len(trim(line))==0)then
!            this%num_branch_root = this%num_branch_root -1
!            exit
!        else
!            this%num_branch_root = this%num_branch_root  + 1
!            this%num_branch_root_node = this%num_branch_root_node + fint(line)
!            cycle
!        endif
!    enddo


    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif
    this%num_stem = this%mainstem_node 
    !this%num_root =this%num_branch_root_node + this%mainroot_node

    !allocate(this%leaf_list(this%num_leaf))
    !allocate(this%stem_list(this%num_stem))
    !allocate(this%root_list(this%num_root))

    allocate(this%leaf(this%num_leaf))
    allocate(this%stem(this%num_stem))
    !allocate(this%root(this%num_root))



    this%leaf2stem = zeros( this%num_leaf , this%num_stem) 
    this%stem2stem = zeros( this%num_stem, this%num_stem)
    this%panicle2stem = zeros( this%num_panicle, this%num_stem) 
    !this%root2stem = zeros( this%num_root , this%num_stem) 
    !this%root2root = zeros( this%num_root , this%num_root ) 

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    ! set mainstem
    do i=1,this%mainstem_node

        call this%stem(i)%init(&
            x_num = this%stem_division(1),&
            y_num = this%stem_division(2),&
            z_num = this%stem_division(3) &
            )
             
        call this%stem(i)%resize(&
            x = this%mainstem_width, &
            y = this%mainstem_width, &
            z = this%mainstem_length/dble(this%mainstem_node) &
            )
        call this%stem(i)%rotate(&
            x = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig)),  &
            y = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig)),  &
            z = radian(random%gauss(mu=this%ms_angle_ave,sigma=this%ms_angle_sig))   &
            )
    enddo

    do i=1,this%mainstem_node-1
        call this%stem(i+1)%connect("=>",this%stem(i))
        this%stem2stem(i+1,i) = 1
    enddo
    

    if(debug_log)then
        cpid = cpid + 1
        call print("createRice #" + str(cpid) )
    endif

    !set leaf
    num_leaf = 0
    do i=1,this%num_leaf
        ! 1葉/1節
        ! add leaves
        
        num_leaf=num_leaf+1
        
        call this%leaf(num_leaf)%init(species=PF_RICE,&
            x_num = this%leaf_division(1),&
            y_num = this%leaf_division(2),&
            z_num = this%leaf_division(3) &
            )
        call this%leaf(num_leaf)%femdomain%resize(&
                y = random%gauss(mu=this%leaf_thickness_ave(i),sigma=this%leaf_thickness_sig(i))  , &
                z = random%gauss(mu=this%leaf_length_ave(i)   ,sigma=this%leaf_length_sig(i)) , &
                x = random%gauss(mu=this%leaf_width_ave(i)    ,sigma=this%leaf_width_sig(i)) &
            )
        call this%leaf(num_leaf)%curve(curvature=this%leaf_curvature(i) )
        
        call this%leaf(num_leaf)%femdomain%rotate(&
                x = radian(random%gauss(mu=this%leaf_angle_ave_x(i),sigma=this%leaf_angle_sig_x(i))), &
                y = 0.0d0, &
                z = radian(random%gauss(mu=this%leaf_angle_ave_z(i),sigma=this%leaf_angle_sig_z(i)) ) &
            )
        call this%leaf(num_leaf)%connect("=>",this%stem( this%Leaf_From(i) ))
            this%leaf2stem(num_leaf, this%Leaf_From(i) ) = 1
    enddo

    ! set panicle

    ! set panicles
    ! get number of panicles
    this%num_panicle=0
    do
        if(this%num_panicle ==  this%mainstem_node)then
            this%num_panicle = this%num_panicle -1
            exit
        endif

        line = Riceconfig%parse_json(config,to_list(&
            "Shoot_"+str(ShootIdx)+"_","Panicle_"//str(this%num_panicle+1)//"_","From"))
        
        
        if ("not found" .in. line)then
            exit
        else
            this%num_panicle = this%num_panicle  + 1
        endif
        
    enddo

    allocate(this%panicle(this%num_panicle) )
    this%panicle2stem = zeros(this%num_panicle, size(this%stem2stem,1 ) )
    do i=1,this%num_panicle
        call this%panicle(i)%init(&
            x_num = this%panicle_division(1),&
            y_num = this%panicle_division(2),&
            z_num = this%panicle_division(3),&
            rice=.true.,&
            Length=freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","Length"))),&
            Width= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","Width"))),&
            rice_seed_interval= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_seed_interval"))),&
            rice_seed_branch_length= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_seed_branch_length"))),&
            rice_seed_length= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_seed_length"))),&
            rice_seed_width= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_seed_width"))),&
            rice_seed_thickness= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_seed_thickness"))),&
            rice_panicle_curvature= freal(Riceconfig%parse_json(&
                config,to_list("Shoot_"+str(ShootIdx)+"_",&
                    "Panicle_"//str(i)//"_","rice_panicle_curvature"))) &
            )
        n = fint(Riceconfig%parse_json(&
            config,to_list("Shoot_"+str(ShootIdx)+"_","Panicle_"//str(i)//"_","From")))
        this%panicle2stem(i,n) = 1
        
        ! こいつを実装する．
        call this%panicle(i)%connect("=>",this%stem(n) )

    enddo


    
    call this%update()
    


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


! ########################################
subroutine mshRice(this,name,num_threads)
    class(Rice_),intent(inout) :: this
    character(*),intent(in) :: name
    integer(int32),optional,intent(in) :: num_threads
    integer(int32) :: i,n

    n = input(default=1,option=num_threads)
    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%stem)
        if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%stem(i)%msh(name=name//"_stem"//str(i))
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%root)
        if(this%root(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%root(i)%msh(name=name//"_root"//str(i))
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%leaf)
        if(this%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%leaf(i)%msh(name=name//"_leaf"//str(i))
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

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


! ########################################
subroutine vtkRice(this,name,num_threads,single_file,&
    scalar_field,vector_field,tensor_field,field_name)
    class(Rice_),intent(inout) :: this
    character(*),intent(in) :: name
    integer(int32),optional,intent(in) :: num_threads
    logical,optional,intent(in) :: single_file
    character(*),optional,intent(in) :: field_name
    integer(int32) :: i,n
    type(FEMDomain_) :: femdomain
    real(real64),optional,intent(in) :: scalar_field(:)
    real(real64),optional,intent(in) :: vector_field(:,:)
    real(real64),optional,intent(in) :: tensor_field(:,:,:)

    if(present(single_file) )then
        if(single_file)then
            ! export mesh for a single file
            if(allocated(this%stem) )then
                do i=1,size(this%stem)
                    if(.not.this%stem(i)%femdomain%empty() )then
                        femdomain = femdomain + this%stem(i)%femdomain
                    endif
                enddo
            endif

            if(allocated(this%leaf) )then
                do i=1,size(this%leaf)
                    if(.not.this%leaf(i)%femdomain%empty() )then
                        femdomain = femdomain + this%leaf(i)%femdomain
                    endif
                enddo
            endif

            if(allocated(this%Panicle) )then
                do i=1,size(this%Panicle)
                    if(.not.this%Panicle(i)%femdomain%empty() )then
                        femdomain = femdomain + this%Panicle(i)%femdomain
                    endif
                enddo
            endif

            if(allocated(this%root) )then
                do i=1,size(this%root)
                    if(.not.this%root(i)%femdomain%empty() )then
                        femdomain = femdomain + this%root(i)%femdomain
                    endif
                enddo
            endif

            if(present(scalar_field) )then
                ! export scalar-valued field 
                ! as a single file
                call femdomain%vtk(field=field_name,name=name,scalar=scalar_field)
            elseif(present(vector_field) )then
                ! export vector-valued field 
                ! as a single file
                call femdomain%vtk(field=field_name,name=name,vector=vector_field)
            elseif(present(tensor_field) )then
                ! export tensor-valued field 
                ! as a single file
                call femdomain%vtk(field=field_name,name=name,tensor=tensor_field)
            else
                call femdomain%vtk(field=field_name,name=name)
            endif
            return
        endif
    endif


    n = input(default=1,option=num_threads)
    if(allocated(this%stem) )then
        !$OMP parallel num_threads(n) private(i)
        !$OMP do 
        do i=1,size(this%stem)
            if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%stem(i)%vtk(name=name//"_stem"//str(i))
            endif
        enddo
        !$OMP end do
        !$OMP end parallel
    endif

    if(allocated(this%root))then
        !$OMP parallel num_threads(n) private(i)
        !$OMP do 
        do i=1,size(this%root)
            if(this%root(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%root(i)%vtk(name=name//"_root"//str(i))
            endif
        enddo
        !$OMP end do
        !$OMP end parallel
    endif


    if(allocated(this%leaf))then
        !$OMP parallel num_threads(n) private(i)
        !$OMP do 
        do i=1,size(this%leaf)
            if(this%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%leaf(i)%vtk(name=name//"_leaf"//str(i))
            endif
        enddo
        !$OMP end do
        !$OMP end parallel
    endif


    if(allocated(this%panicle))then
        !$OMP parallel num_threads(n) private(i)
        !$OMP do 
        do i=1,size(this%panicle)
            if(this%panicle(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%panicle(i)%vtk(name=name//"_panicle"//str(i))
            endif
        enddo
        !$OMP end do
        !$OMP end parallel
    endif

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


! ########################################
subroutine jsonRice(this,name)
    class(Rice_),intent(inout) :: this
    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(this%stem)
        if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            countnum=countnum+1
            call f%write('"'//"stem"//str(i)//'":')
            call this%stem(i)%femdomain%json(name=name//"_stem"//str(i),fh=f%fh,endl=.false.)
        endif
    enddo
    call f%write('"num_stem":'//str(countnum)//',' )

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

! ########################################
subroutine stlRice(this,name,num_threads,single_file)
    class(Rice_),intent(inout) :: this
    character(*),intent(in) :: name
    integer(int32),optional,intent(in) :: num_threads
    integer(int32) :: i,n
    logical,optional,intent(in) :: single_file
    logical :: save_as_single_file
    n = input(default=1,option=num_threads)
    save_as_single_file = input(default=.false.,option=single_file)

    if(save_as_single_file)then
        call execute_command_line("echo ' ' > "//name//".stl")
    endif
    
    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%stem)
        if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%stem(i)%stl(name=name//"_stem"//str(i))
            if(save_as_single_file)then
                call execute_command_line("cat "//name//"_stem"//str(i)//".stl >> "//name//".stl")
                call execute_command_line("rm "//name//"_stem"//str(i)//".stl")
            endif
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%root)
        if(this%root(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%root(i)%stl(name=name//"_root"//str(i))
            if(save_as_single_file)then
                call execute_command_line("cat "//name//"_root"//str(i)//".stl >> "//name//".stl")
                call execute_command_line("rm "//name//"_root"//str(i)//".stl")
            endif
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%leaf)
        if(this%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%leaf(i)%stl(name=name//"_leaf"//str(i))
            if(save_as_single_file)then
                call execute_command_line("cat "//name//"_leaf"//str(i)//".stl >> "//name//".stl")
                call execute_command_line("rm "//name//"_leaf"//str(i)//".stl")
            endif
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(this%Panicle)
        if(this%Panicle(i)%femdomain%mesh%empty() .eqv. .false. )then
            call this%Panicle(i)%stl(name=name//"_Panicle"//str(i))
            if(save_as_single_file)then
                call execute_command_line("cat "//name//"_Panicle"//str(i)//".stl >> "//name//".stl")
                call execute_command_line("rm "//name//"_Panicle"//str(i)//".stl")
            endif
        endif
    enddo
    !$OMP end do
    !$OMP end parallel


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

! ########################################
subroutine moveRice(this,x,y,z)
    class(Rice_),intent(inout) :: this
    real(real64),optional,intent(in) :: x,y,z
    integer(int32) :: i
    if(allocated(this%stem) )then
        do i=1,size(this%stem)
            if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%stem(i)%move(x=x,y=y,z=z)
            endif
        enddo
    endif
    if(allocated(this%root) )then
        do i=1,size(this%root)
            if(this%root(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%root(i)%move(x=x,y=y,z=z)
            endif
        enddo
    endif
    if(allocated(this%leaf) )then
        do i=1,size(this%leaf)
            if(this%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%leaf(i)%move(x=x,y=y,z=z)
            endif
        enddo
    endif
    
    if(allocated(this%panicle) )then
        do i=1,size(this%panicle)
            if(this%panicle(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%panicle(i)%move(x=x,y=y,z=z)
            endif
        enddo
    endif

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



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

    if(allocated(this%stem) )then
        do i=1,size(this%stem)
            if(this%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%stem(i)%rotate(x=x,y=y,z=z)
            endif
        enddo
    endif


    if(allocated(this%leaf) )then
        do i=1,size(this%leaf)
            if(this%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%leaf(i)%rotate(x=x,y=y,z=z)
            endif
        enddo
    endif

    if(allocated(this%root) )then
        do i=1,size(this%root)
            if(this%root(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%root(i)%rotate(x=x,y=y,z=z)
            endif
        enddo
    endif


    if(allocated(this%panicle) )then
        do i=1,size(this%panicle)
            if(this%panicle(i)%femdomain%mesh%empty() .eqv. .false. )then
                call this%panicle(i)%rotate(x=x,y=y,z=z)
            endif
        enddo
    endif

    call this%update()

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



! ########################################
recursive subroutine updateRice(this,stem_id, root_id, leaf_id, overset_margin,debug)
    class(Rice_),intent(inout) :: this
    integer(int32),optional,intent(in) :: stem_id, root_id, leaf_id    
    real(real64),optional,intent(in) :: overset_margin
    integer(int32) :: i,j,this_stem_id,next_stem_id,A_id,B_id,itr_tol,itr,k,kk
    integer(int32) :: this_leaf_id,next_leaf_id
    integer(int32) :: this_root_id,next_root_id,InterNodeID,PetioleID,StemID,LeafID
    real(real64) :: x_A(3),x_B(3),diff(3),error,last_error,mgn,overset_m,error_tol
    logical,optional,intent(in) :: debug
    integer(int32) :: next_panicle_id

!    if(this%default_Leaf_growth_ratio > 0.0d0)then
!        do i=1,size(this%leaf)
!            if(this%leaf(i)%empty() ) cycle
!            this%leaf(i)%length_growth_ratio = this%default_Leaf_growth_ratio
!            this%leaf(i)%Width_growth_ratio = this%default_Leaf_growth_ratio
!        enddo
!    endif
!
!
!    if(this%default_stem_growth_ratio > 0.0d0)then
!        do i=1,size(this%stem)
!            if(this%stem(i)%empty() ) cycle
!            this%stem(i)%length_growth_ratio = this%default_stem_growth_ratio
!            this%stem(i)%Width_growth_ratio = this%default_stem_growth_ratio
!        enddo
!    endif

    ! if Rice_internode_info_ is active
    ! update parameters
!    if(allocated(this%InterNodeInfo) )then
!        do i=0,this%MaxBranchNum
!            
!            if(allocated(this%InterNodeInfo(i)%FinalInterNodeLength ) )then
!                do j=1,this%maxInterNodeID(StemID=i)
!                    InterNodeID = this%searchStem(StemID=i,InterNodeID=j)
!                    if(size(this%InterNodeInfo(i)%FinalInterNodeLength) < j ) then
!                        print *, "ERROR :: updateRice >> "
!                        print *, "size(this%InterNodeInfo(i)%FinalInterNodeLength) is not enough"
!                        stop
!                    endif
!                    if(InterNodeID<1)then
!                        cycle
!                    endif
!                    this%stem(InterNodeID)%final_length = this%InterNodeInfo(i)%FinalInterNodeLength(j)
!                enddo
!            endif
!            
!            if(allocated(this%InterNodeInfo(i)%FinalPetioleLength) )then
!                do j=1,this%maxInterNodeID(StemID=i)
!                    do k=1,this%maxPetioleID(StemID=i,InterNodeID=j)
!                        if(size(this%InterNodeInfo(i)%FinalPetioleLength) < j ) then
!                            print *, "ERROR :: updateRice >> "
!                            print *, "size(this%InterNodeInfo(i)%FinalInterNodeLength) is not enough"
!                            stop
!                        endif
!
!                        PetioleID = this%searchPetiole(StemID=i,InterNodeID=j,PetioleID=k)
!                        
!                        this%stem(PetioleID)%final_length = this%InterNodeInfo(i)%FinalPetioleLength(j)
!                    enddo
!                enddo
!            endif
!
!                if(allocated(this%InterNodeInfo(i)%FinalLeafLength) )then
!                    do j=1,this%maxInterNodeID(StemID=i)
!                        do k=1,this%maxPetioleID(StemID=i,InterNodeID=j)
!                            do kk = 1, this%maxleafID(StemID=i,InterNodeID=j,PetioleID=k)
!                                if(size(this%InterNodeInfo(i)%FinalLeafLength) < j ) then
!                                    print *, "ERROR :: updateRice >> "
!                                    print *, "size(this%InterNodeInfo(i)%FinalInterNodeLength) is not enough"
!                                    stop
!                                endif
!                                LeafID = this%searchleaf(StemID=i,InterNodeID=j,PetioleID=k,LeafID=kk)
!                                this%leaf(LeafID)%final_length = this%InterNodeInfo(i)%FinalLeafLength(j)
!                            enddo
!                        enddo
!                    enddo
!                endif
!
!
!                if(allocated(this%InterNodeInfo(i)%FinalLeafWidth) )then
!                    do j=1,this%maxInterNodeID(StemID=i)
!                        do k=1,this%maxPetioleID(StemID=i,InterNodeID=j)
!                            do kk = 1, this%maxleafID(StemID=i,InterNodeID=j,PetioleID=k)
!                                if(size(this%InterNodeInfo(i)%FinalLeafWidth) < j ) then
!                                    print *, "ERROR :: updateRice >> "
!                                    print *, "size(this%InterNodeInfo(i)%FinalInterNodeLength) is not enough"
!                                    stop
!                                endif
!                                LeafID = this%searchleaf(StemID=i,InterNodeID=j,PetioleID=k,LeafID=kk)
!                                this%leaf(LeafID)%final_Width = this%InterNodeInfo(i)%FinalLeafWidth(j)
!                            enddo
!                        enddo
!                    enddo
!                endif
!            
!        enddo
!    endif

    ! update connectivity
    if(.not. allocated(this%stem2stem ))then
        print *, "updateRice >> ERROR :: .not. allocated(this%stem2stem )"
        return
    endif



    error_tol = dble(1.0e-14)

    ! margin between subdomains
    overset_m = input(default=0.03d0, option=overset_margin)
    
    itr_tol = 100
    itr=0

    ! if debug
    !if(present(debug) )then
    !    if(debug)then
    !        print *, "this%stem2stem"
    !        call print(this%stem2stem)
    !    endif
    !endif

    ! stem to stem
    last_error = 1.0d0
    if(maxval(this%stem2stem)/=0) then
        

        do 
            itr=itr+1
            error = 0.0d0
            do i=1, size(this%stem2stem,1)
                do j=1, size(this%stem2stem,2)
                    this_stem_id = j
                    next_stem_id = i
                    if(this%stem2stem(i,j)/=0 .and. i /= j)then
                        ! this_stem_id ===>>> next_stem_id, connected!

                        !x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
                        !x_A(:) = this%stem(next_stem_id)%getCoordinate("A")
                        ! Overset分食い込ませる
                        x_B(:) = (1.0d0-overset_m)*this%stem(this_stem_id)%getCoordinate("B")&
                            + overset_m*this%stem(this_stem_id)%getCoordinate("A")
                        ! Overset分食い込ませる
                        x_A(:) = (1.0d0-overset_m)*this%stem(next_stem_id)%getCoordinate("A") &
                            + overset_m*this%stem(next_stem_id)%getCoordinate("B")


                        diff(:) = x_B(:) - x_A(:)
                        
                        error = error + dot_product(diff,diff)
                        call this%stem(next_stem_id)%move(x=diff(1),y=diff(2),z=diff(3) )

                    endif
                enddo
            enddo
            if(present(debug) )then
                if(debug)then
                    print *, "Rice % update s2s >> error :: ",error
                endif
            endif
            if(itr > itr_tol) then
                print *, "Rice % update s2s >> ERROR :: not converged"
                stop
            endif

            if( abs(error) + abs(last_error) < error_tol) exit
            last_error = error
        enddo
    endif

    ! root to stem
    if(allocated(this%root2stem) )then
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(this%root2stem,1)
            do j=1, size(this%root2stem,2)
                this_stem_id = j
                next_root_id = i
                if(this%root2stem(i,j)==1)then
                    ! this_stem_id ===>>> next_root_id, connected!
                    !x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
                    !x_A(:) = this%root(next_root_id)%getCoordinate("A")

                    ! Overset分食い込ませる
                    x_B(:) = (1.0d0-overset_m)*this%stem(this_stem_id)%getCoordinate("A")&
                        + overset_m*this%stem(this_stem_id)%getCoordinate("B")
                    ! Overset分食い込ませる
                    x_A(:) = (1.0d0-overset_m)*this%root(next_root_id)%getCoordinate("A") &
                        + overset_m*this%root(next_root_id)%getCoordinate("B")
                        

                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call this%root(next_root_id)%move(x=diff(1),y=diff(2),z=diff(3) )
                endif
            enddo
        enddo
        if(present(debug) )then
            if(debug)then
                print *, "Rice % update r2s >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "Rice % update r2s  >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error)  < error_tol) exit
        last_error = error
    enddo
    endif

    
    if(allocated(this%root2root) )then
    ! root to root
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(this%root2root,1)
            do j=1, size(this%root2root,2)
                this_root_id = j
                next_root_id = i
                if(next_root_id==1)then
                    cycle
                endif
                if(this%root2root(i,j)/=0 .and. i /= j)then
                    ! this_root_id ===>>> next_root_id, connected!
                    !x_B(:) = this%root(this_root_id)%getCoordinate("B")
                    !x_A(:) = this%root(next_root_id)%getCoordinate("A")
                    
                    ! Overset分食い込ませる
                    x_B(:) = (1.0d0-overset_m)*this%root(this_root_id)%getCoordinate("B")&
                        + overset_m*this%root(this_root_id)%getCoordinate("A")
                    ! Overset分食い込ませる
                    x_A(:) = (1.0d0-overset_m)*this%root(next_root_id)%getCoordinate("A") &
                        + overset_m*this%root(next_root_id)%getCoordinate("B")


                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call this%root(next_root_id)%move(x=diff(1),y=diff(2),z=diff(3) )
                endif
            enddo
        enddo
        if(present(debug) )then
            if(debug)then
                print *, "Rice % update r2r >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "Rice % update r2r >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error)  < error_tol) exit
        last_error = error
    enddo
    endif

    ! leaf to stem
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(this%leaf2stem,1)
            do j=1, size(this%leaf2stem,2)
                this_stem_id = j
                next_leaf_id = i
                if(this%leaf2stem(i,j)==1)then
                    ! this_stem_id ===>>> next_leaf_id, connected!
                    !x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
                    !x_A(:) = this%leaf(next_leaf_id)%getCoordinate("A")

                    ! Overset分食い込ませる
                    x_B(:) = (1.0d0-overset_m)*this%stem(this_stem_id)%getCoordinate("B")&
                        + overset_m*this%stem(this_stem_id)%getCoordinate("A")
                    ! Overset分食い込ませる
                    x_A(:) = this%leaf(next_leaf_id)%getCoordinate("A")
                        

                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call this%leaf(next_leaf_id)%move(x=diff(1),y=diff(2),z=diff(3) )
                endif
            enddo
        enddo
        if(present(debug) )then
            if(debug)then
                print *, "Rice % update l2s >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "Rice % update l2s  >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error)  < error_tol) exit
        last_error = error
    enddo

    

    
    ! panicle to stem
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(this%panicle2stem,1)
            do j=1, size(this%panicle2stem,2)
                this_stem_id = j
                next_panicle_id = i
                if(this%panicle2stem(i,j)==1)then
                    ! this_stem_id ===>>> next_panicle_id, connected!
                    !x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
                    !x_A(:) = this%panicle(next_panicle_id)%getCoordinate("A")

                    ! Overset分食い込ませる
                    x_B(:) = (1.0d0-overset_m)*this%stem(this_stem_id)%getCoordinate("B")&
                        + overset_m*this%stem(this_stem_id)%getCoordinate("A")
                    ! Overset分食い込ませる
                    x_A(:) = this%panicle(next_panicle_id)%getCoordinate("A")
                        

                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call this%panicle(next_panicle_id)%move(x=diff(1),y=diff(2),z=diff(3) )
                endif
            enddo
        enddo
        
        if(present(debug) )then
            if(debug)then
                print *, "Rice % update l2s >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "Rice % update l2s  >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error)  < error_tol) exit
        last_error = error
    enddo

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



! ########################################
!recursive subroutine updateRice(this,debug)
!    class(Rice_),intent(inout) :: this
!    integer(int32) :: i,j,this_stem_id,next_stem_id,A_id,B_id,itr_tol,itr
!    integer(int32) :: this_leaf_id,next_leaf_id
!    integer(int32) :: this_root_id,next_root_id
!    real(real64) :: x_A(3),x_B(3),diff(3),error,last_error
!    logical,optional,intent(in) :: debug
!    
!    
!    ! update connectivity
!    if(.not. allocated(this%stem2stem ))then
!        print *, "updateRice >> ERROR :: .not. allocated(this%stem2stem )"
!        return
!    endif
!
!    itr_tol = 100
!    itr=0
!
!    ! stem to stem
!    last_error = 1.0d0
!    do 
!        itr=itr+1
!        error = 0.0d0
!        do i=1, size(this%stem2stem,1)
!            do j=1, size(this%stem2stem,2)
!                this_stem_id = j
!                next_stem_id = i
!                if(this%stem2stem(i,j)/=0 .and. i /= j)then
!                    ! this_stem_id ===>>> next_stem_id, connected!
!                    x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
!                    x_A(:) = this%stem(next_stem_id)%getCoordinate("A")
!                    diff(:) = x_B(:) - x_A(:)
!                    error = error + dot_product(diff,diff)
!                    call this%stem(next_stem_id)%move(x=diff(1),y=diff(2),z=diff(3) )
!                endif
!            enddo
!        enddo
!        if(present(debug) )then
!            if(debug)then
!                print *, "Rice % update >> error :: ",error
!            endif
!        endif
!        if(itr > itr_tol) then
!            print *, "Rice % update >> ERROR :: not converged"
!            stop
!        endif
!        
!        if( abs(error) + abs(last_error) == 0.0d0) exit
!        last_error = error
!    enddo
!
!    ! leaf to stem
!    last_error = 1.0d0
!    do 
!        itr=itr+1
!        error = 0.0d0
!        do i=1, size(this%leaf2stem,1)
!            do j=1, size(this%leaf2stem,2)
!                this_stem_id = j
!                next_leaf_id = i
!                if(this%leaf2stem(i,j)==1)then
!                    ! this_stem_id ===>>> next_leaf_id, connected!
!                    x_B(:) = this%stem(this_stem_id)%getCoordinate("B")
!                    x_A(:) = this%leaf(next_leaf_id)%getCoordinate("A")
!                    diff(:) = x_B(:) - x_A(:)
!                    error = error + dot_product(diff,diff)
!                    call this%leaf(next_leaf_id)%move(x=diff(1),y=diff(2),z=diff(3) )
!                endif
!            enddo
!        enddo
!        if(present(debug) )then
!            if(debug)then
!                print *, "Rice % update >> error :: ",error
!            endif
!        endif
!        
!        if(itr > itr_tol) then
!            print *, "Rice % update >> ERROR :: not converged"
!            stop
!        endif
!        
!        if( abs(error) + abs(last_error) == 0.0d0) exit
!        last_error = error
!    enddo
!
!
!    return
!
!
!    ! root to root
!    last_error = 1.0d0
!    do 
!        itr=itr+1
!        error = 0.0d0
!        do i=1, size(this%root2root,1)
!            do j=1, size(this%root2root,2)
!                this_root_id = j
!                next_root_id = i
!                if(this%root2root(i,j)/=0 .and. i /= j)then
!                    ! this_root_id ===>>> next_root_id, connected!
!                    x_B(:) = this%root(this_root_id)%getCoordinate("B")
!                    x_A(:) = this%root(next_root_id)%getCoordinate("A")
!                    diff(:) = x_B(:) - x_A(:)
!                    error = error + dot_product(diff,diff)
!                    call this%root(next_root_id)%move(x=diff(1),y=diff(2),z=diff(3) )
!                endif
!            enddo
!        enddo
!        if(present(debug) )then
!            if(debug)then
!                print *, "Rice % update >> error :: ",error
!            endif
!        endif
!        if(itr > itr_tol) then
!            print *, "Rice % update >> ERROR :: not converged"
!            stop
!        endif
!        
!        if( abs(error) + abs(last_error) == 0.0d0) exit
!        last_error = error
!    enddo
!
!    
!    
!end subroutine
! ########################################


subroutine removeRice(this,root)
    class(Rice_),intent(inout) :: this
    logical,optional,intent(in) :: root

    if(present(root) )then
        if(root)then
            
            if (allocated(this%Root) ) deallocate(this%Root)
            !if (allocated(this%rootYoungModulus) ) deallocate(this%rootYoungModulus)
            !if (allocated(this%rootPoissonRatio) ) deallocate(this%rootPoissonRatio)
            !if (allocated(this%rootDensity) ) deallocate(this%rootDensity)


            if (allocated(this%root2stem) ) deallocate(this%root2stem)
            if (allocated(this%root2root) ) deallocate(this%root2root)
            !if (allocated(this%root_list) ) deallocate(this%root_list)

            !if (allocated(this%root_angle) ) deallocate(this%root_angle)
            !this%rootconfig=" "
            !this%Num_Of_Root = 0


        endif
        return
    endif

    
    
!    this%leaf_angle_ave(:)=0.0d0
!    this%leaf_angle_sig(:)=0.0d0
!    this%leaf_length_ave(:)=0.0d0
!    this%leaf_length_sig(:)=0.0d0
!    this%leaf_width_ave(:)=0.0d0
!    this%leaf_width_sig(:)=0.0d0
!    this%leaf_thickness_ave(:)=0.0d0
!    this%leaf_thickness_sig(:)=0.0d0
!    
!    this%Stage="" ! VE, CV, V1,V2, ..., R1, R2, ..., R8
!    this%name=""
!    this%stage_id=0
!    this%dt=0.0d0
!    call this%Seed%remove()
!    if (allocated(this%NodeSystem) ) deallocate(this%NodeSystem)
!    if (allocated(this%RootSystem) ) deallocate(this%RootSystem)

    if (allocated(this%Stem) ) deallocate(this%Stem)
    if (allocated(this%Leaf) ) deallocate(this%Leaf)
    if (allocated(this%Root) ) deallocate(this%Root)
    if (allocated(this%Panicle) ) deallocate(this%Panicle)
    

    ! material info
!    if (allocated(this%stemYoungModulus) ) deallocate(this%stemYoungModulus)
!    if (allocated(this%leafYoungModulus) ) deallocate(this%leafYoungModulus)
!    if (allocated(this%rootYoungModulus) ) deallocate(this%rootYoungModulus)
!
!    if (allocated(this%stemPoissonRatio) ) deallocate(this%stemPoissonRatio)
!    if (allocated(this%leafPoissonRatio) ) deallocate(this%leafPoissonRatio)
!    if (allocated(this%rootPoissonRatio) ) deallocate(this%rootPoissonRatio)
!
!    if (allocated(this%stemDensity) ) deallocate(this%stemDensity)
!    if (allocated(this%leafDensity) ) deallocate(this%leafDensity)
!    if (allocated(this%rootDensity) ) deallocate(this%rootDensity)
    

!    if(allocated(this%NodeID_MainStem)) deallocate(this%NodeID_MainStem)
!    if(allocated(this%NodeID_Branch)) deallocate(this%NodeID_Branch)
    ! 節-節点データ構造
    ! call this%struct%remove(all=.true.)
    if (allocated(this%leaf2stem) ) deallocate(this%leaf2stem)
    if (allocated(this%stem2stem) ) deallocate(this%stem2stem)
    if (allocated(this%root2stem) ) deallocate(this%root2stem)
    if (allocated(this%Panicle2stem) ) deallocate(this%Panicle2stem)
    if (allocated(this%root2root) ) deallocate(this%root2root)
    
    ! シミュレータ
    call this%contact%remove()
!    this%time=0.0d0
!    this%seed_length=0.0d0
!    this%seed_width=0.0d0
!    this%seed_height=0.0d0
!    if (allocated(this%stem_angle) ) deallocate(this%stem_angle)
!    if (allocated(this%root_angle) ) deallocate(this%root_angle)
!    if (allocated(this%leaf_angle) ) deallocate(this%leaf_angle)
!
!    this%stemconfig=" "
!    this%rootconfig=" "
!    this%leafconfig=" "

    this%TYPE_STEM    = 1
    this%TYPE_LEAF    = 2
    this%TYPE_ROOT    = 3
    this%TYPE_PANICLE = 4
    ! 節-節点データ構造
    ! call this % struct % remove()
    if(allocated(this%leaf2stem)) deallocate(this%leaf2stem)! (:,:)
    if(allocated(this%stem2stem)) deallocate(this%stem2stem)! (:,:)
    if(allocated(this%panicle2stem)) deallocate(this%panicle2stem)! (:,:)
    if(allocated(this%root2stem)) deallocate(this%root2stem)! (:,:)
    if(allocated(this%root2root)) deallocate(this%root2root)! (:,:)

    this%mainstem_length = 0.0d0
    this%mainstem_width = 0.0d0
    this%mainstem_node=0

    this%mainroot_length = 0.0d0
    this%mainroot_width = 0.0d0
    this%mainroot_node=0


    this%num_branch_root=0
    this%num_branch_root_node=0

    this%ms_angle_ave = 0.0d0
    this%ms_angle_sig = 0.0d0

    if(allocated(this%Leaf_From)) deallocate(this%Leaf_From)! (:)

    !real(real64),allocatable :: leaf_Length(:)
    !real(real64),allocatable :: leaf_Width(:)

    if(allocated(this%leaf_curvature) )deallocate(this%leaf_curvature)! (:)

    if(allocated(this%leaf_thickness_ave) )deallocate(this%leaf_thickness_ave)! (:)
    if(allocated(this%leaf_thickness_sig) )deallocate(this%leaf_thickness_sig)! (:)

    if(allocated(this%leaf_angle_ave_x) )deallocate(this%leaf_angle_ave_x)! (:)
    if(allocated(this%leaf_angle_sig_x) )deallocate(this%leaf_angle_sig_x)! (:)
    if(allocated(this%leaf_angle_ave_z) )deallocate(this%leaf_angle_ave_z)! (:)
    if(allocated(this%leaf_angle_sig_z) )deallocate(this%leaf_angle_sig_z)! (:)

    if(allocated(this%leaf_length_ave) )deallocate(this%leaf_length_ave)! (:)
    if(allocated(this%leaf_length_sig) )deallocate(this%leaf_length_sig)! (:)
    if(allocated(this%leaf_width_ave) )deallocate(this%leaf_width_ave)! (:)
    if(allocated(this%leaf_width_sig) )deallocate(this%leaf_width_sig)! (:)
    
    this%num_leaf=0
    this%num_stem=0
    this%num_panicle = 1
    this%num_root=0


    ! 器官オブジェクト配列
    !if(allocated(this%leaf_list)) deallocate(this%leaf_list)! (:)
    !if(allocated(this%stem_list)) deallocate(this%stem_list)! (:)
    !if(allocated(this%root_list)) deallocate(this%root_list)! (:)

    !this%LeafSurfaceData = ""
    if(allocated(this % Leaf)) deallocate(this % Leaf)! (:)
    if(allocated(this % Stem)) deallocate(this % Stem)! (:)
    if(allocated(this % Panicle)) deallocate(this % Panicle)! (:)
    if(allocated(this % Root)) deallocate(this % Root)! (:)



    
    !if(allocated(this%NodeID_MainStem)) deallocate(this%NodeID_MainStem)! (:)
    !if(allocated(this%NodeID_Branch)) deallocate(this%NodeID_Branch)! (:)
    
    this%inLoop = .false.
    this%hours = 0.0d0
    
    ! growth simulation
    this%FullyExpanded_stem_threshold = 0.10d0
    this%MaxBranchNum = 20
    if(allocated(this%InterNodeInfo)) deallocate(this%InterNodeInfo)! (:)
    this%default_Leaf_growth_ratio = 1.0d0/3.0d0
    this%default_Stem_growth_ratio = 1.0d0/3.0d0
    if(allocated(this%MainStem_num_branch)) deallocate(this%MainStem_num_branch)! (:)
    this%apical_dominance_distance = 1.0d0
    

    ! シミュレータ
    call this%contact%remove()
    this%Gravity_acceralation = 9.810d0
    this%PenaltyParameter = 100000.0d0
    this%GaussPointProjection = .false.

    ! setting
    this%stem_division(1:3) = [10, 10, 10]
    this%leaf_division(1:3) = [10, 10, 10]
    this%panicle_division(1:3) = [2, 2, 300]

end subroutine
! ################################################################
subroutine checkMemoryRequirementRice(this)
    class(Rice_),intent(in) :: this
    real(real64) :: re_val
    integer(int64) :: val
    
    print *, "===================================="
    print *, "checking Memory (RAM) Requirement..."
    print *, "------------------------------------"
    print *, "| thisect type                     | Rice"
    print *, "| Number of points                | "+str(this%nn())
    print *, "| Degree of freedom | Deformation | "+str(this%nn()*3)
    print *, "|                   | ModeAnalysis| ",dble(this%nn())*3*dble(this%nn())*3
    print *, "|                   | Diffusion   | "+str(this%nn())
    print *, "|                   | Reaction    | "+str(this%nn())

    print *, "| DRAM requirement  | Deformation | "+str(this%nn()*3*40*30/1000/1000)+" (MB)"
    val = this%nn()*3*30
    val = val * this%nn()*3/1000/1000
    print *, "|                   | ModeAnalysis| ",str(val)," (MB)"
    print *, "|                   | Diffusion   | "+str(this%nn()*1*20*10/1000/1000)+" (MB)"
    print *, "|                   | Reaction    | "+str(this%nn()*1*20*10/1000/1000)+" (MB)"
    print *, "===================================="


end subroutine

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


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

    ! get number of subdomain
    ret = 0
    do i=1,size(this%stem)
        if( .not.this%stem(i)%femdomain%mesh%empty() ) then
            ret = ret + 1
        endif
    enddo
    do i=1,size(this%leaf)
        if( .not.this%leaf(i)%femdomain%mesh%empty() ) then
            ret = ret + 1
        endif
    enddo
    do i=1,size(this%root)
        if( .not.this%root(i)%femdomain%mesh%empty() ) then
            ret = ret + 1
        endif
    enddo


    do i=1,size(this%Panicle)
        if( .not.this%Panicle(i)%femdomain%mesh%empty() ) then
            ret = ret + 1
        endif
    enddo
end function


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

! ##################################################################
pure function nnRice(this) result(ret)
    class(Rice_) ,intent(in) :: this
    integer(int32) :: ret, i

    ! get number of node (point)
    ret = 0
    
    if(allocated(this%stem) )then
        do i=1,size(this%stem)
            if( .not.this%stem(i)%femdomain%mesh%empty() ) then
                ret = ret + this%stem(i)%femdomain%nn()
            endif
        enddo
    endif

    if(allocated(this%leaf) )then
        do i=1,size(this%leaf)
            if( .not.this%leaf(i)%femdomain%mesh%empty() ) then
                ret = ret + this%leaf(i)%femdomain%nn()
            endif
        enddo
    endif

    if(allocated(this%root) )then
        do i=1,size(this%root)
            if( .not.this%root(i)%femdomain%mesh%empty() ) then
                ret = ret + this%root(i)%femdomain%nn()
            endif
        enddo
    endif



    if(allocated(this%Panicle) )then
        do i=1,size(this%Panicle)
            if( .not.this%Panicle(i)%femdomain%mesh%empty() ) then
                ret = ret + this%Panicle(i)%femdomain%nn()
            endif
        enddo
    endif
end function
! ##################################################################



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

    ! get number of element
    ret = 0
    if(allocated(this%stem) )then
        do i=1,size(this%stem)
            if( .not.this%stem(i)%femdomain%mesh%empty() ) then
                ret = ret + this%stem(i)%femdomain%ne()
            endif
        enddo
    endif

    if(allocated(this%leaf) )then
        do i=1,size(this%leaf)
            if( .not.this%leaf(i)%femdomain%mesh%empty() ) then
                ret = ret + this%leaf(i)%femdomain%ne()
            endif
        enddo
    endif

    if(allocated(this%root) )then
        do i=1,size(this%root)
            if( .not.this%root(i)%femdomain%mesh%empty() ) then
                ret = ret + this%root(i)%femdomain%ne()
            endif
        enddo
    endif


    if(allocated(this%Panicle) )then
        do i=1,size(this%Panicle)
            if( .not.this%Panicle(i)%femdomain%mesh%empty() ) then
                ret = ret + this%Panicle(i)%femdomain%ne()
            endif
        enddo
    endif
end function
! ##################################################################


! ################################################################
recursive subroutine setYoungModulusRice(this,YoungModulus,stem,root,leaf,panicle,ElementList)
    class(Rice_),intent(inout) :: this
    logical,optional,intent(in) :: stem, root, leaf, panicle
    
    ! ElementList(Idx, [TYPE, DOMAIN, ELEMENT]) 
    integer(int32),optional,intent(in) :: ElementList(:,:)

    real(real64),intent(in) :: YoungModulus
    integer(int32) :: i, j, n, domain_idx, elem_idx



    n = 0
    if(present(stem) )then
        if(stem)then
            n=n+1
            if(allocated(this%stem) )then
                do i=1,size(this%stem)
                    if(this%stem(i)%femdomain%empty() )then
                        cycle
                    elseif(present(ElementList) )then
                        do j=1, size(ElementList,1)
                            if(ElementList(j,1) == this%TYPE_STEM )then
                                domain_idx = ElementList(j,2)  
                                elem_idx = ElementList(j,3) 
                                this%stem(domain_idx)%YoungModulus(elem_idx) = YoungModulus
                            endif
                        enddo
                    else
                        this%stem(i)%YoungModulus = YoungModulus*eyes(this%stem(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(leaf) )then
        if(leaf)then
            n=n+10
            if(allocated(this%leaf) )then
                do i=1,size(this%leaf)
                    if(this%leaf(i)%femdomain%empty() )then
                        cycle
                    elseif(present(ElementList) )then
                        do j=1, size(ElementList,1)
                            if(ElementList(j,1) == this%TYPE_LEAF )then
                                domain_idx = ElementList(j,2)  
                                elem_idx = ElementList(j,3) 
                                this%LEAF(domain_idx)%YoungModulus(elem_idx) = YoungModulus
                            endif
                        enddo
                    else
                        this%leaf(i)%YoungModulus = YoungModulus*eyes(this%leaf(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(root) )then
        if(root)then
            n=n+100
            if(allocated(this%root) )then
                do i=1,size(this%root)
                    if(this%root(i)%femdomain%empty() )then
                        cycle
                    elseif(present(ElementList) )then
                        do j=1, size(ElementList,1)
                            if(ElementList(j,1) == this%TYPE_ROOT )then
                                domain_idx = ElementList(j,2)  
                                elem_idx = ElementList(j,3) 
                                this%ROOT(domain_idx)%YoungModulus(elem_idx) = YoungModulus
                            endif
                        enddo
                    else
                        this%root(i)%YoungModulus = YoungModulus*eyes(this%root(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif
    if(present(panicle) )then
        if(panicle)then
            n=n+100
            if(allocated(this%panicle) )then
                do i=1,size(this%panicle)
                    if(this%panicle(i)%femdomain%empty() )then
                        cycle
                    elseif(present(ElementList) )then
                        do j=1, size(ElementList,1)
                            if(ElementList(j,1) == this%TYPE_PANICLE )then
                                domain_idx = ElementList(j,2)  
                                elem_idx = ElementList(j,3) 
                                this%PANICLE(domain_idx)%YoungModulus(elem_idx) = YoungModulus
                            endif
                        enddo
                    else
                        this%panicle(i)%YoungModulus = YoungModulus*eyes(this%panicle(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(n==0)then
        call this%setYoungModulus(YoungModulus=YoungModulus,stem=.true.,root=.true.,leaf=.true.&
        ,Panicle=.true.,ElementList=ElementList)
    endif
    
end subroutine
! ################################################################


! ################################################################
recursive subroutine setPoissonRatioRice(this,PoissonRatio,stem,root,leaf,panicle)
    class(Rice_),intent(inout) :: this
    logical,optional,intent(in) :: stem, root, leaf, panicle
    real(real64),intent(in) :: PoissonRatio
    integer(int32) :: i, n

    n = 0
    if(present(stem) )then
        if(stem)then
            n=n+1
            if(allocated(this%stem) )then
                do i=1,size(this%stem)
                    if(this%stem(i)%femdomain%empty() )then
                        cycle
                    else
                        this%stem(i)%PoissonRatio = PoissonRatio*eyes(this%stem(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(leaf) )then
        if(leaf)then
            n=n+10
            if(allocated(this%leaf) )then
                do i=1,size(this%leaf)
                    if(this%leaf(i)%femdomain%empty() )then
                        cycle
                    else
                        this%leaf(i)%PoissonRatio = PoissonRatio*eyes(this%leaf(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(root) )then
        if(root)then
            n=n+100
            if(allocated(this%root) )then
                do i=1,size(this%root)
                    if(this%root(i)%femdomain%empty() )then
                        cycle
                    else
                        this%root(i)%PoissonRatio = PoissonRatio*eyes(this%root(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(panicle) )then
        if(panicle)then
            n=n+100
            if(allocated(this%panicle) )then
                do i=1,size(this%panicle)
                    if(this%panicle(i)%femdomain%empty() )then
                        cycle
                    else
                        this%panicle(i)%PoissonRatio = PoissonRatio*eyes(this%panicle(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(n==0)then
        call this%setPoissonRatio(PoissonRatio=PoissonRatio,stem=.true.,root=.true.,leaf=.true. &
        ,Panicle=.true.)
    endif
    
end subroutine
! ################################################################


! ################################################################
recursive subroutine setDensityRice(this,Density,stem,root,leaf,panicle)
    class(Rice_),intent(inout) :: this
    logical,optional,intent(in) :: stem, root, leaf, panicle
    real(real64),intent(in) :: Density
    integer(int32) :: i, n

    n = 0
    if(present(stem) )then
        if(stem)then
            n=n+1
            if(allocated(this%stem) )then
                do i=1,size(this%stem)
                    if(this%stem(i)%femdomain%empty() )then
                        cycle
                    else
                        this%stem(i)%Density = Density*eyes(this%stem(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(leaf) )then
        if(leaf)then
            n=n+10
            if(allocated(this%leaf) )then
                do i=1,size(this%leaf)
                    if(this%leaf(i)%femdomain%empty() )then
                        cycle
                    else
                        this%leaf(i)%Density = Density*eyes(this%leaf(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(root) )then
        if(root)then
            n=n+100
            if(allocated(this%root) )then
                do i=1,size(this%root)
                    if(this%root(i)%femdomain%empty() )then
                        cycle
                    else
                        this%root(i)%Density = Density*eyes(this%root(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(present(panicle) )then
        if(panicle)then
            n=n+100
            if(allocated(this%panicle) )then
                do i=1,size(this%panicle)
                    if(this%panicle(i)%femdomain%empty() )then
                        cycle
                    else
                        this%panicle(i)%Density = Density*eyes(this%panicle(i)%femdomain%ne())
                    endif
                enddo
            endif
        endif
    endif

    if(n==0)then
        call this%setDensity(Density=Density,stem=.true.,root=.true.,leaf=.true.&
        ,Panicle=.true.)
    endif
    
end subroutine
! ################################################################


! ################################################################
function getEigenModeRice(this, ground_level,penalty,debug,Frequency,EbOM_Algorithm,num_mode) result(EigenVectors)
    class(Rice_),target,intent(inout) :: this
    real(real64),intent(in) :: ground_level
    real(real64),optional,intent(in) :: penalty
    logical,optional,intent(in) :: debug
    real(real64),allocatable,intent(inout) :: Frequency(:)
    character(*),optional,intent(in) :: EbOM_Algorithm

    integer(int32),optional,intent(in) :: num_mode
    integer(int32) :: num_freq
    !integer(int32),optional,intent(in) :: num_mode
    

    type(FEMDomainp_),allocatable :: FEMDomainPointers(:)
    type(FEMSolver_) :: solver
    type(Math_) :: math

    real(real64),allocatable :: EigenVectors(:,:),buf(:,:),buf_vec(:)

    integer(int32) :: stem_id, leaf_id, root_id,DomainID,ElementID,i,n
    integer(int32) :: myStemID, yourStemID, myLeafID,myRootID, yourRootID
    integer(int32),allocatable :: FixBoundary(:)
    integer(int32) :: nn_domains,EbOM_Algorithm_id
    real(real64) :: vec_norm

    integer(int32) ::  myPanicleID
    real(real64),allocatable :: all_frequency(:),All_EigenVectors(:,:)

    num_freq = input(default=10,option=num_mode)
    
    EbOM_Algorithm_id = FEMDomain_Overset_GPP
    if(present(EbOM_Algorithm) )then
        if(EbOM_Algorithm=="P2P")then
            EbOM_Algorithm_id=FEMDomain_Overset_P2P
        elseif(EbOM_Algorithm=="GPP")then
            EbOM_Algorithm_id=FEMDomain_Overset_P2P
        endif
    endif

    ! linear elasticity with infinitesimal strain theory
    n = this%numStem() + this%numLeaf() + this%numRoot() + this%numPanicle()
    
    allocate(FEMDomainPointers(n) )
    ! ORDER 
    ! [STEM] => [LEAF] => [ROOT] => => [PANICLE]
    !(1) >> compute overset
    ! For stems
    if(allocated(this%stem2stem) )then
        if(allocated(this%stem) )then
            do myStemID = 1,size(this%stem2stem,1)
                do yourStemID = 1, size(this%stem2stem,2)
                    if(this%stem2stem(myStemID,yourStemID)>=1 )then
                        ! connected
                        call this%stem(myStemID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = myStemID  ,&
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%stem(myStemID)%femdomain,&
                            DomainID   = myStemID    ,& 
                            MyDomainID = yourStemID  ,&
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif

    if(allocated(this%leaf2stem) )then
        if(allocated(this%leaf) .and. allocated(this%stem) )then
            do myLeafID = 1,size(this%leaf2stem,1)
                do yourStemID = 1, size(this%leaf2stem,2)
                    if(this%leaf2stem(myLeafID,yourStemID)>=1 )then
                        ! connected
                        call this%leaf(myLeafID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() + myLeafID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%leaf(myLeafID)%femdomain,&
                            DomainID   = this%numStem() + myLeafID    ,& 
                            MyDomainID = yourStemID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif
    



    if(allocated(this%root2stem) )then
        if(allocated(this%stem) .and. allocated(this%root) )then
            do myRootID = 1,size(this%root2stem,1)
                do yourStemID = 1, size(this%root2stem,2)
                    if(this%root2stem(myRootID,yourStemID)>=1 )then
                        ! connected
                        call this%root(myRootID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() +this%numLeaf() + myRootID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=  this%root(myRootID)%femdomain,&
                            DomainID   = this%numStem() +this%numLeaf() + myRootID    ,& 
                            MyDomainID =  yourStemID , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif


    if(allocated(this%root2root) )then
        if(allocated(this%root) )then
            do myRootID = 1,size(this%root2root,1)
                do yourrootID = 1, size(this%root2root,2)
                    if(this%root2root(myRootID,yourrootID)>=1 )then
                        ! connected
                        call this%root(myRootID)%femdomain%overset(&
                            FEMDomain=this%root(yourrootID)%femdomain,&
                            DomainID   = this%numroot() +this%numLeaf() + yourrootID    ,& 
                            MyDomainID = this%numroot() +this%numLeaf() + myRootID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    
                        call this%root(yourrootID)%femdomain%overset(&
                            FEMDomain=this%root(myRootID)%femdomain,&
                            DomainID   = this%numroot() +this%numLeaf() + myRootID    ,& 
                            MyDomainID =  this%numroot() +this%numLeaf() + yourrootID , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif

    
    if(allocated(this%Panicle2stem) )then
        if(allocated(this%Panicle) .and. allocated(this%stem) )then
            do myPanicleID = 1,size(this%Panicle2stem,1)
                do yourStemID = 1, size(this%Panicle2stem,2)
                    if(this%Panicle2stem(myPanicleID,yourStemID)>=1 )then
                        ! connected
                        call this%Panicle(myPanicleID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() + this%numLeaf() + this%numRoot() + &
                                 myPanicleID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%Panicle(myPanicleID)%femdomain,&
                            DomainID   = this%numStem() + this%numLeaf() + this%numRoot() + &
                                 myPanicleID    ,& 
                            MyDomainID = yourStemID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif
    


    if(present(debug) )then
        if(debug)then
            print *, "[ok] overset >> done."        
        endif
    endif



    call solver%init(NumDomain=this%numStem() +this%numLeaf() + this%numRoot() &
        + this%numPanicle() )
    
    FEMDomainPointers = this%getFEMDomainPointers()
    call solver%setDomain(FEMDomainPointers=FEMDomainPointers )
    
    if(present(debug) )then
        if(debug)then
            print *, "[ok] initSolver >> done."        
        endif
    endif

    call solver%setCRS(DOF=3,debug=debug)

    ! CRS ready!

    if( .not. this%checkYoungModulus())then
        print *, "[ERROR] YoungModulus(:) is not ready."
        stop
    endif
    if( .not. this%checkPoissonRatio())then
        print *, "[ERROR] PoissonRatio(:) is not ready."
        stop
    endif
    if( .not. this%checkDensity())then
        print *, "[ERROR] Density(:) is not ready."
        stop
    endif


    if(present(debug) )then
        if(debug)then
            print *, "[ok] setCRS >> done."        
        endif
    endif
    
    !$OMP parallel 
    !$OMP do
    do DomainID=1,size(FEMDomainPointers)
        do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne()
            call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
               Matrix=FEMDomainPointers(DomainID)%femdomainp%StiffnessMatrix(&
               ElementID=ElementID,&
               E=this%getYoungModulus(DomainID=DomainID,ElementID=ElementID), &
               v=this%getPoissonRatio(DomainID=DomainID,ElementID=ElementID)  ) )
        enddo
    enddo
    !$OMP end do
    !$OMP end parallel
    
    
    if(present(debug) )then
        if(debug)then
            print *, "[ok] set Matrix & vectors >> done."        
        endif
    endif
    
    call solver%setEbOM(penalty=input(default=10000000.0d0,option=penalty), DOF=3)

    if(present(debug) )then
        if(debug)then
            print *, "[ok] set EbOM >> done."        
        endif
    endif
    call solver%keepThisMatrixAs("A")
    !call solver%saveMatrix(name="A",CRS_as_dense=.true.,zero_or_nonzero=.true)
    call solver%zeros()
    
    ! mass matrix
    !$OMP parallel 
    !$OMP do
    do DomainID=1,size(FEMDomainPointers)
        do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne()
            call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
               Matrix=FEMDomainPointers(DomainID)%femdomainp%massMatrix(&
                ElementID=ElementID,&
                Density=this%getDensity(DomainID=DomainID,ElementID=ElementID), &
                DOF=3 ) )
        enddo
    enddo
    !$OMP end do
    !$OMP end parallel
    call solver%keepThisMatrixAs("B")
    
    ! fix-boundary conditions
    nn_domains = 0
    do i=1,size(FEMDomainPointers)
        if(FEMDomainPointers(i)%FEMDomainp%z_min() <= ground_level )then
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-2 + nn_domains*3
            call solver%fix_eig(IDs=FixBoundary)
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-1 + nn_domains*3
            call solver%fix_eig(IDs=FixBoundary)
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-0 + nn_domains*3
            call solver%fix_eig(IDs=FixBoundary)
        endif
        nn_domains = nn_domains + FEMDomainPointers(i)%FEMDomainp%nn()
    enddo

    if(present(debug) )then
        if(debug)then
            print *, "[ok] FixBoundary >> done."        
        endif
    endif

    if(present(debug) )then
        solver%debug = debug
    endif



    call solver%eig(eigen_value=All_Frequency,eigen_vectors=All_EigenVectors)
    call solver%remove()

    ! simplify this part
    ! normalize EigenVectors
    do i=1,size(All_EigenVectors,2)
        vec_norm = norm(All_EigenVectors(:,i) )
        
        All_EigenVectors(:,i) = All_EigenVectors(:,i)/vec_norm
    enddo

    Frequency = zeros(num_freq)
    EigenVectors = zeros(size(All_EigenVectors,1),num_freq)

    do i=1,num_freq
        n = minvalID(All_Frequency)
        EigenVectors(:,i) = All_EigenVectors(:,n)
        Frequency(i)      = All_Frequency(n)
        All_Frequency(n) = maxval(All_Frequency) 
    enddo

    do i=1,size(Frequency)
        if(Frequency(i)<0.0d0)then
            Frequency(i)=0.0d0
        endif
    enddo
    Frequency = sqrt((Frequency))/(2.0d0*math%PI)



    
    if(present(debug) )then
        if(debug)then
            print *, "[ok] Solve >> done."        
        endif
    endif
    


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


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

    ret=0
    if(.not.allocated(this%leaf) )then
        return
    endif
    
    do i=1,size(this%leaf)
        if(this%leaf(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################

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

    ret=0
    if(.not.allocated(this%Stem) )then
        return
    endif
    
    do i=1,size(this%Stem)
        if(this%Stem(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################

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

    ret=0
    if(.not.allocated(this%Root) )then
        return
    endif
    
    do i=1,size(this%Root)
        if(this%Root(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################

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

    ret=0
    if(.not.allocated(this%Panicle) )then
        return
    endif
    
    do i=1,size(this%Panicle)
        if(this%Panicle(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################


! ################################################################
function getFEMDomainPointersRice(this) result(FEMDomainPointers)
    class(Rice_),target,intent(in) :: this
    type(FEMDomainp_),allocatable :: FEMDomainPointers(:)
    integer(int32) :: num_FEMDomain,i, n

    num_FEMDomain = this%numStem() + this%numLeaf() + this%numRoot() &
        + this%numPanicle() 
    allocate(FEMDomainPointers(num_FEMDomain) )
    n = 0
    do i=1,this%numStem()
        if(.not.this%stem(i)%femdomain%empty() )then
            n = n + 1
            FEMDomainPointers(n)%femdomainp => this%stem(i)%femdomain
        endif
    enddo

    do i=1,this%numLeaf()
        if(.not.this%leaf(i)%femdomain%empty() )then
            n = n + 1
            FEMDomainPointers(n)%femdomainp => this%leaf(i)%femdomain
        endif
    enddo
    do i=1,this%numRoot()
        if(.not.this%root(i)%femdomain%empty() )then
            n = n + 1
            FEMDomainPointers(n)%femdomainp => this%root(i)%femdomain
        endif
    enddo


    do i=1,this%numPanicle()
        if(.not.this%Panicle(i)%femdomain%empty() )then
            n = n + 1
            FEMDomainPointers(n)%femdomainp => this%Panicle(i)%femdomain
        endif
    enddo
end function
! ################################################################


! ################################################################
function checkYoungModulusRice(this) result(all_young_modulus_is_set)
    class(Rice_),intent(in) :: this
    logical :: all_young_modulus_is_set 
    integer(int32) :: i
    ! order: stem -> leaf -> root

    all_young_modulus_is_set = .true.
    do i=1,this%numStem()
        if(.not.allocated(this%stem(i)%YoungModulus) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkYoungModulusRice >> Young Modulus is not set"
            print *, "@ Stem ID:",i
            print *, "check it by: allocated(this%stem("+str(i)+")%YoungModulus)"
            return
        endif
    enddo

    do i=1,this%numLeaf()
        if(.not.allocated(this%Leaf(i)%YoungModulus) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkYoungModulusRice >> Young Modulus is not set"
            print *, "@ Leaf ID:",i
            print *, "check it by: allocated(this%Leaf("+str(i)+")%YoungModulus)"
            return
        endif
    enddo

    do i=1,this%numRoot()
        if(.not.allocated(this%Root(i)%YoungModulus) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkYoungModulusRice >> Young Modulus is not set"
            print *, "@ Root ID:",i
            print *, "check it by: allocated(this%Root("+str(i)+")%YoungModulus)"
            return
        endif
    enddo

    do i=1,this%numPanicle()
        if(.not.allocated(this%Panicle(i)%YoungModulus) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkYoungModulusRice >> Young Modulus is not set"
            print *, "@ Panicle ID:",i
            print *, "check it by: allocated(this%Panicle("+str(i)+")%YoungModulus)"
            return
        endif
    enddo
end function
! ################################################################




! ################################################################
function checkPoissonRatioRice(this) result(all_young_modulus_is_set)
    class(Rice_),intent(in) :: this
    logical :: all_young_modulus_is_set 
    integer(int32) :: i
    ! order: stem -> leaf -> root

    all_young_modulus_is_set = .true.
    do i=1,this%numStem()
        if(.not.allocated(this%stem(i)%PoissonRatio) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkPoissonRatioRice >> Young Modulus is not set"
            print *, "@ Stem ID:",i
            print *, "check it by: allocated(this%stem("+str(i)+")%PoissonRatio)"
            return
        endif
    enddo

    do i=1,this%numLeaf()
        if(.not.allocated(this%Leaf(i)%PoissonRatio) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkPoissonRatioRice >> Young Modulus is not set"
            print *, "@ Leaf ID:",i
            print *, "check it by: allocated(this%Leaf("+str(i)+")%PoissonRatio)"
            return
        endif
    enddo

    do i=1,this%numRoot()
        if(.not.allocated(this%Root(i)%PoissonRatio) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkPoissonRatioRice >> Young Modulus is not set"
            print *, "@ Root ID:",i
            print *, "check it by: allocated(this%Root("+str(i)+")%PoissonRatio)"
            return
        endif
    enddo

    do i=1,this%numPanicle()
        if(.not.allocated(this%Panicle(i)%PoissonRatio) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkPoissonRatioRice >> Young Modulus is not set"
            print *, "@ Panicle ID:",i
            print *, "check it by: allocated(this%Panicle("+str(i)+")%PoissonRatio)"
            return
        endif
    enddo
end function
! ################################################################

! ################################################################
function checkDensityRice(this) result(all_young_modulus_is_set)
    class(Rice_),intent(in) :: this
    logical :: all_young_modulus_is_set 
    integer(int32) :: i
    ! order: stem -> leaf -> root

    all_young_modulus_is_set = .true.
    do i=1,this%numStem()
        if(.not.allocated(this%stem(i)%Density) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkDensityRice >> Young Modulus is not set"
            print *, "@ Stem ID:",i
            print *, "check it by: allocated(this%stem("+str(i)+")%Density)"
            return
        endif
    enddo

    do i=1,this%numLeaf()
        if(.not.allocated(this%Leaf(i)%Density) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkDensityRice >> Young Modulus is not set"
            print *, "@ Leaf ID:",i
            print *, "check it by: allocated(this%Leaf("+str(i)+")%Density)"
            return
        endif
    enddo

    do i=1,this%numRoot()
        if(.not.allocated(this%Root(i)%Density) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkDensityRice >> Young Modulus is not set"
            print *, "@ Root ID:",i
            print *, "check it by: allocated(this%Root("+str(i)+")%Density)"
            return
        endif
    enddo



    do i=1,this%numPanicle()
        if(.not.allocated(this%Panicle(i)%Density) )then
            all_young_modulus_is_set = .false.
            print *, "[!Warning!] checkDensityRice >> Young Modulus is not set"
            print *, "@ Panicle ID:",i
            print *, "check it by: allocated(this%Panicle("+str(i)+")%Density)"
            return
        endif
    enddo

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

function getYoungModulusFieldRice(this) result(YoungModulus)
    class(Rice_),intent(inout) :: this
    real(real64),allocatable :: YoungModulus(:)
    integer(int32),allocatable :: ElementList(:,:)
    integer(int32) :: TYPE_IDX, DOMAIN_IDX, ELEMENT_IDX, i

    ElementList = this%getElementList()
    YoungModulus = zeros(size(ElementList,1))

    do i=1,size(ElementList,1)
        TYPE_IDX = ElementList(i,1)
        DOMAIN_IDX = ElementList(i,2)
        ELEMENT_IDX = ElementList(i,3)
        if(TYPE_IDX == this%TYPE_STEM)then
            YoungModulus(i) = this%stem(DOMAIN_IDX)%YoungModulus(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_LEAF)then
            YoungModulus(i) = this%LEAF(DOMAIN_IDX)%YoungModulus(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_ROOT)then
            YoungModulus(i) = this%ROOT(DOMAIN_IDX)%YoungModulus(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_PANICLE)then
            YoungModulus(i) = this%PANICLE(DOMAIN_IDX)%YoungModulus(ELEMENT_IDX)
        endif
    enddo
    
end function

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

! ################################################################
function getPoissonRatioFieldRice(this) result(PoissonRatio)
    class(Rice_),intent(inout) :: this
    real(real64),allocatable :: PoissonRatio(:)
    integer(int32),allocatable :: ElementList(:,:)
    integer(int32) :: TYPE_IDX, DOMAIN_IDX, ELEMENT_IDX, i

    ElementList = this%getElementList()
    PoissonRatio = zeros(size(ElementList,1))

    do i=1,size(ElementList,1)
        TYPE_IDX = ElementList(i,1)
        DOMAIN_IDX = ElementList(i,2)
        ELEMENT_IDX = ElementList(i,3)
        if(TYPE_IDX == this%TYPE_STEM)then
            PoissonRatio(i) = this%stem(DOMAIN_IDX)%PoissonRatio(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_LEAF)then
            PoissonRatio(i) = this%LEAF(DOMAIN_IDX)%PoissonRatio(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_ROOT)then
            PoissonRatio(i) = this%ROOT(DOMAIN_IDX)%PoissonRatio(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_PANICLE)then
            PoissonRatio(i) = this%PANICLE(DOMAIN_IDX)%PoissonRatio(ELEMENT_IDX)
        endif
    enddo
    
end function

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

! ################################################################
function getDensityFieldRice(this) result(Density)
    class(Rice_),intent(inout) :: this
    real(real64),allocatable :: Density(:)
    integer(int32),allocatable :: ElementList(:,:)
    integer(int32) :: TYPE_IDX, DOMAIN_IDX, ELEMENT_IDX, i

    ElementList = this%getElementList()
    Density = zeros(size(ElementList,1))

    do i=1,size(ElementList,1)
        TYPE_IDX = ElementList(i,1)
        DOMAIN_IDX = ElementList(i,2)
        ELEMENT_IDX = ElementList(i,3)
        if(TYPE_IDX == this%TYPE_STEM)then
            Density(i) = this%stem(DOMAIN_IDX)%Density(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_LEAF)then
            Density(i) = this%LEAF(DOMAIN_IDX)%Density(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_ROOT)then
            Density(i) = this%ROOT(DOMAIN_IDX)%Density(ELEMENT_IDX)
        elseif(TYPE_IDX == this%TYPE_PANICLE)then
            Density(i) = this%PANICLE(DOMAIN_IDX)%Density(ELEMENT_IDX)
        endif
    enddo
    
end function

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

function getYoungModulusRice(this,DomainID,ElementID) result(YoungModulus)
    class(Rice_),intent(in) :: this
    integer(int32),optional,intent(in) :: DomainID, ElementID
    real(real64) :: YoungModulus 
    integer(int32) :: i, n
    
    if(DomainID > this%numStem() + this%numLeaf() + this%numRoot() &
        +  this%numPanicle() )then
        print *, "ERROR :: getYoungModulusRice >>  DomainID exceeds max_domain_size"
        return
    endif

    ! default >> search @ all domains
    ! order: stem -> leaf -> root
    if(1 <= DomainID &
        .and. DomainID <= this%numStem() )then
        n = DomainID - 0    
        YoungModulus = this%stem(n)%YoungModulus(ElementID)
        return
    elseif(this%numStem() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  )then
        n = DomainID - this%numStem()
        YoungModulus = this%leaf(n)%YoungModulus(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot()  )then
        n = DomainID - (this%numStem() + this%numLeaf())
        YoungModulus = this%root(n)%YoungModulus(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + this%numRoot() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot() &
        + this%numPanicle()  )then
        
        n = DomainID - (this%numStem() + this%numLeaf() + this%numRoot())
        YoungModulus = this%panicle(n)%YoungModulus(ElementID)
        return
    else
        print *, "[ERROR] >> getYoungModulusRice >> Invalid DomainID",DomainID
        return
    endif

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


! ################################################################
function getPoissonRatioRice(this,DomainID,ElementID) result(PoissonRatio)
    class(Rice_),intent(in) :: this
    integer(int32),intent(in) :: DomainID, ElementID
    real(real64) :: PoissonRatio 
    integer(int32) :: i, n
    
    if(DomainID > this%numStem() + this%numLeaf() + this%numRoot() &
        + this%numPanicle() )then
        print *, "ERROR :: getPoissonRatioRice >>  DomainID exceeds max_domain_size"
        return
    endif

    ! default >> search @ all domains
    ! order: stem -> leaf -> root
    if(1 <= DomainID &
        .and. DomainID <= this%numStem() )then
        n = DomainID - 0    
        PoissonRatio = this%stem(n)%PoissonRatio(ElementID)
        return
    elseif(this%numStem() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  )then
        n = DomainID - this%numStem()
        PoissonRatio = this%leaf(n)%PoissonRatio(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot()  )then
        n = DomainID - (this%numStem() + this%numLeaf())
        PoissonRatio = this%root(n)%PoissonRatio(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + this%numRoot() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot() &
        + this%numPanicle()  )then
        
        n = DomainID - (this%numStem() + this%numLeaf() + this%numRoot())
        PoissonRatio = this%panicle(n)%PoissonRatio(ElementID)
        return
    else
        print *, "[ERROR] >> getPoissonRatioRice >> Invalid DomainID",DomainID
        return
    endif

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


! ################################################################
function getDensityRice(this,DomainID,ElementID) result(Density)
    class(Rice_),intent(in) :: this
    integer(int32),intent(in) :: DomainID, ElementID
    real(real64) :: Density 
    integer(int32) :: i, n
    
    if(DomainID > this%numStem() + this%numLeaf() + this%numRoot() &
        +  this%numPanicle() )then
        print *, "ERROR :: getDensityRice >>  DomainID exceeds max_domain_size"
        return
    endif

    ! default >> search @ all domains
    ! order: stem -> leaf -> root
    if(1 <= DomainID &
        .and. DomainID <= this%numStem() )then
        n = DomainID - 0    
        Density = this%stem(n)%Density(ElementID)
        return
    elseif(this%numStem() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  )then
        n = DomainID - this%numStem()
        Density = this%leaf(n)%Density(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot()  )then
        n = DomainID - (this%numStem() + this%numLeaf())
        Density = this%root(n)%Density(ElementID)
        return
    elseif(this%numStem() + this%numLeaf() + this%numRoot() + 1 <= DomainID &
        .and. DomainID <= this%numStem() + this%numLeaf()  + this%numRoot() &
        + this%numPanicle()  )then
        
        n = DomainID - (this%numStem() + this%numLeaf() + this%numRoot())
        Density = this%panicle(n)%Density(ElementID)
        return
    else
        print *, "[ERROR] >> getDensityRice >> Invalid DomainID",DomainID
        return
    endif

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


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

subroutine deformRice(this,displacement) 
    class(Rice_),target,intent(inout) :: this

    ! deform Rice by displacement
    real(real64),optional,intent(in) :: displacement(:)

    type(FEMDomainp_),allocatable :: domainsp(:)
    integer(int32),allocatable :: contactList(:,:)
    integer(int32) :: i,j,numDomain,stemDomain,leafDomain,rootDomain,from,to,nd,nn
    real(real64) :: penalty,GLevel

    if(present(displacement) )then
        if(size(displacement)/=this%nn()*3 )then
            print *, "ERROR :: deformRice >> size(displacement) should be (this%numStem() + this%numLeaf() + this%numRoot())*3"
            return
        endif

        ! order :: stem -> leaf -> root
        from = 1
        to   = 0
        if(allocated(this%stem) )then
            do i=1,size(this%stem)
                if(.not. this%stem(i)%femdomain%Mesh%empty() )then
                    nn = this%stem(i)%femdomain%nn()
                    nd = this%stem(i)%femdomain%nd()

                    to = from + this%stem(i)%femdomain%nn()*this%stem(i)%femdomain%nd() -1

                    this%stem(i)%femdomain%mesh%nodcoord(:,:) = &
                    this%stem(i)%femdomain%mesh%nodcoord(:,:) + &
                    reshape(displacement(from:to),nn,nd )

                    from = to + 1
                endif
            enddo
        endif

        if(allocated(this%leaf) )then
            do i=1,size(this%leaf)
                if(.not. this%leaf(i)%femdomain%Mesh%empty() )then
                    nn = this%leaf(i)%femdomain%nn()
                    nd = this%leaf(i)%femdomain%nd()

                    to = from + this%leaf(i)%femdomain%nn()*this%leaf(i)%femdomain%nd() -1

                    this%leaf(i)%femdomain%mesh%nodcoord(:,:) = &
                    this%leaf(i)%femdomain%mesh%nodcoord(:,:) + &
                    reshape(displacement(from:to),nn,nd )

                    from = to + 1
                endif
            enddo
        endif

        if(allocated(this%root) )then
            do i=1,size(this%root)
                if(.not. this%root(i)%femdomain%Mesh%empty() )then
                    nn = this%root(i)%femdomain%nn()
                    nd = this%root(i)%femdomain%nd()

                    to = from + this%root(i)%femdomain%nn()*this%root(i)%femdomain%nd() -1

                    this%root(i)%femdomain%mesh%nodcoord(:,:) = &
                    this%root(i)%femdomain%mesh%nodcoord(:,:) + &
                    reshape(displacement(from:to),nn,nd )

                    from = to + 1
                endif
            enddo
        endif


        if(allocated(this%panicle) )then
            do i=1,size(this%panicle)
                if(.not. this%panicle(i)%femdomain%Mesh%empty() )then
                    nn = this%panicle(i)%femdomain%nn()
                    nd = this%panicle(i)%femdomain%nd()

                    to = from + this%panicle(i)%femdomain%nn()*this%panicle(i)%femdomain%nd() -1

                    this%panicle(i)%femdomain%mesh%nodcoord(:,:) = &
                    this%panicle(i)%femdomain%mesh%nodcoord(:,:) + &
                    reshape(displacement(from:to),nn,nd )

                    from = to + 1
                endif
            enddo
        endif
        
        
        return
    endif


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

function getElementListRice(this,x_min,x_max,y_min,y_max,z_min,z_max,debug) result(ElementList)
    class(Rice_),intent(inout) :: this
    integer(int32),allocatable :: ElementList(:,:)
    integer(int32),allocatable :: obj_type(:),obj_idx(:),elem_idx(:)
    real(real64),optional,intent(in) :: x_min,x_max,y_min,y_max,z_min,z_max
    logical,optional,intent(in) :: debug
    logical :: do_debug
    integer(int32) :: idx,n,m

    do_debug = input(default=.false.,option=debug)

    !ElementList(idx, [ObjType, ObjID, ElementID] )
    allocate(elem_idx(0) )
    allocate(obj_type(0) )
    allocate(obj_idx(0) )

    if(allocated(this%stem) )then
        do idx=1,size(this%stem)
            if(this%stem(idx)%femdomain%empty() )cycle
            m = size(elem_idx)
            elem_idx = &
                elem_idx // this%stem(idx)%femdomain%mesh%getElementList(&
                xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,zmin=z_min,zmax=z_max)

            obj_idx = obj_idx // idx*int(eyes( size(elem_idx)- m))
        enddo

        if(do_debug)then
            print *, "[o] STEM"
        endif
    else
        if(do_debug)then
            print *, "NO STEM"
        endif
    endif
    
    ! debug

    
    obj_type = obj_type // this%TYPE_STEM*int(eyes(size(elem_idx)))

    if(allocated(this%leaf) )then
        do idx=1,size(this%leaf)
            if(this%leaf(idx)%femdomain%empty() )cycle
            m = size(elem_idx)
            elem_idx = &
                elem_idx // this%leaf(idx)%femdomain%mesh%getElementList(&
                xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,zmin=z_min,zmax=z_max)
            obj_idx = obj_idx // idx*int(eyes( size(elem_idx)- m))
        enddo

        if(do_debug)then
            print *, "[o] LEAF"
        endif
    else
        if(do_debug)then
            print *, "NO LEAF"
        endif
    endif

    n = size(obj_type)
    obj_type = obj_type // this%TYPE_LEAF*int(eyes(size(elem_idx)-n))
    

    if(allocated(this%root) )then
        do idx=1,size(this%root)
            if(this%root(idx)%femdomain%empty() )cycle
            m = size(elem_idx)
            elem_idx = &
                elem_idx // this%root(idx)%femdomain%mesh%getElementList(&
                xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,zmin=z_min,zmax=z_max)
            obj_idx = obj_idx // idx*int(eyes( size(elem_idx)- m))
        enddo

        if(do_debug)then
            print *, "[o] ROOT"
        endif
    else
        if(do_debug)then
            print *, "NO ROOT"
        endif
    endif
    n = size(obj_type)
    obj_type = obj_type // this%TYPE_ROOT*int(eyes(size(elem_idx)-n))
    
    
    if(allocated(this%Panicle) )then
        do idx=1,size(this%Panicle)
            if(this%Panicle(idx)%femdomain%empty() )cycle
            m = size(elem_idx)
            elem_idx = &
                elem_idx // this%Panicle(idx)%femdomain%mesh%getElementList(&
                xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,zmin=z_min,zmax=z_max)
            obj_idx = obj_idx // idx*int(eyes( size(elem_idx)- m))
        enddo

        if(do_debug)then
            print *, "[o] PANICLE"
        endif
    else
        if(do_debug)then
            print *, "NO PANICLE"
        endif
    endif
    n = size(obj_type)
    obj_type = obj_type // this%TYPE_PANICLE*int(eyes(size(elem_idx)-n))
    
    ElementList = zeros( size(elem_idx),3 )
    ElementList(:,1) = obj_type
    ElementList(:,2) = obj_idx
    ElementList(:,3) = elem_idx

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


! ################################################################
function getDisplacementRice(this, ground_level,penalty,debug,EbOM_Algorithm) result(displacement)
    class(Rice_),target,intent(inout) :: this
    real(real64),intent(in) :: ground_level
    real(real64),optional,intent(in) :: penalty
    logical,optional,intent(in) :: debug
    real(real64),allocatable :: Frequency(:)
    character(*),optional,intent(in) :: EbOM_Algorithm
    !integer(int32),optional,intent(in) :: num_mode
    
    real(real64),allocatable :: displacement(:)

    type(FEMDomainp_),allocatable :: FEMDomainPointers(:)
    type(FEMSolver_) :: solver
    type(Math_) :: math

    real(real64),allocatable :: EigenVectors(:,:),buf(:,:),buf_vec(:)

    integer(int32) :: stem_id, leaf_id, root_id,DomainID,ElementID,i,n
    integer(int32) :: myStemID, yourStemID, myLeafID,myRootID, yourRootID
    integer(int32),allocatable :: FixBoundary(:)
    integer(int32) :: nn_domains,EbOM_Algorithm_id
    real(real64) :: vec_norm

    integer(int32) :: myPanicleID

    EbOM_Algorithm_id = FEMDomain_Overset_GPP
    if(present(EbOM_Algorithm) )then
        if(EbOM_Algorithm=="P2P")then
            EbOM_Algorithm_id=FEMDomain_Overset_P2P
        elseif(EbOM_Algorithm=="GPP")then
            EbOM_Algorithm_id=FEMDomain_Overset_P2P
        endif
    endif

    ! linear elasticity with infinitesimal strain theory
    n = this%numStem() + this%numLeaf() + this%numRoot() + this%numPanicle()
    
    allocate(FEMDomainPointers(n) )
    ! ORDER 
    ! [STEM] => [LEAF] => [ROOT] => [PANICLE]
    !(1) >> compute overset
    ! For stems
    if(allocated(this%stem2stem) )then
        if(allocated(this%stem) )then
            do myStemID = 1,size(this%stem2stem,1)
                do yourStemID = 1, size(this%stem2stem,2)
                    if(this%stem2stem(myStemID,yourStemID)>=1 )then
                        ! connected
                        call this%stem(myStemID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = myStemID  ,&
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%stem(myStemID)%femdomain,&
                            DomainID   = myStemID    ,& 
                            MyDomainID = yourStemID  ,&
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif

    if(allocated(this%leaf2stem) )then
        if(allocated(this%leaf) .and. allocated(this%stem) )then
            do myLeafID = 1,size(this%leaf2stem,1)
                do yourStemID = 1, size(this%leaf2stem,2)
                    if(this%leaf2stem(myLeafID,yourStemID)>=1 )then
                        ! connected
                        call this%leaf(myLeafID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() + myLeafID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%leaf(myLeafID)%femdomain,&
                            DomainID   = this%numStem() + myLeafID    ,& 
                            MyDomainID = yourStemID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif
    



    if(allocated(this%root2stem) )then
        if(allocated(this%stem) .and. allocated(this%root) )then
            do myRootID = 1,size(this%root2stem,1)
                do yourStemID = 1, size(this%root2stem,2)
                    if(this%root2stem(myRootID,yourStemID)>=1 )then
                        ! connected
                        call this%root(myRootID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() +this%numLeaf() + myRootID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=  this%root(myRootID)%femdomain,&
                            DomainID   = this%numStem() +this%numLeaf() + myRootID    ,& 
                            MyDomainID =  yourStemID , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif


    if(allocated(this%root2root) )then
        if(allocated(this%root) )then
            do myRootID = 1,size(this%root2root,1)
                do yourrootID = 1, size(this%root2root,2)
                    if(this%root2root(myRootID,yourrootID)>=1 )then
                        ! connected
                        call this%root(myRootID)%femdomain%overset(&
                            FEMDomain=this%root(yourrootID)%femdomain,&
                            DomainID   = this%numroot() +this%numLeaf() + yourrootID    ,& 
                            MyDomainID = this%numroot() +this%numLeaf() + myRootID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    
                        call this%root(yourrootID)%femdomain%overset(&
                            FEMDomain=this%root(myRootID)%femdomain,&
                            DomainID   = this%numroot() +this%numLeaf() + myRootID    ,& 
                            MyDomainID =  this%numroot() +this%numLeaf() + yourrootID , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif

    
    if(allocated(this%Panicle2stem) )then
        if(allocated(this%Panicle) .and. allocated(this%stem) )then
            do myPanicleID = 1,size(this%Panicle2stem,1)
                do yourStemID = 1, size(this%Panicle2stem,2)
                    if(this%Panicle2stem(myPanicleID,yourStemID)>=1 )then
                        ! connected
                        call this%Panicle(myPanicleID)%femdomain%overset(&
                            FEMDomain=this%stem(yourStemID)%femdomain,&
                            DomainID   = yourStemID    ,& 
                            MyDomainID = this%numStem() + this%numLeaf() + this%numRoot() + &
                                 myPanicleID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                        
                        call this%stem(yourStemID)%femdomain%overset(&
                            FEMDomain=this%Panicle(myPanicleID)%femdomain,&
                            DomainID   = this%numStem() + this%numLeaf() + this%numRoot() + &
                             myPanicleID    ,& 
                            MyDomainID = yourStemID  , &
                            algorithm=EbOM_Algorithm_id ) ! or "P2P"
                    endif
                enddo
            enddo
        endif
    endif
    


    if(present(debug) )then
        if(debug)then
            print *, "[ok] overset >> done."        
        endif
    endif



    call solver%init(NumDomain=this%numStem() +this%numLeaf() + this%numRoot() &
        + this%numPanicle() )
    
    FEMDomainPointers = this%getFEMDomainPointers()
    call solver%setDomain(FEMDomainPointers=FEMDomainPointers )
    
    if(present(debug) )then
        if(debug)then
            print *, "[ok] initSolver >> done."        
        endif
    endif

    call solver%setCRS(DOF=3,debug=debug)

    ! CRS ready!

    if( .not. this%checkYoungModulus())then
        print *, "[ERROR] YoungModulus(:) is not ready."
        stop
    endif
    if( .not. this%checkPoissonRatio())then
        print *, "[ERROR] PoissonRatio(:) is not ready."
        stop
    endif
    if( .not. this%checkDensity())then
        print *, "[ERROR] Density(:) is not ready."
        stop
    endif


    if(present(debug) )then
        if(debug)then
            print *, "[ok] setCRS >> done."        
        endif
    endif
    
    !$OMP parallel 
    !$OMP do
    do DomainID=1,size(FEMDomainPointers)
        do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne()
            call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
                Matrix=FEMDomainPointers(DomainID)%femdomainp%StiffnessMatrix(&
                ElementID=ElementID,&
                E=this%getYoungModulus(DomainID=DomainID,ElementID=ElementID), &
                v=this%getPoissonRatio(DomainID=DomainID,ElementID=ElementID)  ) )
            call solver%setVector(DomainID=DomainID,ElementID=ElementID,DOF=3,&
                Vector=FEMDomainPointers(DomainID)%femdomainp%massVector(&
                ElementID=ElementID,&
                Density=this%getDensity(DomainID=DomainID,ElementID=ElementID), &
                Accel=[0.0d0, 0.0d0, -this%Gravity_acceralation], &
                DOF=3 ) )
        enddo
    enddo
    !$OMP end do
    !$OMP end parallel
    
    
    if(present(debug) )then
        if(debug)then
            print *, "[ok] set Matrix & vectors >> done."        
        endif
    endif
    
    call solver%setEbOM(penalty=input(default=10000000.0d0,option=penalty), DOF=3)

    if(present(debug) )then
        if(debug)then
            print *, "[ok] set EbOM >> done."        
        endif
    endif
    
    ! mass matrix
    !$OMP parallel 
    !$OMP do
    do DomainID=1,size(FEMDomainPointers)
        do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne()
            call solver%setMatrix(DomainID=DomainID,ElementID=ElementID,DOF=3,&
               Matrix=FEMDomainPointers(DomainID)%femdomainp%massMatrix(&
                ElementID=ElementID,&
                Density=this%getDensity(DomainID=DomainID,ElementID=ElementID), &
                DOF=3 ) )
        enddo
    enddo
    !$OMP end do
    !$OMP end parallel
    
    ! fix-boundary conditions
    nn_domains = 0
    do i=1,size(FEMDomainPointers)
        if(FEMDomainPointers(i)%FEMDomainp%z_min() <= ground_level )then
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-2 !+ nn_domains*3
            call solver%fix(DomainID=i,IDs=FixBoundary,FixValue=0.0d0)
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-1 !+ nn_domains*3
            call solver%fix(DomainID=i,IDs=FixBoundary,FixValue=0.0d0)
            FixBoundary = FEMDomainPointers(i)%FEMDomainp%select(z_max = ground_level )*3-0 !+ nn_domains*3
            call solver%fix(DomainID=i,IDs=FixBoundary,FixValue=0.0d0)
        endif
        !nn_domains = nn_domains + FEMDomainPointers(i)%FEMDomainp%nn()
    enddo

    if(present(debug) )then
        if(debug)then
            print *, "[ok] FixBoundary >> done."        
        endif
    endif

    if(present(debug) )then
        solver%debug = debug
    endif



    displacement = solver%solve(algorithm="BiCGSTAB")
    call solver%remove()

    if(present(debug) )then
        if(debug)then
            print *, "[ok] Solve >> done."        
        endif
    endif
    

end function

! ################################################################
function getStressFieldRice(this,displacement,i,j,option) result(StressField)
    class(Rice_),intent(inout) :: this
    real(real64),intent(in) :: displacement(:)
    integer(int32),optional,intent(in) :: i,j
    character(*),optional,intent(in) :: option

    real(real64),allocatable :: StressField(:)
    integer(int32) :: ii,jj, n, obj_idx

    StressField = zeros(0)
    n = 1
    if(allocated(this%stem) )then
        do obj_idx=1,size(this%stem)
            if(this%stem(obj_idx)%femdomain%mesh%empty() ) cycle
            StressField = StressField // &
                this%stem(obj_idx)%femdomain%getElementCauchyStress(&
                displacement=displacement(n:n+this%stem(obj_idx)%femdomain%nn()&
                    *this%stem(obj_idx)%femdomain%nd()-1 ),&
                E = this%stem(obj_idx)%YoungModulus(:),&
                v = this%stem(obj_idx)%PoissonRatio(:) ,i=i,j=j,option=option)
                n = n + this%stem(obj_idx)%femdomain%nn()&
                *this%stem(obj_idx)%femdomain%nd()
        enddo
    endif


    if(allocated(this%leaf) )then
        do obj_idx=1,size(this%leaf)
            if(this%leaf(obj_idx)%femdomain%mesh%empty() ) cycle
            StressField = StressField // &
                this%leaf(obj_idx)%femdomain%getElementCauchyStress(&
                displacement=displacement(n:n+this%leaf(obj_idx)%femdomain%nn()&
                *this%leaf(obj_idx)%femdomain%nd()-1 ),&
                E = this%leaf(obj_idx)%YoungModulus(:),&
                v = this%leaf(obj_idx)%PoissonRatio(:) ,i=i,j=j,option=option)
                n = n + this%leaf(obj_idx)%femdomain%nn()&
                *this%leaf(obj_idx)%femdomain%nd()
        enddo
    endif

    if(allocated(this%root) )then
        do obj_idx=1,size(this%root)
            if(this%root(obj_idx)%femdomain%mesh%empty() ) cycle
            StressField = StressField // &
                this%root(obj_idx)%femdomain%getElementCauchyStress(&
                displacement=displacement(n:n+this%root(obj_idx)%femdomain%nn()&
                *this%root(obj_idx)%femdomain%nd()-1 ),&
                E = this%root(obj_idx)%YoungModulus(:),&
                v = this%root(obj_idx)%PoissonRatio(:) ,i=i,j=j,option=option)
                n = n + this%root(obj_idx)%femdomain%nn()&
                *this%root(obj_idx)%femdomain%nd()
        enddo
    endif


    if(allocated(this%panicle) )then
        do obj_idx=1,size(this%panicle)
            if(this%panicle(obj_idx)%femdomain%mesh%empty() ) cycle
            StressField = StressField // &
                this%panicle(obj_idx)%femdomain%getElementCauchyStress(&
                displacement=displacement(n:n+this%panicle(obj_idx)%femdomain%nn()&
                *this%panicle(obj_idx)%femdomain%nd()-1 ),&
                E = this%panicle(obj_idx)%YoungModulus(:),&
                v = this%panicle(obj_idx)%PoissonRatio(:) ,i=i,j=j,option=option)
                n = n + this%panicle(obj_idx)%femdomain%nn()&
                *this%panicle(obj_idx)%femdomain%nd()
        enddo
    endif




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

subroutine export_eigRice(this,name,Frequency,ModeVectors,stress_type)
    class(Rice_),intent(inout) :: this
    character(*),intent(in) :: Name
    character(*),optional,intent(in) :: stress_type
    real(real64),intent(in) :: Frequency(:),ModeVectors(:,:)
    real(real64),allocatable :: displacement(:),stress(:)
    integer(int32) :: i,j 
    type(IO_) :: f

    call f%open(name + ".csv","w")
    call f%write("# Mode, Eigenfrequency (Hz)")
    do i=1,10
        displacement = ModeVectors(:,i)
        do j=1,36
            call this%deform(displacement =  cos(radian(j*10.0d0) ) * displacement )

            if(present(stress_type) )then
                stress = this%getStressField(Displacement=cos(radian(j*10.0d0) ) * displacement,option=stress_type)
                call this%vtk(name + zfill(i,3)+"_"+zfill(j,4),single_file = .true.,scalar_field=stress)
            else
                call this%vtk(name + zfill(i,3)+"_"+zfill(j,4),single_file = .true.)
            endif
            call this%deform(displacement = -cos(radian(j*10.0d0) ) *  displacement)
        enddo
        write(f%fh,*) str(i) +" , " , Frequency(i)
    enddo
    call f%close()
    
end subroutine

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


function add_rice_to_rice(rice1,rice2) result(ret)
    type(Rice_),intent(in) :: rice1,rice2
    type(Rice_) :: ret

    ret%num_leaf = rice1%num_leaf + rice2%num_leaf
    ret%num_stem = rice1%num_stem + rice2%num_stem
    ret%num_panicle = rice1%num_panicle + rice2%num_panicle
    ret%num_root = rice1%num_root + rice2%num_root


    ret%leaf = rice1%leaf // rice2%leaf
    ret%stem = rice1%stem // rice2%stem
    ret%panicle = rice1%panicle // rice2%panicle
    !ret%root = rice1%root // rice2%root
    
    ret%leaf2stem = rice1%leaf2stem .diag. rice2%leaf2stem
    ret%stem2stem = rice1%stem2stem .diag. rice2%stem2stem
    ret%panicle2stem = rice1%panicle2stem .diag. rice2%panicle2stem
    !ret%root2stem = rice1%root2stem .diag. rice2%root2stem
    !et%root2root = rice1%root2root .diag. rice2%root2root

end function



subroutine substitute_rice_into_rice(ret,rice1)
    type(Rice_),intent(out) :: ret
    type(Rice_),intent(in) :: rice1
    

    ret%num_leaf = rice1%num_leaf 
    ret%num_stem = rice1%num_stem 
    ret%num_panicle = rice1%num_panicle 
    ret%num_root = rice1%num_root 

    if(allocated(rice1%leaf) )then
        ret%leaf = rice1%leaf
    endif
    if(allocated(rice1%stem) )then
        ret%stem = rice1%stem
    endif
    if(allocated(rice1%panicle) )then
        ret%panicle = rice1%panicle
    endif
    if(allocated(rice1%root) )then
        ret%root = rice1%root
    endif
    
    if(allocated(rice1%leaf2stem ))then
        ret%leaf2stem = rice1%leaf2stem 
    endif
    if(allocated(rice1%stem2stem ))then
        ret%stem2stem = rice1%stem2stem 
    endif
    if(allocated(rice1%panicle2stem ))then
        ret%panicle2stem = rice1%panicle2stem 
    endif
    if(allocated(rice1%root2stem ))then
        ret%root2stem = rice1%root2stem 
    endif
    if(allocated(rice1%root2root ))then
        ret%root2root = rice1%root2root 
    endif

end subroutine

end module RiceClass