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

    type :: soybean_NodeID_Branch_
        integer(int32),allocatable :: ID(:)
    end type
    
    integer(int32),parameter :: PF_DEFAULT_SOYBEAN_ASIZE=300

    type :: soybean_
        ! growth_habit = determinate, indeterminate, semi-indeterminate, or vine
        character*20 :: growth_habit
        character*2  :: growth_stage
        integer(int32) :: Num_Of_Node
        integer(int32) :: num_leaf
        integer(int32) :: num_stem_node
        integer(int32) :: Num_Of_Root
        
        integer(int32) :: MaxLeafNum= PF_DEFAULT_SOYBEAN_ASIZE
        integer(int32) :: MaxRootNum= PF_DEFAULT_SOYBEAN_ASIZE
        integer(int32) :: MaxStemNum= PF_DEFAULT_SOYBEAN_ASIZE

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

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

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

        real(real64)    :: leaf_angle_ave(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_angle_sig(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_length_ave(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_length_sig(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_width_ave(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_width_sig(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_thickness_ave(PF_DEFAULT_SOYBEAN_ASIZE*3)
        real(real64)    :: leaf_thickness_sig(PF_DEFAULT_SOYBEAN_ASIZE*3)
        
        character(3) :: Stage ! VE, CV, V1,V2, ..., R1, R2, ..., R8
        character(200) :: name
        integer(int32)::stage_id=0
        real(real64) :: dt
        type(Seed_) :: Seed
        type(PlantNode_),allocatable :: NodeSystem(:)
        type(PlantRoot_),allocatable :: RootSystem(:)

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

        ! material info
        real(real64),allocatable :: stemYoungModulus(:)
        real(real64),allocatable :: leafYoungModulus(:)
        real(real64),allocatable :: rootYoungModulus(:)

        real(real64),allocatable :: stemPoissonRatio(:)
        real(real64),allocatable :: leafPoissonRatio(:)
        real(real64),allocatable :: rootPoissonRatio(:)

        real(real64),allocatable :: stemDensity(:)
        real(real64),allocatable :: leafDensity(:)
        real(real64),allocatable :: rootDensity(:)
        
        ! 節-節点データ構造
        type(Mesh_) :: struct 
        integer(int32),allocatable :: leaf2stem(:,:)
        integer(int32),allocatable :: stem2stem(:,:)
        integer(int32),allocatable :: root2stem(:,:)
        integer(int32),allocatable :: root2root(:,:)
        
        ! 器官オブジェクト配列
        type(FEMDomain_),allocatable :: leaf_list(:)
        type(FEMDomain_),allocatable :: stem_list(:)
        type(FEMDomain_),allocatable :: root_list(:)

        ! シミュレータ
        type(ContactMechanics_) :: contact
        real(real64) :: time
        real(real64) :: seed_length
        real(real64) :: seed_width
        real(real64) :: seed_height
        real(real64),allocatable :: stem_angle(:,:)
        real(real64),allocatable :: root_angle(:,:)
        real(real64),allocatable :: leaf_angle(:,:)

        character(200) :: stemconfig=" "
        character(200) :: rootconfig=" "
        character(200) :: leafconfig=" "

        
        integer(int32),allocatable :: NodeID_MainStem(:)
        type(soybean_NodeID_Branch_),allocatable :: NodeID_Branch(:)

    contains
        procedure,public :: addStem => addStemSoybean
        !procedure,public :: addRoot => addRootSoybean
        !procedure,public :: addLeaf => addLeafSoybean

        ! creation
        procedure,public :: Init => initsoybean
        procedure,public :: remove => removeSoybean
        procedure,public :: create => initsoybean
        procedure,public :: new => initsoybean
        procedure,public :: sowing => initsoybean
        procedure,public :: export => exportSoybean

        ! observation
        procedure,public :: stemlength => stemlengthSoybean
        procedure,public :: NumberOfBranch => NumberOfBranchSoybean
        

        ! operation
        procedure,public :: findApical => findApicalSoybean
        
        procedure,public :: grow => growSoybean
        procedure,public :: getVolume => getVolumeSoybean
        procedure,public :: getBioMass => getBioMassSoybean
        procedure,public :: getTotalWeight => getTotalWeightSoybean
        procedure,public :: getSubDomain => getSubDomainSoybean
        procedure,public :: getSubDomainType => getSubDomainTypeSoybean
        procedure,public :: setSubDomain => setSubDomainSoybean
        
        
        procedure,public :: resize => resizeSoybean
        procedure,public :: deform => deformSoybean
        
        ! visualization
        procedure,public :: show => showSoybean
        procedure,public :: gmsh => gmshSoybean
        procedure,public :: msh => mshSoybean
        procedure,public :: vtk => vtkSoybean
        procedure,public :: stl => stlSoybean
        procedure,public :: json => jsonSoybean
        
        ! get info
        !procedure,public :: properties => propertiesSoybean
        ! number of subdomain
        procedure,public :: ns => nsSoybean
        ! number of element
        procedure,public :: ne => neSoybean
        ! number of points
        procedure,public :: nn => nnSoybean
        procedure,public :: np => nnSoybean
        

        ! regacy/experimental
        procedure,public :: WaterAbsorption => WaterAbsorptionSoybean
        procedure,public :: move => moveSoybean

        procedure,public :: numleaf => numleafsoybean
        procedure,public :: numstem => numstemsoybean
        procedure,public :: numroot => numrootsoybean

        procedure,public :: laytracing => laytracingsoybean
        procedure,public :: SinkSourceFlow => SinkSourceFlowSoybean

        procedure,public :: update => updateSoybean
        procedure,public :: AddNode => AddNodeSoybean
    end type

    type :: SoybeanCanopy_
        real(real64) :: inter_row, intra_row
        type(soybean_),allocatable :: Canopy(:,:)
    end type

contains

! ########################################
recursive subroutine updateSoybean(obj,stem_id, root_id, leaf_id,debug)
    class(Soybean_),intent(inout) :: obj
    integer(int32),optional,intent(in) :: stem_id, root_id, leaf_id    
    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,mgn
    logical,optional,intent(in) :: debug
    ! update connectivity
    if(.not. allocated(obj%stem2stem ))then
        print *, "updateSoybean >> ERROR :: .not. allocated(obj%stem2stem )"
        return
    endif


    ! margin between subdomains
    
    itr_tol = 100
    itr=0

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

    ! stem to stem
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(obj%stem2stem,1)
            do j=1, size(obj%stem2stem,2)
                this_stem_id = j
                next_stem_id = i
                if(obj%stem2stem(i,j)/=0 .and. i /= j)then
                    ! this_stem_id ===>>> next_stem_id, connected!
                    x_B(:) = obj%stem(this_stem_id)%getCoordinate("B")
                    x_A(:) = obj%stem(next_stem_id)%getCoordinate("A")
                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call obj%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 *, "soybean % update >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "soybean % update >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error) == 0.0d0) exit
        last_error = error
    enddo

    ! root to root
    last_error = 1.0d0
    do 
        itr=itr+1
        error = 0.0d0
        do i=1, size(obj%root2root,1)
            do j=1, size(obj%root2root,2)
                this_root_id = j
                next_root_id = i
                if(obj%root2root(i,j)/=0 .and. i /= j)then
                    ! this_root_id ===>>> next_root_id, connected!
                    x_B(:) = obj%root(this_root_id)%getCoordinate("B")
                    x_A(:) = obj%root(next_root_id)%getCoordinate("A")
                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call obj%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 *, "soybean % update >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "soybean % 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(obj%leaf2stem,1)
            do j=1, size(obj%leaf2stem,2)
                this_stem_id = j
                next_leaf_id = i
                if(obj%leaf2stem(i,j)==1)then
                    ! this_stem_id ===>>> next_leaf_id, connected!
                    x_B(:) = obj%stem(this_stem_id)%getCoordinate("B")
                    x_A(:) = obj%leaf(next_leaf_id)%getCoordinate("A")
                    diff(:) = x_B(:) - x_A(:)
                    error = error + dot_product(diff,diff)
                    call obj%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 *, "soybean % update >> error :: ",error
            endif
        endif
        if(itr > itr_tol) then
            print *, "soybean % update >> ERROR :: not converged"
            stop
        endif
        
        if( abs(error) + abs(last_error) == 0.0d0) exit
        last_error = error
    enddo

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

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

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


    timeOpt = input(default=.false.,option=profiler)

    if(timeOpt) then
        call time%start()
    endif


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

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

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

    obj%brr_angle_ave(:)= 0.0d0
    obj%brr_angle_sig(:)=10.0d0
    obj%brr_angle_ave(1)=30.0d0
    obj%brr_angle_sig(1)=2.0d0
    
    obj%mr_angle_ave=0.0d0
    obj%mr_angle_sig=2.0d0
    ! peti
    ! is also stem
    
    obj%peti_size_ave(:) = 0.20d0
    obj%peti_size_sig(:) = 0.010d0

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

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

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

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

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

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

    
    if(timeOpt) then
        print *, "[1] set default values :: "
        call time%show()
    endif
    

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

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

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


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

    if(timeOpt) then
        print *, "[2] create default config "
        call time%show()
    endif

    line = soyconf%parse(trim(conf),key1="Genotype",key2="Dt1")
    if(index(line,"Dt1")/=0 )then
        obj%determinate=.False.
    else
        obj%determinate=.True.
    endif

    call soyconf%open(trim(conf))
    blcount=0
    do
        read(soyconf%fh,'(a)') line
        if(debug) print *, trim(line)
        if( adjustl(trim(line))=="{" )then
            blcount=1
            cycle
        endif
        if( adjustl(trim(line))=="}" )then
            exit
        endif
        
        if(blcount==1)then
            
            if(index(line,"Name")/=0)then
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%name
            endif

            if(index(line,"Mainstem")/=0)then
                do
                    read(soyconf%fh,'(a)') line
                    if(debug) print *, trim(line)
                    if( index(line,"}")/=0 )then
                        exit
                    endif
                    
                    if(index(line,"Length")/=0 )then
                        rmc=index(line,",")
                        if(rmc /= 0)then
                            line(rmc:rmc)=" "
                        endif
                        id = index(line,":")
                        read(line(id+1:),*) obj%ms_length
                    endif

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

                    
                
                enddo
            endif

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

                do
                    read(soyconf%fh,'(a)') line
                    if(debug) print *, trim(line)
                    if( index(line,"}")/=0 )then
                        exit
                    endif
                    
                    if(index(line,"Length")/=0 )then
                        rmc=index(line,",")
                        if(rmc /= 0)then
                            line(rmc:rmc)=" "
                        endif
                        id = index(line,":")
                        read(line(id+1:),*) obj%br_length(branch_id)
                    endif

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

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

            ! for roots

            if(index(line,"Mainroot")/=0)then
                do
                    read(soyconf%fh,'(a)') line
                    if(debug) print *, trim(line)
                    if( index(line,"}")/=0 )then
                        exit
                    endif
                    
                    if(index(line,"Length")/=0 )then
                        rmc=index(line,",")
                        if(rmc /= 0)then
                            line(rmc:rmc)=" "
                        endif
                        id = index(line,":")
                        read(line(id+1:),*) obj%mr_length
                    endif

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

                    
                
                enddo
            endif

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

                do
                    read(soyconf%fh,'(a)') line
                    if(debug) print *, trim(line)
                    if( index(line,"}")/=0 )then
                        exit
                    endif
                    
                    if(index(line,"Length")/=0 )then
                        rmc=index(line,",")
                        if(rmc /= 0)then
                            line(rmc:rmc)=" "
                        endif
                        id = index(line,":")
                        read(line(id+1:),*) obj%brr_length(branch_id)
                    endif

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

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


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

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

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


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


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


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


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

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

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

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


            ! for version 2020.11.24

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

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

            if(index(line,"ms_angle_ave")/=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%ms_angle_ave = readvalreal
            endif
            
            if(index(line,"ms_angle_sig")/=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%ms_angle_sig = readvalreal
            endif
            ! peti
            ! is also stem
            
            if(index(line,"peti_size_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_size_ave(:) = readvalreal
            endif
            
            if(index(line,"peti_size_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_size_sig(:) = readvalreal
            endif
            
            if(index(line,"peti_width_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_width_ave(:) = readvalreal
            endif
            
            if(index(line,"peti_width_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_width_sig(:) = readvalreal
            endif
            
            if(index(line,"peti_angle_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_angle_ave(:) = readvalreal
            endif
            
            if(index(line,"peti_angle_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%peti_angle_sig(:) = readvalreal
            endif
            ! leaf
            
            if(index(line,"leaf_length_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_length_ave(:) = readvalreal
            endif
            
            if(index(line,"leaf_length_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_length_sig(:) = readvalreal
            endif
            
            if(index(line,"leaf_width_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_width_ave(:) = readvalreal
            endif
            
            if(index(line,"leaf_width_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_width_sig(:) = readvalreal
            endif
            
            if(index(line,"leaf_thickness_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_thickness_ave(:) = readvalreal
            endif
            
            if(index(line,"leaf_thickness_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_thickness_sig(:) = readvalreal
            endif
            
            if(index(line,"leaf_angle_ave")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_angle_ave(:) = readvalreal
            endif
            
            if(index(line,"leaf_angle_sig")  /=0 )then
                
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) readvalreal
                obj%leaf_angle_sig(:) = readvalreal
            endif


            ! added in 2020/12/15
            ! for roots



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

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

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


            cycle

        endif

    enddo
    call soyconf%close()

    

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


    if(timeOpt) then
        print *, "[3] read config "
        call time%show()
    endif

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

        allocate(obj%leafYoungModulus(obj%MaxLeafNum) )
        allocate(obj%rootYoungModulus(obj%MaxrootNum) )
        allocate(obj%stemYoungModulus(obj%MaxstemNum) )
        ! default value
        obj%leafYoungModulus(:) = 1000.0d0
        obj%rootYoungModulus(:) = 1000.0d0
        obj%stemYoungModulus(:) = 1000.0d0
        
        allocate(obj%leafPoissonRatio(obj%MaxLeafNum) )
        allocate(obj%rootPoissonRatio(obj%MaxrootNum) )
        allocate(obj%stemPoissonRatio(obj%MaxstemNum) )
        obj%leafPoissonRatio(:) = 0.30d0
        obj%rootPoissonRatio(:) = 0.30d0
        obj%stemPoissonRatio(:) = 0.30d0
        
        allocate(obj%leafDensity(obj%MaxLeafNum) )
        allocate(obj%rootDensity(obj%MaxrootNum) )
        allocate(obj%stemDensity(obj%MaxstemNum) )

        obj%leafDensity(:) = 0.0d0
        obj%rootDensity(:) = 0.0d0
        obj%stemDensity(:) = 0.0d0

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

        ! set mainstem
        
        allocate(obj%NodeID_MainStem(obj%ms_node) )
        call stem%init(config=obj%stemconfig)
        do i=1,obj%ms_node

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


        if(timeOpt) then
            print *, "[4] created Main stem."
            call time%show()
        endif
        
        do i=1,obj%ms_node-1
            call obj%stem(i+1)%connect("=>",obj%stem(i))
            obj%stem2stem(i+1,i) = 1
        enddo

        ! set branches
        k=obj%ms_node
        allocate(obj%NodeID_Branch( size(obj%br_node) ) )
        do i=1,size(obj%br_node) ! num branch
            allocate( obj%NodeID_Branch(i)%ID(obj%br_node(i))  )
            do j=1, obj%br_node(i)
                k = k + 1
                !call obj%stem(k)%init(config=obj%stemconfig)
                obj%stem(k) = stem
                obj%NodeID_Branch(i)%ID(j) = k

                call obj%stem(k)%resize(&
                    x = obj%br_width(i), &
                    y = obj%br_width(i), &
                    z = obj%br_length(i)/dble(obj%br_node(i) ) &
                    )
                    
                call obj%stem(k)%rotate(&
                    x = radian(random%gauss(mu=obj%br_angle_ave(j),sigma=obj%br_angle_sig(j) )),  &
                    y = 0.0d0,  &
                    z = radian(360.0d0*random%random() )   &
                    )                
                
                if(j==1)then
                    call obj%stem(k)%connect("=>",obj%stem(obj%br_from(i)  ))
                    obj%stem2stem(k,obj%br_from(i) ) = 1
                else
                    call obj%stem(k)%connect("=>",obj%stem(k-1))
                    obj%stem2stem(k,k-1) = 1
                endif
                    
            enddo
        enddo
        
        if(timeOpt) then
            print *, "[4] created Branches."
            call time%show()
        endif




        ! peti and leaf
        obj%num_stem_node = k
        obj%num_leaf = 0
        ! bugfix 2021/08/18
        call leaf%init(config=obj%leafconfig,species=PF_GLYCINE_SOJA)
        do i=1, k
            ! ３複葉
            ! add peti
            obj%num_stem_node = obj%num_stem_node +1
            !call obj%stem(obj%num_stem_node)%init(config=obj%stemconfig)
            obj%stem(obj%num_stem_node) = stem
            call obj%stem(obj%num_stem_node)%resize(&
                x = random%gauss(mu=obj%peti_width_ave(i),sigma=obj%peti_width_sig(i)), &
                y = random%gauss(mu=obj%peti_width_ave(i),sigma=obj%peti_width_sig(i)), &
                z = random%gauss(mu=obj%peti_size_ave(i),sigma=obj%peti_size_sig(i)) &
                )
            
            call obj%stem(obj%num_stem_node)%rotate(&
                x = radian(random%gauss(mu=obj%peti_angle_ave(i),sigma=obj%peti_angle_sig(i) )),  &
                y = 0.0d0,  &
                z = radian(360.0d0*random%random() )   &
                )      
            call obj%stem(obj%num_stem_node)%connect("=>",obj%stem(i))
            !obj%leaf2stem(num_stem_node,i) = 1   
            obj%stem2stem(obj%num_stem_node,i) = 1            

            

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


        if(timeOpt) then
            print *, "[4] created Peti and Leaves."
            call time%show()
        endif

        ! set mainroot
        call root%init(obj%rootconfig)
        do i=1,obj%mr_node

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

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

        ! set branches
        k=obj%mr_node
        do i=1,size(obj%brr_node)
            do j=1, obj%brr_node(i)
                k = k + 1
                !call obj%root(k)%init(config=obj%rootconfig)
                obj%root(k) = root
                call obj%root(k)%resize(&
                    x = obj%mr_width, &
                    y = obj%mr_width, &
                    z = obj%mr_length/dble(obj%mr_node) &
                    )
                    
                call obj%root(k)%rotate(&
                    x = radian(random%gauss(mu=obj%brr_angle_ave(j),sigma=obj%brr_angle_sig(j) )),  &
                    y = 0.0d0,  &
                    z = radian(360.0d0*random%random() )   &
                    )                
                
                if(j==1)then
                    call obj%root(k)%connect("=>",obj%root(obj%brr_from(i)  ))
                    obj%root2root(k,obj%brr_from(i) ) = 1
                else
                    call obj%root(k)%connect("=>",obj%root(k-1))
                    obj%root2root(k,k-1) = 1
                endif
                    
            enddo
        enddo
        

        obj%stage = "V"//trim(str(obj%ms_node))


        if(timeOpt) then
            print *, "[4] create objects."
        call time%show()
    endif
        return
    else
        ! create leaf, root, stem
        allocate(obj%leaf(obj%MaxLeafNum) )
        allocate(obj%root(obj%MaxrootNum) )
        allocate(obj%stem(obj%MaxstemNum) )

        allocate(obj%leafYoungModulus(obj%MaxLeafNum) )
        allocate(obj%rootYoungModulus(obj%MaxrootNum) )
        allocate(obj%stemYoungModulus(obj%MaxstemNum) )
        ! default value
        obj%leafYoungModulus(:) = 1000.0d0
        obj%rootYoungModulus(:) = 1000.0d0
        obj%stemYoungModulus(:) = 1000.0d0
        
        allocate(obj%leafPoissonRatio(obj%MaxLeafNum) )
        allocate(obj%rootPoissonRatio(obj%MaxrootNum) )
        allocate(obj%stemPoissonRatio(obj%MaxstemNum) )
        obj%leafPoissonRatio(:) = 0.30d0
        obj%rootPoissonRatio(:) = 0.30d0
        obj%stemPoissonRatio(:) = 0.30d0
        
        allocate(obj%leafDensity(obj%MaxLeafNum) )
        allocate(obj%rootDensity(obj%MaxrootNum) )
        allocate(obj%stemDensity(obj%MaxstemNum) )

        obj%leafDensity(:) = 0.0d0
        obj%rootDensity(:) = 0.0d0
        obj%stemDensity(:) = 0.0d0

        allocate(obj%stem2stem(obj%MaxstemNum,obj%MaxstemNum) )
        allocate(obj%leaf2stem(obj%MaxstemNum,obj%MaxLeafNum) )
        allocate(obj%root2stem(obj%MaxrootNum,obj%MaxstemNum) )
        allocate(obj%root2root(obj%MaxrootNum,obj%MaxrootNum) )
        
        !allocate(obj%struct%NodCoord(4,3) )
        !allocate(obj%struct%ElemNod(3,2) )
        !allocate(obj%struct%ElemMat(3) )
        ! 子葉結節部=(0,0,0)
        !obj%struct%NodCoord(1,1:3) = 0.0d0
        call obj%leaf(1)%init(obj%leafconfig,species=PF_GLYCINE_SOJA)
        call obj%leaf(1)%rotate(x=radian(90.0d0),y=radian(90.0d0),z=radian(10.0d0) )
        call obj%leaf(2)%init(obj%leafconfig,species=PF_GLYCINE_SOJA)
        call obj%leaf(2)%rotate(x=radian(90.0d0),y=radian(90.0d0),z=radian(-10.0d0) )
        
        call obj%stem(1)%init(obj%stemconfig)
        call obj%stem(1)%rotate(x=radian(40.0d0) )
        
        call obj%stem(2)%init(obj%stemconfig)
        call obj%stem(2)%rotate(x=radian(80.0d0) )
    
        call obj%root(1)%init(obj%rootconfig)
        call obj%root(1)%fix(x=0.0d0,y=0.0d0,z=0.0d0)
        call obj%root(1)%rotate(x=radian(-60.0d0) )
    
        call obj%leaf(1)%connect("=>",obj%stem(1))
        obj%leaf2stem(1,1) = 1
        
        call obj%leaf(2)%connect("=>",obj%stem(1))
        obj%leaf2stem(2,1) = 1
        
        call obj%stem(2)%connect("=>",obj%stem(1))
        obj%stem2stem(2,1) = 1
        
        call obj%root(1)%connect("=>",obj%stem(1))
        obj%root2stem(1,1) = 1
        
        obj%stage = "VE"
        ! 初生葉結節部
        !obj%struct%NodCoord(2,1) = 0.0d0
        !obj%struct%NodCoord(2,2) = 0.0d0
        !obj%struct%NodCoord(2,3) = 1.0d0/20.0d0*obj%seed_height
        ! 地際部
        !obj%struct%NodCoord(3,1) = 1.0d0/4.0d0*obj%seed_length
        !obj%struct%NodCoord(3,2) = 0.0d0
        !obj%struct%NodCoord(3,3) = -1.0d0/3.0d0*obj%seed_height
        ! 根冠
        !obj%struct%NodCoord(4,1) = 1.0d0/2.0d0*obj%seed_length
        !obj%struct%NodCoord(4,2) = 0.0d0
        !obj%struct%NodCoord(4,3) = -1.0d0/2.0d0*obj%seed_height
    
        ! 子葉-初生葉節
        !obj%struct%ElemNod(1,1) = 1
        !obj%struct%ElemNod(1,2) = 2
        ! 地際-子葉節
        !obj%struct%ElemNod(2,1) = 3
        !obj%struct%ElemNod(2,2) = 1
        ! 地際-根冠節
        !obj%struct%ElemNod(3,1) = 3
        !obj%struct%ElemNod(3,2) = 4
    
        ! 子葉-初生葉節 stem: 1
        !obj%struct%ElemMat(1) = 1
        ! 地際-子葉節 stem: 1
        !obj%struct%ElemMat(2) = 1
        ! 地際-根冠節 primary root: -1
        !obj%struct%ElemMat(3) = -1
    
        ! FEメッシュを生成
        ! 領域を確保
    !    n = input(default=80,option=max_leaf_num)
    !    allocate(obj%leaf_list(n) )
    !    n = input(default=80,option=max_stem_num)
    !    allocate(obj%stem_list(n) )
    !    n = input(default=80,option=max_root_num)
    !    allocate(obj%root_list(n) )
    !
    !    ! 子葉のメッシュを生成
    !    call obj%leaf_list(1)%create(meshtype="Sphere3D",x_num=10,y_num=10,z_num=10,&
    !        x_len=obj%seed_length,y_len=obj%seed_width,z_len=obj%seed_height)
    !    call obj%leaf_list(1)%move(x=0.0d0,y=-0.50d0*obj%seed_width,z=-0.50d0*obj%seed_height)
    !
    !    call obj%leaf_list(2)%create(meshtype="Sphere3D",x_num=10,y_num=10,z_num=10,&
    !        x_len=obj%seed_length,y_len=obj%seed_width,z_len=obj%seed_height)
    !    call obj%leaf_list(2)%rotate(x=radian(180.0d0) )
    !    call obj%leaf_list(2)%move(x=0.0d0,y=-0.50d0*obj%seed_width,z=-0.50d0*obj%seed_height)
    !
    !
    !
    !    ! 子葉-初生葉節のメッシュを生成
    !    rot(:) = 0.0d0
    !    call obj%stem_list(1)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
    !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
    !    ! 節基部の節点ID
    !    node_id = obj%struct%ElemNod(1,1)
    !    ! 節先端部の節点ID
    !    node_id2= obj%struct%ElemNod(1,2)
    !    ! 節基部の位置ベクトル
    !    loc(:) = obj%struct%NodCoord( node_id  ,:)
    !    ! 節先端部までの方向ベクトル
    !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)  
    !    
    !    ! structの構造データにメッシュデータを合わせる。
    !    print *, obj%stem_list(1)%Mesh%BottomElemID
    !    print *, obj%stem_list(1)%Mesh%TopElemID
    !
    !    elemid = obj%stem_list(1)%Mesh%BottomElemID
    !    node_id = obj%stem_list(1)%Mesh%ElemNod(elemID,1)
    !    meshloc(:) = obj%stem_list(1)%Mesh%NodCoord(node_id,:)
    !
    !    elemid = obj%stem_list(1)%Mesh%TopElemID
    !    node_id = obj%stem_list(1)%Mesh%ElemNod(elemID,1)
    !    meshvec(:) = obj%stem_list(1)%Mesh%NodCoord(node_id,:)-meshloc(:)
    
        !print *, "loc",loc
        !print *, "meshloc",meshloc
        !print *, "vec",vec
        !print *, "meshvec",meshvec
        
    !    ! 節中央を原点へ
    !    call obj%stem_list(1)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0)
    !    
    !    print *, "loc",loc
    !    print *, "vec",vec
    !    print *, "rot",rot
    !    zaxis(:)=0.0d0
    !    zaxis(3)=obj%seed_length/5.0d0
    !    rot(:) = angles(zaxis,vec)
    !    call obj%stem_list(1)%move(x=loc(1),y=loc(2),z=loc(3) )
    !    call obj%stem_list(1)%rotate(x=0.0d0,y=0.0d0,z=0.0d0 )
    !!    
    !    
    !!    
    !
    !
    !    ! 地際-子葉節のメッシュを生成
    !    rot(:) = 0.0d0
    !    call obj%stem_list(2)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
    !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
    !    ! 節基部の節点ID
    !    node_id = obj%struct%ElemNod(2,1)
    !    ! 節先端部の節点ID
    !    node_id2= obj%struct%ElemNod(2,2)
    !    ! 節基部の位置ベクトル
    !    loc(:) = obj%struct%NodCoord( node_id  ,:)
    !    ! 節先端部までの方向ベクトル
    !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)  
    !    ! 節中央を原点へ
    !    call obj%stem_list(2)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0,&
    !        z=-obj%seed_length/8.0d0)
    !    zaxis(:)=0.0d0
    !    zaxis(3)=obj%seed_length/5.0d0
    !    rot(:) = angles(zaxis,vec)
    !    print *, "loc",loc
    !    print *, "vec",vec
    !    print *, "rot",rot
    !    !call obj%stem_list(2)%rotate(x=rot(1),y=rot(2),z=rot(3) )
    !    call obj%stem_list(2)%move(x=loc(1),y=loc(2),z=loc(3) )
    !    
    !
    !
    !    ! 地際-根冠節のメッシュ生成
    !    rot(:) = 0.0d0
    !    call obj%root_list(1)%create(meshtype="rectangular3D",x_num=5,y_num=5,z_num=10,&
    !        x_len=obj%seed_width/6.0d0,y_len=obj%seed_width/6.0d0,z_len=obj%seed_length/4.0d0)
    !    ! 節基部の節点ID
    !    node_id = obj%struct%ElemNod(3,1)
    !    ! 節先端部の節点ID
    !    node_id2= obj%struct%ElemNod(3,2)
    !    ! 節基部の位置ベクトル
    !    loc(:) = obj%struct%NodCoord( node_id  ,:)
    !    ! 節先端部までの方向ベクトル
    !    vec(:) =  obj%struct%NodCoord( node_id2 ,:) - obj%struct%NodCoord( node_id  ,:)  
    !    ! 節基部へ移動
    !    call obj%root_list(1)%move(x=-obj%seed_width/12.0d0,y=-obj%seed_width/12.0d0,&
    !        z=-obj%seed_length/8.0d0)
    !    call obj%root_list(1)%move(x=loc(1),y=loc(2),z=loc(3) )
    !    zaxis(:)=0.0d0
    !    zaxis(3)=obj%seed_length/5.0d0
    !    rot(:) = angles(zaxis,vec)
    !    !call obj%root_list(1)%rotate(x=rot(1),y=rot(2),z=rot(3) )
    !    print *, "loc",loc
    !    print *, "vec",vec
    !    print *, "rot",rot    
    endif




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

            loc(:)=0.0d0

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

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

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

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

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

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

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

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


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

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

            obj%time=0.0d0
            return
        endif
    endif


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

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

    obj%dt = dt

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

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

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

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

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

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


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

subroutine SinkSourceFlowSoybean(obj)
    class(Soybean_),intent(inout) :: obj
    type(DiffusionEq_) :: DiffusionEq

    !DiffusionEq%femdomain => obj%femdomain

    

end subroutine


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

    obj%time=obj%time+dt


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

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

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


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


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

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

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

    itr=itr+1
    ! export NodeSystem
    do i=1,size(obj%NodeSystem)
            
        if(present(FilePath) )then
            call obj%NodeSystem(i)%export(FileName=trim(FilePath)//"/Node.geo",objID=itr)
        else
            call obj%NodeSystem(i)%export(FileName=trim(FileName)//"_Node.geo",objID=itr)
        endif
        if(i==obj%num_of_node  )then
            exit
        endif
    enddo

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




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



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

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

! ########################################
subroutine showSoybean(obj,name)
    class(Soybean_),intent(inout) :: obj
    character(*),intent(in)::name

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

    call obj%struct%export(name=name)

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



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

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

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

    ret=0
    do i=1,size(obj%stem)
        if(obj%stem(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################

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

    ret=0
    do i=1,size(obj%root)
        if(obj%root(i)%femdomain%Mesh%empty() .eqv. .false. )then
            ret=ret+1
        endif
    enddo
    
end function
! ########################################


! ########################################
subroutine gmshSoybean(obj,name,num_threads)
    class(Soybean_),intent(inout) :: obj
    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(obj%stem)
        !if(obj%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            call obj%stem(i)%gmsh(name=trim(name)//"_stem"//trim(str(i)))
        !endif
    enddo
    !$OMP end do
    !$OMP end parallel

    !$OMP parallel num_threads(n) private(i)
    !$OMP do 

    do i=1,size(obj%root)
        !if(obj%root(i)%femdomain%mesh%empty() .eqv. .false. )then
            call obj%root(i)%gmsh(name=trim(name)//"_root"//trim(str(i)))
        !endif
    enddo
    !$OMP end do
    !$OMP end parallel

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


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


! ########################################
subroutine mshSoybean(obj,name,num_threads)
    class(Soybean_),intent(inout) :: obj
    character(*),intent(in) :: name
    integer(int32),optional,intent(in) :: num_threads
    integer(int32) :: i,n
    type(IO_) :: f
    ! index file
    call f%open(trim(name)//"_index.txt","w")

    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_stem"//trim(str(i))//".msh")
            endif
        enddo
    endif
    
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_root"//trim(str(i))//".msh")
            endif
        enddo
    endif

    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_leaf"//trim(str(i))//".msh")
            endif
        enddo
    endif
    call f%close()
    
    n = input(default=1,option=num_threads)
    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(obj%stem)
        !if(obj%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            call obj%stem(i)%msh(name=trim(name)//"_stem"//trim(str(i)))
        !endif
    enddo
    !$OMP end do
    !$OMP end parallel

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

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

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


! ########################################
subroutine vtkSoybean(obj,name,num_threads)
    class(Soybean_),intent(inout) :: obj
    character(*),intent(in) :: name
    type(IO_) :: f
    integer(int32),optional,intent(in) :: num_threads
    integer(int32) :: i, n

    n = input(default=1,option=num_threads)
    
    ! index file
    call f%open(trim(name)//"_index.txt","w")

    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_stem"//trim(str(i))//".vtk")
            endif
        enddo
    endif
    
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_root"//trim(str(i))//".vtk")
            endif
        enddo
    endif

    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_leaf"//trim(str(i))//".vtk")
            endif
        enddo
    endif
    call f%close()

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


    if(allocated(obj%root))then

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

        !$OMP end do
        !$OMP end parallel
    endif


    if(allocated(obj%leaf))then

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

    

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


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

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

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

! ########################################
subroutine stlSoybean(obj,name,num_threads)
    class(Soybean_),intent(inout) :: obj
    character(*),intent(in) :: name
    integer(int32),optional,intent(in) :: num_threads
    integer(int32) :: i,n

    type(IO_) :: f
    ! index file
    call f%open(trim(name)//"_index.txt","w")

    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_stem"//trim(str(i))//".stl")
            endif
        enddo
    endif
    
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_root"//trim(str(i))//".stl")
            endif
        enddo
    endif

    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                call f%write(trim(name)//"_leaf"//trim(str(i))//".stl")
            endif
        enddo
    endif
    call f%close()
    
    n = input(default=1,option=num_threads)
    !call execute_command_line("echo ' ' > "//trim(name)//".stl")
    !$OMP parallel num_threads(n) private(i)
    !$OMP do 
    do i=1,size(obj%stem)
        if(obj%stem(i)%femdomain%mesh%empty() .eqv. .false. )then
            call obj%stem(i)%stl(name=trim(name)//"_stem"//trim(str(i)))
            !call execute_command_line("cat "//trim(name)//"_stem"//trim(str(i))//"_000001.stl >> "//trim(name)//".stl")
        endif
    enddo
    !$OMP end do
    !$OMP end parallel

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

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

    call execute_command_line("cat "//trim(name)//"*_leaf*.stl > "//trim(name)//"_leaf.stl" )
    call execute_command_line("cat "//trim(name)//"*_stem*.stl > "//trim(name)//"_stem.stl" )
    call execute_command_line("cat "//trim(name)//"*_root*.stl > "//trim(name)//"_root.stl" )
    call execute_command_line("cat "//trim(name)//"_leaf.stl "//trim(name)//"_stem.stl "&
        //trim(name)//"_root.stl > "//trim(name)//".stl" )

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

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

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

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

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

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

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

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

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

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


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

    ! DEBUG
    call f%open("leaf.txt")
    do i=1,size(leafcenter,1)
        write(f%fh,*) leafcenter(i,:)
    enddo
    call f%close()
    
    call f%open("stem.txt")
    do i=1,size(stemcenter,1)
        write(f%fh,*) stemcenter(i,:)
    enddo
    call f%close()
    
    allocate(x(3),x2(3) )
    
    
    num_particle = 0
    totcount = 0
    do i=1,size(obj%leaf)
        if(obj%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then
            ! 葉あり
            obj%leaf(i)%PPFD(:) = max_PPFD
            do j=1,size(obj%leaf(i)%PPFD)
                totcount = totcount + 1
                num_particle = num_particle + 1
                ! それぞれの要素について、遮蔽particleを探索
                ! 茎：全減衰
                ! 葉：半減衰
                ! 簡単のため上からのみ
                ! x-yのみについて見て、上方かつx-y平面距離が半径以内で覆陰判定
                x(:) = leafcenter(num_particle,:)
                r0   = leafradius(num_particle)
                ! 枝による覆陰判定
                
                do k=1, size(stemcenter,1)
                    x2(:) = stemcenter(k,:)
                    r     = stemradius(k)
                    rc    = ( x(1)-x2(1) )**(2.0d0) + ( x(2)-x2(2) )**(2.0d0) 
                    rc    = sqrt(rc)
                    if(rc <= r0 + r .and. x(3) < x2(3) )then
                        ! 茎により覆陰されてる
                        obj%leaf(i)%PPFD(j) = 0.0d0
                        exit
                    endif
                enddo
                if(obj%leaf(i)%PPFD(j) == 0.0d0)then
                    cycle
                endif

                do k=1, size(leafcenter,1)
                    ! もし自信だったら除外
                    if(totcount == k)then
                        cycle
                    endif
                    
                    x2(:) = leafcenter(k,:)
                    r     = leafradius(k)
                    rc    = ( x(1)-x2(1) )**(2.0d0) + ( x(2)-x2(2) )**(2.0d0) 
                    rc    = sqrt(rc)
                    if(rc <= (r0 + r)/2.0d0 .and. x(3) < x2(3) )then
                        ! 茎により覆陰されてる
                        obj%leaf(i)%PPFD(j) = &
                        obj%leaf(i)%PPFD(j)*(1.0d0-extinction_ratio*2.0d0*r)
                        if( obj%leaf(i)%PPFD(j) <= 0.0d0 )then
                            obj%leaf(i)%PPFD(j) = 0.0d0
                        endif
                    endif
                enddo

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


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

subroutine addNodeSoybean(obj,StemNodeID,peti_width_ave,peti_width_sig,peti_size_ave &
    ,peti_size_sig,peti_angle_ave,peti_angle_sig,leaf_thickness_ave,leaf_thickness_sig &
    ,leaf_length_ave,leaf_length_sig,leaf_width_ave,leaf_width_sig,leaf_angle_sig &
    ,leaf_angle_ave)
    class(Soybean_),intent(inout) :: obj
    integer(int32),intent(in) :: StemNodeID
    real(real64),optional,intent(in) :: peti_width_ave,peti_width_sig,peti_size_ave &
    ,peti_size_sig,peti_angle_ave,peti_angle_sig,leaf_thickness_ave,leaf_thickness_sig &
    ,leaf_length_ave,leaf_length_sig,leaf_width_ave,leaf_width_sig,leaf_angle_sig &
    ,leaf_angle_ave
    type(Random_) :: random
    integer(int32) :: i,j

    i = StemNodeID
    obj%num_stem_node = obj%num_stem_node +1
    
    obj%leaf_thickness_ave(obj%num_leaf)=input(&
        default=obj%leaf_thickness_ave(obj%num_leaf),&
        option=leaf_thickness_ave)
    obj%leaf_thickness_sig(obj%num_leaf)=input(&
        default=obj%leaf_thickness_sig(obj%num_leaf),&
        option=leaf_thickness_sig)
    
    obj%leaf_length_ave(obj%num_leaf)=input(&
        default=obj%leaf_length_ave(obj%num_leaf),&
        option=leaf_length_ave)
    obj%leaf_length_sig(obj%num_leaf)=input(&
        default=obj%leaf_length_sig(obj%num_leaf),&
        option=leaf_length_sig)
    obj%leaf_width_ave(obj%num_leaf)=input(&
        default=obj%leaf_width_ave(obj%num_leaf),&
        option=leaf_width_ave)
    obj%leaf_width_sig(obj%num_leaf)=input(&
        default=obj%leaf_width_sig(obj%num_leaf),&
        option=leaf_width_sig)
    obj%leaf_angle_sig(obj%num_leaf)=input(&
        default=obj%leaf_angle_sig(obj%num_leaf),&
        option=leaf_angle_sig)
    
    obj%leaf_angle_ave(obj%num_leaf)=input(&
        default=obj%leaf_angle_ave(obj%num_leaf),&
        option=leaf_angle_ave)

    
    call obj%stem(obj%num_stem_node)%init(config=obj%stemconfig)

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

    

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


! ########################################
subroutine addStemSoybean(obj,stemid,rotx,roty,rotz,json)
    class(Soybean_),intent(inout) :: obj
    integer(int32),intent(in) :: stemid
    character(*),optional,intent(in) :: json
    real(real64),optional,intent(in) :: rotx,roty,rotz
    integer(int32) :: i

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



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

subroutine deformSoybean(obj,penaltyparameter,groundLevel,disp,x_min,x_max,y_min,y_max,z_min,z_max) 
    class(Soybean_),target,intent(inout) :: obj
    real(real64),optional,intent(in) :: groundLevel,disp(3)
    real(real64),optional,intent(in) :: penaltyparameter,x_min,x_max,y_min,y_max,z_min,z_max
    type(FEMDomainp_),allocatable :: domainsp(:)
    integer(int32),allocatable :: contactList(:,:)
    integer(int32) :: i,j,numDomain,stemDomain,leafDomain,rootDomain
    real(real64) :: penalty,displacement(3),GLevel


    if(.not. allocated(obj%Stem) )then
        print *, "ERROR :: deformSoybean >> no soybean is found!"
        return
    endif
    numDomain = 0
    
    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if(obj%stem(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
            endif
        enddo
    endif
    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if(obj%leaf(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
            endif
        enddo
    endif
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if(obj%root(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
            endif
        enddo
    endif
    
    allocate(domainsp(numDomain) )
    numDomain=0
    stemDomain=0
    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if(obj%stem(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                stemDomain = stemDomain + 1
                domainsp(numDomain)%femdomainp =>  obj%stem(i)%femdomain
            endif
        enddo
    endif

    leafDomain = 0
    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if(obj%leaf(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                leafDomain = leafDomain + 1
                domainsp(numDomain)%femdomainp =>  obj%leaf(i)%femdomain
            endif
        enddo
    endif

    rootDomain = 0
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if(obj%root(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                rootDomain = rootDomain + 1
                domainsp(numDomain)%femdomainp =>  obj%root(i)%femdomain
            endif
        enddo
    endif

    contactlist = zeros(numDomain,numDomain)
    if(allocated(obj%stem2stem))then
        do i=1,stemDomain
            do j=1,stemDomain
                contactlist( i, j  ) = obj%stem2stem(i,j)
            enddo
        enddo
    endif

    if(allocated(obj%leaf2stem) )then
        do i=1,leafDomain
            do j=1,stemDomain
                contactlist( i + stemDomain, j  ) = obj%leaf2stem(i,j)
            enddo
        enddo
    endif

    if(allocated(obj%root2stem) )then
        do i=1,rootDomain
            do j=1,stemDomain
                contactlist( i + stemDomain + leafDomain, j  ) = obj%root2stem(i,j)
            enddo
        enddo
    endif

    if(allocated(obj%root2root) )then
        do i=1,rootDomain
            do j=1,rootDomain
                contactlist( i + stemDomain + leafDomain, j+ stemDomain + leafDomain  ) = obj%root2root(i,j)
            enddo
        enddo
    endif
    !call print(contactlist)
    !stop
    call obj%contact%init(femdomainsp=domainsp,contactlist=contactlist)

    ! load material info
    numDomain = 0
    if(allocated(obj%stem) )then
        do i=1,size(obj%stem)
            if(obj%stem(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                call obj%contact%setYoungModulus(YoungModulus=obj%stemYoungModulus(i),DomainID=numDomain) 
                call obj%contact%setPoissonRatio(PoissonRatio=obj%stemPoissonRatio(i),DomainID=numDomain) 
                call obj%contact%setDensity(density=obj%stemDensity(i),DomainID=numDomain) 
            endif
        enddo
    endif
    if(allocated(obj%leaf) )then
        do i=1,size(obj%leaf)
            if(obj%leaf(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                call obj%contact%setYoungModulus(YoungModulus=obj%leafYoungModulus(i),DomainID=numDomain) 
                call obj%contact%setPoissonRatio(PoissonRatio=obj%leafPoissonRatio(i),DomainID=numDomain) 
                call obj%contact%setDensity(density=obj%leafDensity(i),DomainID=numDomain) 
            endif
        enddo
    endif
    if(allocated(obj%root) )then
        do i=1,size(obj%root)
            if(obj%root(i)%femdomain%mesh%empty() )then
                cycle
            else
                numDomain = numDomain + 1
                call obj%contact%setYoungModulus(YoungModulus=obj%rootYoungModulus(i),DomainID=numDomain) 
                call obj%contact%setPoissonRatio(PoissonRatio=obj%rootPoissonRatio(i),DomainID=numDomain) 
                call obj%contact%setDensity(density=obj%rootDensity(i),DomainID=numDomain) 
            endif
        enddo
    endif
    !


    penalty = input(default=1000.0d0, option=penaltyparameter)
    
    call obj%contact%setup(penaltyparameter=penalty)

    ! if displacement is set, load displacement
    if(present(disp) )then
        do i=1,numDomain
            call obj%contact%fix(direction="x",disp=disp(1), DomainID=i,&
                x_min=x_min,x_max=x_max,&
                y_min=y_min,y_max=y_max,&
                z_min=z_min,z_max=z_max)
            call obj%contact%fix(direction="y",disp=disp(2), DomainID=i,&
                x_min=x_min,x_max=x_max,&
                y_min=y_min,y_max=y_max,&
                z_min=z_min,z_max=z_max)
            call obj%contact%fix(direction="z",disp=disp(3), DomainID=i,&
                x_min=x_min,x_max=x_max,&
                y_min=y_min,y_max=y_max,&
                z_min=z_min,z_max=z_max)
        enddo    
    endif
    

    Glevel = input(default=0.0d0,option=groundLevel)
    ! under-ground parts are fixed.
    do i=1,numDomain
        call obj%contact%fix(direction="x",disp=0.0d0, DomainID=i,&
            z_max=Glevel)
        call obj%contact%fix(direction="y",disp=0.0d0, DomainID=i,&
            z_max=Glevel)
        call obj%contact%fix(direction="z",disp=0.0d0, DomainID=i,&
            z_max=Glevel)
    enddo

    ! solve > get displacement
    call obj%contact%solver%solve("BiCGSTAB")
    ! update mesh
    call obj%contact%updateMesh()


end subroutine

function getVolumeSoybean(obj,stem,leaf,root) result(ret)
    class(Soybean_),intent(in) :: obj
    logical,optional,intent(in) :: stem, leaf, root
    logical :: all
    integer(int32) :: i,j
    real(real64) :: ret

    all = .false.
    if(.not.present(stem) .and..not.present(leaf)  )then
        if(.not. present(root) )then
            all = .true.
        endif
    endif

    ret =0.0d0
    if(all)then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                do j=1,obj%stem(i)%femdomain%ne()
                    ret = ret + obj%stem(i)%femdomain%getVolume(elem=j)
                enddo
            endif
        enddo
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                do j=1,obj%leaf(i)%femdomain%ne()
                    ret = ret + obj%leaf(i)%femdomain%getVolume(elem=j)
                enddo
            endif
        enddo
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                do j=1,obj%root(i)%femdomain%ne()
                    ret = ret + obj%root(i)%femdomain%getVolume(elem=j)
                enddo
            endif
        enddo
        return
    endif

    if(present(stem))then
        if(stem  .or. all)then
            do i=1,size(obj%stem)
                if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                    do j=1,obj%stem(i)%femdomain%ne()
                        ret = ret + obj%stem(i)%femdomain%getVolume(elem=j)
                    enddo
                endif
            enddo
        endif
    endif
    if(present(leaf) )then
        if(leaf )then
            do i=1,size(obj%leaf)
                if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                    do j=1,obj%leaf(i)%femdomain%ne()
                        ret = ret + obj%leaf(i)%femdomain%getVolume(elem=j)
                    enddo
                endif
            enddo
        endif
    endif
    if(present(root))then
        if(root)then
            do i=1,size(obj%root)
                if( .not.obj%root(i)%femdomain%mesh%empty() )then
                    do j=1,obj%root(i)%femdomain%ne()
                        ret = ret + obj%root(i)%femdomain%getVolume(elem=j)
                    enddo
                endif
            enddo
        endif
    endif

end function

function getBiomassSoybean(obj,stem,leaf,root) result(ret)
    class(Soybean_),intent(in) :: obj
    logical,optional,intent(in) :: stem, leaf, root
    logical :: all
    integer(int32) :: i,j
    real(real64) :: ret,volume

    all = .false.
    if(.not.present(stem) .and..not.present(leaf)  )then
        if(.not. present(root) )then
            all = .true.
        endif
    endif

    ret =0.0d0
    if(all)then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                do j=1,obj%stem(i)%femdomain%ne()
                    volume = obj%stem(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume)
                    ret = ret + volume*obj%stem(i)%drydensity(j) 
                enddo
                
            endif
        enddo
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                do j=1,obj%leaf(i)%femdomain%ne()
                    volume = obj%leaf(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume) 
                    ret = ret + volume*obj%leaf(i)%drydensity(j) 
                enddo
                
            endif
        enddo
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                do j=1,obj%root(i)%femdomain%ne()
                    volume = obj%root(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume) 
                    ret = ret + volume*obj%root(i)%drydensity(j) 
                enddo
                
            endif
        enddo
        return
    endif

    if(present(stem))then
        if(stem  .or. all)then
            do i=1,size(obj%stem)
                if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                    do j=1,obj%stem(i)%femdomain%ne()
                        volume = obj%stem(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) 
                        ret = ret + volume*obj%stem(i)%drydensity(j) 
                    enddo 
                endif
            enddo
        endif
    endif
    if(present(leaf) )then
        if(leaf )then
            do i=1,size(obj%leaf)
                if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                    do j=1,obj%leaf(i)%femdomain%ne()
                        volume = obj%leaf(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) 
                        ret = ret + volume*obj%leaf(i)%drydensity(j) 
                    enddo
                endif
            enddo
        endif
    endif
    if(present(root))then
        if(root)then
            do i=1,size(obj%root)
                if( .not.obj%root(i)%femdomain%mesh%empty() )then
                    do j=1,obj%root(i)%femdomain%ne()
                        volume = obj%root(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) 
                        ret = ret + volume*obj%root(i)%drydensity(j) 
                    enddo
                endif
            enddo
        endif
    endif

end function


function getTotalWeightSoybean(obj,stem,leaf,root,waterDensity) result(ret)
    class(Soybean_),intent(in) :: obj
    logical,optional,intent(in) :: stem, leaf, root
    real(real64),optional,intent(in) :: waterDensity
    logical :: all
    integer(int32) :: i,j
    real(real64) :: ret,volume,water_density

    ! kg, m
    water_density=input(default=1000.0d0,option=waterDensity)

    all = .false.
    if(.not.present(stem) .and..not.present(leaf)  )then
        if(.not. present(root) )then
            all = .true.
        endif
    endif

    ret =0.0d0
    if(all)then
        do i=1,size(obj%stem)
            if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                do j=1,obj%stem(i)%femdomain%ne()
                    volume = obj%stem(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                    ret = ret + volume*obj%stem(i)%drydensity(j) + volume*obj%stem(i)%watercontent(j)*water_density
                enddo
                
            endif
        enddo
        do i=1,size(obj%leaf)
            if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                do j=1,obj%leaf(i)%femdomain%ne()
                    volume = obj%leaf(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                    ret = ret + volume*obj%leaf(i)%drydensity(j) + volume*obj%leaf(i)%watercontent(j)*water_density
                enddo
                
            endif
        enddo
        do i=1,size(obj%root)
            if( .not.obj%root(i)%femdomain%mesh%empty() )then
                do j=1,obj%root(i)%femdomain%ne()
                    volume = obj%root(i)%femdomain%getVolume(elem=j)
                    ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                    ret = ret + volume*obj%root(i)%drydensity(j) + volume*obj%root(i)%watercontent(j)*water_density
                enddo
                
            endif
        enddo
        return
    endif

    if(present(stem))then
        if(stem  .or. all)then
            do i=1,size(obj%stem)
                if( .not.obj%stem(i)%femdomain%mesh%empty() )then
                    do j=1,obj%stem(i)%femdomain%ne()
                        volume = obj%stem(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                        ret = ret + volume*obj%stem(i)%drydensity(j) + volume*obj%stem(i)%watercontent(j)*water_density
                    enddo 
                endif
            enddo
        endif
    endif
    if(present(leaf) )then
        if(leaf )then
            do i=1,size(obj%leaf)
                if( .not.obj%leaf(i)%femdomain%mesh%empty() )then
                    do j=1,obj%leaf(i)%femdomain%ne()
                        volume = obj%leaf(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                        ret = ret + volume*obj%leaf(i)%drydensity(j) + volume*obj%leaf(i)%watercontent(j)*water_density
                    enddo
                endif
            enddo
        endif
    endif
    if(present(root))then
        if(root)then
            do i=1,size(obj%root)
                if( .not.obj%root(i)%femdomain%mesh%empty() )then
                    do j=1,obj%root(i)%femdomain%ne()
                        volume = obj%root(i)%femdomain%getVolume(elem=j)
                        ! total = total + solid(=drydensity * volume) + fluid (=watercontent * volume)
                        ret = ret + volume*obj%root(i)%drydensity(j) + volume*obj%root(i)%watercontent(j)*water_density
                    enddo
                endif
            enddo
        endif
    endif

end function


!function getBioMassSoybean(obj,stemDensity,leafDensity,rootDensity) result(ret)
!    class(Soybean_),intent(in) :: obj
!    real(real64),optional,intent(in) :: stemDensity,leafDensity,rootDensity
!    logical :: all
!    integer(int32) :: i,j
!    real(real64) :: ret
!
!    ret = 0.0d0
!
!    if(present(stemDensity))then
!        ret = ret + obj%getVolume(stem=.true.) * stemDensity
!    endif
!
!    if(present(leafDensity))then
!        ret = ret + obj%getVolume(leaf=.true.) * leafDensity
!    endif
!
!    if(present(rootDensity))then
!        ret = ret + obj%getVolume(root=.true.) * rootDensity
!    endif
!
!
!
!end function

subroutine removeSoybean(obj)
    class(Soybean_),intent(inout) :: obj

    obj%growth_habit = " "
    obj%growth_stage = " "
    obj%Num_Of_Node = 0
    obj%Num_Of_Root = 0
    
    obj%MaxLeafNum= 300
    obj%MaxRootNum=300
    obj%MaxStemNum= 300

    
    
    obj%ms_node=0
    obj%br_node(:) = 0
    obj%br_from(:) = 0
    obj%ms_length=0.0d0
    obj%br_length(:)=0.0d0
    obj%ms_width=0.0d0
    obj%br_width(:)=0.0d0
    obj%ms_angle_ave=0.0d0
    obj%br_angle_ave(:)=0.0d0
    obj%ms_angle_sig=0.0d0
    obj%br_angle_sig(:)=0.0d0
    

    obj%mr_node=0
    obj%brr_node(:)=0
    obj%brr_from(:)=0
    obj%mr_length=0.0d0
    obj%brr_length(:)=0.0d0
    obj%mr_width=0.0d0
    obj%brr_width(:)=0.0d0
    obj%mr_angle_ave=0.0d0
    obj%brr_angle_ave(:)=0.0d0
    obj%mr_angle_sig=0.0d0
    obj%brr_angle_sig(:)=0.0d0

    obj%peti_size_ave(:)=0.0d0
    obj%peti_size_sig(:)=0.0d0
    obj%peti_width_ave(:)=0.0d0
    obj%peti_width_sig(:)=0.0d0
    obj%peti_angle_ave(:)=0.0d0
    obj%peti_angle_sig(:)=0.0d0

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

    if (allocated(obj%Stem) ) deallocate(obj%Stem)
    if (allocated(obj%Leaf) ) deallocate(obj%Leaf)
    if (allocated(obj%Root) ) deallocate(obj%Root)
    if (allocated(obj%Soil) ) deallocate(obj%Soil)

    ! material info
    if (allocated(obj%stemYoungModulus) ) deallocate(obj%stemYoungModulus)
    if (allocated(obj%leafYoungModulus) ) deallocate(obj%leafYoungModulus)
    if (allocated(obj%rootYoungModulus) ) deallocate(obj%rootYoungModulus)

    if (allocated(obj%stemPoissonRatio) ) deallocate(obj%stemPoissonRatio)
    if (allocated(obj%leafPoissonRatio) ) deallocate(obj%leafPoissonRatio)
    if (allocated(obj%rootPoissonRatio) ) deallocate(obj%rootPoissonRatio)

    if (allocated(obj%stemDensity) ) deallocate(obj%stemDensity)
    if (allocated(obj%leafDensity) ) deallocate(obj%leafDensity)
    if (allocated(obj%rootDensity) ) deallocate(obj%rootDensity)
    

    if(allocated(obj%NodeID_MainStem)) deallocate(obj%NodeID_MainStem)
    if(allocated(obj%NodeID_Branch)) deallocate(obj%NodeID_Branch)
    ! 節-節点データ構造
    call obj%struct%remove(all=.true.)
    if (allocated(obj%leaf2stem) ) deallocate(obj%leaf2stem)
    if (allocated(obj%stem2stem) ) deallocate(obj%stem2stem)
    if (allocated(obj%root2stem) ) deallocate(obj%root2stem)
    if (allocated(obj%root2root) ) deallocate(obj%root2root)
    
    ! 器官オブジェクト配列
    if (allocated(obj%leaf_list) ) deallocate(obj%leaf_list)
    if (allocated(obj%stem_list) ) deallocate(obj%stem_list)
    if (allocated(obj%root_list) ) deallocate(obj%root_list)

    ! シミュレータ
    call obj%contact%remove()
    obj%time=0.0d0
    obj%seed_length=0.0d0
    obj%seed_width=0.0d0
    obj%seed_height=0.0d0
    if (allocated(obj%stem_angle) ) deallocate(obj%stem_angle)
    if (allocated(obj%root_angle) ) deallocate(obj%root_angle)
    if (allocated(obj%leaf_angle) ) deallocate(obj%leaf_angle)

    obj%stemconfig=" "
    obj%rootconfig=" "
    obj%leafconfig=" "

end subroutine

function stemlengthSoybean(obj,StemID) result(ret)
    class(Soybean_),intent(inout) :: obj
    integer(int32),intent(in) :: StemID ! 0, 1, 2...
    integer(int32) :: num_snode,i
    real(real64),allocatable :: ret(:)

    if(StemID==0)then
        ! main stem
        num_snode = size(obj%NodeID_MainStem)
        allocate(ret(num_snode) )
        do i=1,num_snode
            ret(i) = obj%stem(obj%NodeID_MainStem(i))%getLength()
        enddo
    else
        if(StemID >=size(obj%NodeID_Branch)) then
            print *, "ERROR :: stemlengthSoybean >> StemID >=size(obj%NodeID_Branch)"
            ret = zeros(1)
            return
        endif
        ! main stem
        num_snode = size(obj%NodeID_Branch(StemID)%ID)
        allocate(ret(num_snode) )
        do i=1,num_snode
            ret(i) = obj%stem(obj%NodeID_Branch(StemID)%ID(i))%getLength()
        enddo
    endif



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

! ###################################################################
subroutine resizeSoybean(obj,StemID,StemLength)
    class(Soybean_),intent(inout) :: obj
    integer(int32),optional,intent(in) :: StemID
    real(real64),optional,intent(in) :: StemLength(:)
    integer(int32) :: num_snode,i

    if(present(StemID) )then
        if(.not.present(StemLength) )then
            print *, "ERROR :: resizeSoybean >> needs StemLength(:) "
            stop
        endif

        if(StemID==0)then
            ! main stem
            num_snode = size(obj%NodeID_MainStem)
            
            do i=1,num_snode
                call obj%stem(obj%NodeID_MainStem(i))%grow(length=StemLength(i))
            enddo
        else
            if(StemID >=size(obj%NodeID_Branch)) then
                print *, "ERROR :: resizeSoybean >> StemID >=size(obj%NodeID_Branch)"
                
                return
            endif
            ! main stem
            num_snode = size(obj%NodeID_Branch(StemID)%ID)
            
            do i=1,num_snode
                call obj%stem(obj%NodeID_Branch(StemID)%ID(i))%grow(length=StemLength(i) )
            enddo
        endif
        call obj%update()
    

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


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

    ret = 0
    do i=1,size(obj%br_length)
        if(obj%br_length(i)/=0.0d0 )then
            ret = ret + 1
        endif
    enddo

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

! ###################################################################
function findApicalSoybean(obj) result(ret)
    class(Soybean_),intent(in) :: obj
    integer(int32),allocatable :: ret(:)
    !integer(int32),optional,intent(in) :: StemID
    integer(int32),allocatable :: stem
    integer(int32) :: i,j

    ret = zeros( obj%NumberOfBranch()+1 )

    ret(1)=maxval(obj%NodeID_MainStem(:) )
    do i=1,obj%NumberOfBranch()
        ret(i+1) = maxval(obj%NodeID_Branch(i)%ID(:) )
    enddo
    
!    if(present(StemID) )then
!        if(StemID > size(obj%br_node) )then
!            print *, "ERROR >> findApicalSoybean >> number of branch is ",size(obj%br_node)
!            print *, "StemID=",StemID,"is larger than it."
!            return
!        endif
!
!
!        return
!    endif

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

!function propertiesSoybean(obj) result(ret)
!    class(Soybean_) ,intent(in) :: obj
!
!    
!end function

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

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

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



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

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

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


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

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

end function


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

function getSubDomainSoybean(obj,id) result(ret)
    class(Soybean_),intent(in) :: obj
    integer(int32),intent(in) :: id
    type(FEMDomain_) :: ret
    integer(int32) :: i,ret_id
    
    ! get number of subdomain
    ret_id = 0
    do i=1,size(obj%stem)
        if( .not.obj%stem(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = obj%stem(i)%femdomain
                return
            endif
        endif
    enddo
    do i=1,size(obj%leaf)
        if( .not.obj%leaf(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = obj%stem(i)%femdomain
                return
            endif
        endif
    enddo
    do i=1,size(obj%root)
        if( .not.obj%root(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = obj%stem(i)%femdomain
                return
            endif
        endif
    enddo

    print *, "Caution >> getSubDomainSoybean >> exceed total number of subdomains",ret_id
    return

    

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

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

subroutine setSubDomainSoybean(obj,domain,id) 
    class(Soybean_),intent(inout) :: obj
    type(FEMDomain_),intent(in) :: domain
    integer(int32),intent(in) :: id
    integer(int32) :: i,domain_id
    
    ! get number of subdomain
    domain_id = 0
    do i=1,size(obj%stem)
        if( .not.obj%stem(i)%femdomain%mesh%empty() ) then
            domain_id = domain_id + 1
            if(id==domain_id)then
                obj%stem(i)%femdomain = domain
                return
            endif
        endif
    enddo
    do i=1,size(obj%leaf)
        if( .not.obj%leaf(i)%femdomain%mesh%empty() ) then
            domain_id = domain_id + 1
            if(id==domain_id)then
                obj%stem(i)%femdomain = domain
                return
            endif
        endif
    enddo
    do i=1,size(obj%root)
        if( .not.obj%root(i)%femdomain%mesh%empty() ) then
            domain_id = domain_id + 1
            if(id==domain_id)then
                obj%stem(i)%femdomain = domain
                return
            endif
        endif
    enddo

    print *, "Caution >> getSubDomainSoybean >> exceed total number of subdomains",domain_id
    return

    

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


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

function getSubDomainTypeSoybean(obj,id) result(ret)
    class(Soybean_),intent(in) :: obj
    integer(int32),intent(in) :: id
    character(:),allocatable :: ret
    integer(int32) :: i,ret_id
    
    ! get number of subdomain
    ret_id = 0
    do i=1,size(obj%stem)
        if( .not.obj%stem(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = "stem"
                return
            endif
        endif
    enddo
    do i=1,size(obj%leaf)
        if( .not.obj%leaf(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = "leaf"
                return
            endif
        endif
    enddo
    do i=1,size(obj%root)
        if( .not.obj%root(i)%femdomain%mesh%empty() ) then
            ret_id = ret_id + 1
            if(id==ret_id)then
                ret = "root"
                return
            endif
        endif
    enddo

    print *, "Caution >> getSubDomainSoybean >> exceed total number of subdomains",ret_id
    return

    

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