module LeafClass
   use, intrinsic :: iso_fortran_env
   use KinematicClass
   use FEMDomainClass
   use PetiClass
   use StemClass
   use LightClass
   use AirClass
   implicit none

   integer(int32) :: PF_SOYBEAN_CV = 100
   type :: Leaf_
      type(FEMDomain_)    ::  FEMDomain
      real(real64), allocatable ::  LeafSurfaceNode2D(:, :)
      real(real64)             ::  ShapeFactor, Thickness, length, width, center(3)
      real(real64)             ::  MaxThickness, Maxlength, Maxwidth
      real(real64)             ::  center_bottom(3), center_top(3)
      real(real64)             ::  outer_normal_bottom(3), outer_normal_top(3)
      real(real64), allocatable ::  source(:), ppfd(:), A(:)
      integer(int32)             ::  Division
      type(leaf_), pointer ::  pleaf
      type(Peti_), pointer ::  pPeti
      real(real64)             ::  rot_x = 0.0d0
      real(real64)             ::  rot_y = 0.0d0
      real(real64)             ::  rot_z = 0.0d0
      real(real64)             ::  disp_x = 0.0d0
      real(real64)             ::  disp_y = 0.0d0
      real(real64)             ::  disp_z = 0.0d0
      real(real64)             ::  shaperatio = 0.30d0
      real(real64)             ::  minwidth, minlength, MinThickness

      ! id in multi-leaf
      integer(int32) :: LeafID = -1
      logical :: already_grown = .false.

      integer(int32), allocatable  :: I_planeNodeID(:)
      integer(int32), allocatable  :: I_planeElementID(:)
      integer(int32), allocatable  :: II_planeNodeID(:)
      integer(int32), allocatable  :: II_planeElementID(:)
      integer(int32)  :: A_PointNodeID
      integer(int32)  :: B_PointNodeID
      integer(int32)  :: C_PointNodeID
      integer(int32)  :: D_PointNodeID
      integer(int32)  :: A_PointElementID
      integer(int32)  :: B_PointElementID
      integer(int32)  :: C_PointElementID
      integer(int32)  :: D_PointElementID
      integer(int32)  :: xnum = 10
      integer(int32)  :: ynum = 10
      integer(int32)  :: znum = 10

      ! phisiological parameters

      real(real64) :: V_cmax = 100.0d0 ! 最大カルボキシル化反応速度, mincro-mol/m-2/s
      real(real64) :: V_omax = 100.0d0 ! 最大酸素化反応速度, mincro-mol/m-2/s, lambdaから推定
      real(real64) :: CO2 = 380.0d0! 二酸化炭素濃度, ppm
      real(real64) :: O2 = 202000.0d0! 酸素濃度, ppm
      real(real64) :: R_d = 1.0d0 ! 暗呼吸速度, mincro-mol/m-2/s

!        real(real64) :: K_c=272.380d0 ! CO2に対するミカエリス定数
!        real(real64) :: K_o=165820.0d0 ! O2に対するミカエリス定数

      real(real64) :: K_c = 272.380d0  ! CO2に対するミカエリス定数
      real(real64) :: K_o = 165820.0d0 ! O2に対するミカエリス定数

      real(real64) :: J_ = 0.0d0 ! 電子伝達速度
      real(real64) :: I_ = 0.0d0 ! 光強度

      real(real64) :: phi = 0.0d0 ! I-J曲線の初期勾配
      real(real64) :: J_max = 180.0d0 !最大電子伝達速度,mincro-mol/m-2/s
      real(real64) :: theta_r = 0.0d0 ! 曲線の凸度

      real(real64) :: maxPPFD = 1000.0d0 ! micro-mol/m^2/s

      real(real64) :: Lambda = 37.430d0 ! 暗呼吸速度を無視した時のCO2補償点ppm
      real(real64) :: temp = 303.0d0 ! temp

      ! physical parameter
      real(real64), allocatable :: DryDensity(:)  ! element-wise
      real(real64), allocatable :: WaterContent(:)! element-wise

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

      real(real64), allocatable :: BoundaryTractionForce(:, :) ! node-wise, three dimensional
      real(real64), allocatable :: BoundaryDisplacement(:, :) ! node-wise, three dimensional

      ! growth parameters
      real(real64)  :: my_time = 0.0d0
      real(real64)  :: initial_width = 0.0010d0 ! 1.0 mm
      real(real64)  :: initial_length = 0.0010d0 ! 1.0 mm
      real(real64)  :: final_width = 0.060d0   !  60.0 mm
      real(real64)  :: final_length = 0.100d0   ! 100.0 mm
      real(real64)  :: width_growth_ratio = 1.0d0/4.0d0   !
      real(real64)  :: length_growth_ratio = 1.0d0/4.0d0   !

      !Sink-source flow parameter
      real(real64) :: default_CarbonDiffusionCoefficient = 0.0010d0 ! ソースの拡散係数 mincro-mol/m^2/m/s

   contains
      procedure, public :: Init => initLeaf
      procedure, public :: change_length_or_width => change_length_or_width_Leaf

      procedure, public :: rotate => rotateleaf
      procedure, public :: move => moveleaf
      procedure, public :: curve => curveleaf
      procedure, public :: create => createLeaf

      procedure, pass :: connectLeafLeaf => connectLeafLeaf
      procedure, pass :: connectLeafStem => connectLeafStem

      generic :: connect => connectLeafLeaf, connectLeafStem

      procedure, public :: photosynthesis => photosynthesisLeaf
      procedure, public :: getPhotosynthesisSpeedPerVolume => getPhotosynthesisSpeedPerVolumeLeaf

      procedure, public :: rescale => rescaleleaf
      procedure, public :: adjust => adjustLeaf
      procedure, public :: resize => resizeleaf

      ! check condition
      ! is it empty?
      procedure, public :: empty => emptyLeaf
      procedure, public :: getVolume => getVolumeLeaf
      procedure, public :: getLength => getLengthLeaf
      procedure, public :: getBiomass => getBiomassLeaf
      procedure, public :: getCoordinate => getCoordinateleaf
      procedure, public :: getLeafArea => getLeafAreaLeaf
      procedure, public :: getRadius => getRadiusLeaf
      procedure, public :: getCenter => getCenterLeaf
      procedure, public :: getNormalVector => getNormalVectorLeaf

      procedure, public :: FullyExpanded => FullyExpandedLeaf

      procedure, public :: gmsh => gmshleaf
      procedure, public :: msh => mshleaf
      procedure, public :: vtk => vtkleaf
      procedure, public :: stl => stlleaf
      procedure, public :: ply => plyleaf

      procedure, public :: sync => syncleaf

      procedure, public :: nn => nnLeaf
      procedure, public :: ne => neLeaf

      ! remove
      procedure, public :: remove => removeLeaf
   end type

   interface operator(//)
      module procedure append_leaf_object_vector
   end interface
contains

   subroutine createLeaf(obj, SurfacePoints, filename, x_num, y_num, x_len, y_len)
      class(Leaf_), intent(inout) :: obj
      real(real64), optional, intent(in) :: SurfacePoints(:, :), x_len, y_len
      character(*), optional, intent(in) :: filename
      integer(int32), optional, intent(in) :: x_num, y_num

      type(IO_) :: f
      type(FEMDomain_) :: domain
      type(Math_) :: math
      character(:), allocatable :: line
      real(real64) :: x, y, r, theta, x_sum, y_sum, center(2), max_r, coord(2), ret
      real(real64), allocatable :: r_data(:), theta_data(:), tx(:), tfx(:)
      integer(int32) :: num_ptr, i, id, ids(5), id_n

      if (present(filename)) then
         call f%open(filename, "r")
         ! get brief info
         num_ptr = 0

         x_sum = 0.0d0
         y_sum = 0.0d0

         do
            line = f%readline()
            if (f%EOF) exit
            num_ptr = num_ptr + 1
            ! read x-y
            read (line, *) x, y
            x_sum = x_sum + x
            y_sum = y_sum + y
         end do
         call f%close()

         center(1) = x_sum/dble(num_ptr)
         center(2) = y_sum/dble(num_ptr)

         r_data = zeros(num_ptr)
         theta_data = zeros(num_ptr)

         ! get detail
         call f%open(filename, "r")
         num_ptr = 0
         do
            line = f%readline()
            if (f%EOF) exit
            ! read x-y
            read (line, *) x, y

            coord(1) = x - center(1)
            coord(2) = y - center(2)
            r = sqrt(dot_product(coord, coord))
            theta = angles(coord)

            num_ptr = num_ptr + 1

            r_data(num_ptr) = r
            theta_data(num_ptr) = theta
         end do
         max_r = maxval(r_data)
         r_data = r_data/max_r
         call f%close()
      elseif (present(SurfacePoints)) then
         num_ptr = size(SurfacePoints, 1)
         center(1) = x_sum/dble(num_ptr)
         center(2) = y_sum/dble(num_ptr)

         r_data = zeros(num_ptr)
         theta_data = zeros(num_ptr)

         num_ptr = 0
         do i = 1, size(SurfacePoints)

            ! read x-y
            x = SurfacePoints(i, 1)
            y = SurfacePoints(i, 2)

            coord(1) = x - center(1)
            coord(2) = y - center(2)

            r = sqrt(dot_product(coord, coord))
            theta = angles(coord)

            num_ptr = num_ptr + 1

            r_data(num_ptr) = r
            theta_data(num_ptr) = theta
         end do
         max_r = maxval(r_data)
         r_data = r_data/max_r

      else
         print *, "ERROR :: Leaf%create >> Please import SurfacePoints or Filename"
         stop
      end if

      call obj%femdomain%create("Cylinder3D", x_num=x_num, y_num=y_num)
      call obj%femdomain%resize(x=2.0d0)
      call obj%femdomain%resize(y=2.0d0)
      call obj%femdomain%resize(z=0.010d0)

      ! ####################################
      ! test interpolate

      !tx = [0.0d0, 1.0d0, 2.0d0, 3.0d0]
      !tfx = [0.0d0, 2.0d0, 4.0d0, 8.0d0]
      !ret = interpolate(x =tx,Fx=tfx,x_value = -0.50d0)
      !print *, ret
      !stop
      ! ####################################

      ! adjust shape
      do i = 1, obj%femdomain%nn()
         x = obj%femdomain%mesh%nodcoord(i, 1)
         y = obj%femdomain%mesh%nodcoord(i, 2)
         r = sqrt(x**2 + y**2)
         coord(1:2) = obj%femdomain%mesh%nodcoord(i, 1:2)
         r = norm(coord)
         theta = angles(coord)
         ! find nearest theta
         r = r*interpolate(x=theta_data, Fx=r_data, x_value=theta)
         x = r*x
         y = r*y
         obj%femdomain%mesh%nodcoord(i, 1) = x
         obj%femdomain%mesh%nodcoord(i, 2) = y
      end do

      obj%A_PointNodeID = randi(obj%femdomain%nn())
      obj%B_PointNodeID = randi(obj%femdomain%nn())
      obj%A_PointElementID = randi(obj%femdomain%nn())
      obj%B_PointElementID = randi(obj%femdomain%nn())

      if (present(x_len)) then
         call obj%femdomain%resize(x=x_len)
      end if

      if (present(y_len)) then
         call obj%femdomain%resize(y=y_len)
      end if

      obj%thickness = maxval(obj%femdomain%mesh%nodcoord(:, 3)) &
                      - minval(obj%femdomain%mesh%nodcoord(:, 3))
!    ! export data
!    call f%open("theta_r_relation.txt","w")
!    do i=1,size(r_data)
!        call f%write(theta_data(i),r_data(i) )
!    enddo
!    call f%close()
!    call f%plot("theta_r_relation.txt","w l")
!    call f%plot(filename,"w l")

   end subroutine

! ########################################
   subroutine initLeaf(obj, config, regacy, Thickness, length, width, ShapeFactor, &
                       MaxThickness, Maxlength, Maxwidth, rotx, roty, rotz, location, species, SoyWidthRatio, &
                       curvature, x_num, y_num, z_num)
      class(leaf_), intent(inout) :: obj
      real(real64), optional, intent(in) :: Thickness, length, width, ShapeFactor
      real(real64), optional, intent(in) :: MaxThickness, Maxlength, Maxwidth
      real(real64), optional, intent(in)::  rotx, roty, rotz, location(3), SoyWidthRatio, curvature
      integer(int32), optional, intent(in) :: species, x_num, y_num, z_num
      logical, optional, intent(in) :: regacy
      character(*), optional, intent(in) :: config
      type(IO_) :: leafconf, f
      character(200) :: fn, conf, line
      integer(int32), allocatable :: buf(:)
      integer(int32) :: id, rmc, n, node_id, node_id2, elemid, blcount, i, j

      integer(int32), allocatable :: killElementList(:), killElementIDs(:)
      real(real64) :: loc(3), radius, z, leaf_L
      logical :: debug = .false.

      obj%minlength = 0.005d0
      obj%minwidth = 0.005d0
      obj%minthickness = 0.0001d0
      obj%maxlength = 0.07d0
      obj%maxwidth = 0.045d0
      obj%maxthickness = 0.001d0
      obj%shaperatio = 0.3d0
      !obj%drydensity =  0.0
      !obj%watercontent =  0.0
      obj%xnum = 10
      obj%ynum = 10
      obj%znum = 20

!        ! 節を生成するためのスクリプトを開く
!        if(.not.present(config) .or. index(config,".json")==0 )then
!            ! デフォルトの設定を生成
!            if(debug) print *, "New leaf-configuration >> leafconfig.json"
!            call leafconf%open("leafconfig.json")
!            write(leafconf%fh,*) '{'
!            write(leafconf%fh,*) '   "type": "leaf",'
!            write(leafconf%fh,*) '   "minlength": 0.005,'
!            write(leafconf%fh,*) '   "minwidth": 0.005,'
!            write(leafconf%fh,*) '   "minthickness": 0.0001,'
!            write(leafconf%fh,*) '   "maxlength": 0.07,'
!            write(leafconf%fh,*) '   "maxwidth": 0.045,'
!            write(leafconf%fh,*) '   "maxthickness": 0.001,'
!            write(leafconf%fh,*) '   "shaperatio": 0.3,'
!            write(leafconf%fh,*) '   "drydensity": 0.0,'
!            write(leafconf%fh,*) '   "watercontent": 0.0,'
!            write(leafconf%fh,*) '   "xnum": 10,'
!            write(leafconf%fh,*) '   "ynum": 10,'
!            write(leafconf%fh,*) '   "znum": 20'
!            write(leafconf%fh,*) '}'
!            conf="leafconfig.json"
!            call leafconf%close()
!        else
!            conf = config
!        endif

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

            if (blcount == 1) then

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

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

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

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

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

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

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

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

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

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

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

            end if

         end do
         call leafconf%close()
      end if

      obj%xnum = input(default=obj%xnum, option=x_num)
      obj%ynum = input(default=obj%ynum, option=y_num)
      obj%znum = input(default=obj%znum, option=z_num)

      ! グラフ構造とメッシュ構造を生成する。

      !
      !           D %%%%%%%%%%%%%%%%%%%%%%%%%%%  B
      !         %%                        %   %
      !        %%                    %      %%
      !      %%                 %          %%
      !     %%            %              %%
      !     %%      %                  %%
      !     %%                       %%
      !   A   %%                  %%
      !      <I> %%%%%%%%%%%%%%%% C

      ! メッシュを生成
      call obj%FEMdomain%create(meshtype="rectangular3D", x_num=obj%xnum, y_num=obj%ynum, z_num=obj%znum, &
                                x_len=obj%minwidth/2.0d0, y_len=obj%minwidth/2.0d0, z_len=obj%minlength, shaperatio=obj%shaperatio)

      ! physical parameters
      allocate (obj%A(size(obj%FEMDomain%Mesh%ElemNod, 1)))
      obj%A(:) = 0.0d0
      allocate (obj%source(size(obj%FEMDomain%Mesh%ElemNod, 1)))
      obj%source(:) = 0.0d0
      allocate (obj%ppfd(size(obj%FEMDomain%Mesh%ElemNod, 1)))
      obj%ppfd(:) = 0.0d0

      ! initialize physical parameter
      obj%DryDensity = zeros(obj%FEMDomain%ne())
      obj%watercontent = zeros(obj%FEMDomain%ne())
      if (present(config)) then
         obj%DryDensity(:) = freal(leafconf%parse(config, key1="drydensity"))
         obj%watercontent(:) = freal(leafconf%parse(config, key1="watercontent"))
      end if
      ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
      obj%I_planeNodeID = obj%FEMdomain%mesh%getNodeList(zmax=0.0d0)
      obj%I_planeElementID = obj%FEMdomain%mesh%getElementList(zmax=0.0d0)

      ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
      obj%II_planeNodeID = obj%FEMdomain%mesh%getNodeList(zmin=obj%minlength)
      obj%II_planeElementID = obj%FEMdomain%mesh%getElementList(zmin=obj%minlength)

      buf = obj%FEMDomain%mesh%getNodeList( &
            xmin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%xnum)/2.0d0, &
            xmax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%xnum)/2.0d0, &
            ymin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%ynum)/2.0d0, &
            ymax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%ynum)/2.0d0, &
            zmax=0.0d0)
      !obj%A_PointNodeID = buf(1)
      obj%A_PointNodeID = median(buf)

      buf = obj%FEMDomain%mesh%getNodeList( &
            xmin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%xnum)/2.0d0, &
            xmax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%xnum)/2.0d0, &
            ymin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%ynum)/2.0d0, &
            ymax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%ynum)/2.0d0, &
            zmin=obj%minlength)
      !obj%B_PointNodeID = buf(1)
      obj%B_PointNodeID = median(buf)

      buf = obj%FEMDomain%mesh%getElementList( &
            xmin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%xnum)/2.0d0, &
            xmax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%xnum)/2.0d0, &
            ymin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%ynum)/2.0d0, &
            ymax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%ynum)/2.0d0, &
            zmax=0.0d0)

      !obj%A_PointElementID = buf(1)
      obj%A_PointElementID = median(buf)

      buf = obj%FEMDomain%mesh%getElementList( &
            xmin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%xnum)/2.0d0, &
            xmax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%xnum)/2.0d0, &
            ymin=obj%minwidth/2.0d0 - obj%minwidth/dble(obj%ynum)/2.0d0, &
            ymax=obj%minwidth/2.0d0 + obj%minwidth/dble(obj%ynum)/2.0d0, &
            zmin=obj%minlength)
      !obj%B_PointElementID = buf(1)
      obj%B_PointElementID = median(buf)

      !print *, obj%C_PointNodeID
      !print *, obj%D_PointNodeID
      !print *, obj%C_PointElementID
      !print *, obj%D_PointElementID
!
      call obj%FEMdomain%remove()
      
      if (present(species)) then
         if (species == PF_SOYBEAN_CV) then
            call obj%FEMdomain%create(meshtype="Sphere3D", x_num=obj%xnum, y_num=obj%ynum, z_num=obj%xnum, &
                                      x_len=input(default=0.050d0,option=length), &
                                      y_len=input(default=0.050d0,option=thickness), &
                                      z_len=input(default=0.050d0,option=width))
            call obj%FEMDomain%move(&
               x=-(obj%FEMDomain%xmax() + obj%FEMDomain%xmin())/2.0d0,&
               y=-(obj%FEMDomain%ymax() + obj%FEMDomain%ymin())/2.0d0,&
               z=-(obj%FEMDomain%zmax() + obj%FEMDomain%zmin())/2.0d0 )
            
            do i=1,obj%FEMDomain%nn()
               if(obj%FEMDomain%mesh%nodcoord(i,2) >=0.0d0)then
                  obj%FEMDomain%mesh%nodcoord(i,2) = obj%FEMDomain%mesh%nodcoord(i,2)*0.20d0
               endif
            enddo
            
            buf = obj%FEMDomain%mesh%getNodeList( &
               xmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )
            obj%A_PointNodeID = median(buf)
            obj%I_planeNodeID = buf
            
            buf = obj%FEMDomain%mesh%getNodeList( &
               xmax=minval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )
            obj%B_PointNodeID = median(buf)
            obj%II_planeNodeID = buf
            
            buf = obj%FEMDomain%mesh%getNodeList( &
               zmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 3)) &
            )
            obj%C_PointNodeID = median(buf)

            buf = obj%FEMDomain%mesh%getNodeList( &
               zmax=minval(obj%FEMdomain%mesh%nodcoord(:, 3)) &
            )
            obj%D_PointNodeID = median(buf)




            buf = obj%FEMDomain%mesh%getElementList( &
               xmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )
            obj%A_PointElementID = median(buf)
            
            buf = obj%FEMDomain%mesh%getElementList( &
               xmax=maxval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )
            obj%B_PointElementID = median(buf)
            
            buf = obj%FEMDomain%mesh%getElementList( &
               zmax=maxval(obj%FEMdomain%mesh%nodcoord(:, 3)) &
            )
            obj%C_PointElementID = median(buf)

            buf = obj%FEMDomain%mesh%getElementList( &
               zmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 3)) &
            )
            obj%D_PointElementID = median(buf)
            
            return
            !
            !call obj%resize(x=input(default=(obj%minwidth/2.0d0),option=width),&
            !                y=input(default=obj%minthickness/2.0d0,option=thickness), &
            !                z=input(default=obj%minlength,option=thickness))
            !call obj%FEMdomain%move(z=-obj%FEMdomain%zmin())
            !call obj%FEMdomain%move(x=-(obj%FEMdomain%xmax() + obj%FEMdomain%xmin()))
            !call obj%FEMdomain%move(y=-(obj%FEMdomain%ymax() + obj%FEMdomain%ymin()))
         else
            call obj%FEMdomain%create(meshtype="Leaf3D", x_num=obj%xnum, y_num=obj%ynum, z_num=obj%znum, &
                                      x_len=obj%minwidth/2.0d0, y_len=obj%minthickness/2.0d0, &
                                      z_len=obj%minlength, species=species, SoyWidthRatio=SoyWidthRatio)
         end if
      else
         call obj%FEMdomain%create(meshtype="Leaf3D", x_num=obj%xnum, y_num=obj%ynum, z_num=obj%znum, &
                             x_len=obj%minwidth/2.0d0, y_len=obj%minthickness/2.0d0, z_len=obj%minlength, shaperatio=obj%shaperatio)
      end if

      obj%thickness = maxval(obj%femdomain%mesh%nodcoord(:, 2)) &
                      - minval(obj%femdomain%mesh%nodcoord(:, 2))

      buf = obj%FEMDomain%mesh%getNodeList( &
            xmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )
      !obj%A_PointNodeID = buf(1)
      obj%C_PointNodeID = median(buf)

      buf = obj%FEMDomain%mesh%getNodeList( &
            xmax=minval(obj%FEMdomain%mesh%nodcoord(:, 1)) &
            )    !obj%B_PointNodeID = median(buf)
      obj%D_PointNodeID = median(buf)

      buf = obj%FEMDomain%mesh%getElementList( &
            xmin=maxval(obj%FEMdomain%mesh%nodcoord(:, 1)) - obj%minwidth/2.0d0/2.0d0/dble(obj%xnum) &
            )

      !obj%A_PointElementID = median(buf)
      obj%C_PointElementID = median(buf)

      buf = obj%FEMDomain%mesh%getElementList( &
            xmax=minval(obj%FEMdomain%mesh%nodcoord(:, 1)) + obj%minwidth/2.0d0/2.0d0/dble(obj%xnum) &
            )
      !obj%B_PointElementID = median(buf)
      obj%D_PointElementID = median(buf)

      ! Aについて、要素番号、節点番号、要素座標、節点座標のリストを生成

      if (present(regacy)) then
         if (regacy .eqv. .true.) then
            loc(:) = 0.0d0
            if (present(location)) then
               loc(:) = location(:)
            end if
            obj%ShapeFactor = input(default=0.30d0, option=ShapeFactor)
            obj%Thickness = input(default=0.10d0, option=Thickness)
            obj%length = input(default=0.10d0, option=length)
            obj%width = input(default=0.10d0, option=width)

            obj%MaxThickness = input(default=0.10d0, option=MaxThickness)
            obj%Maxlength = input(default=10.0d0, option=Maxlength)
            obj%Maxwidth = input(default=2.0d0, option=Maxwidth)

            obj%outer_normal_bottom(:) = 0.0d0
            obj%outer_normal_bottom(1) = 1.0d0

            obj%outer_normal_top(:) = 0.0d0
            obj%outer_normal_top(1) = 1.0d0

            ! rotate
            obj%outer_normal_Bottom(:) = Rotation3D(vector=obj%outer_normal_bottom, rotx=rotx, roty=roty, rotz=rotz)
            obj%outer_normal_top(:) = Rotation3D(vector=obj%outer_normal_top, rotx=rotx, roty=roty, rotz=rotz)

            obj%center_bottom(:) = loc(:)
            obj%center_top(:) = obj%center_bottom(:) + obj%length*obj%outer_normal_bottom(:)
         end if
      end if
      obj%CarbonDiffusionCoefficient = obj%default_CarbonDiffusionCoefficient*ones(obj%femdomain%ne())

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

   subroutine curveleaf(obj, curvature)
      ! deform by curvature
      class(leaf_), intent(inout) :: obj
      real(real64), intent(in) :: curvature
      real(real64) :: leaf_L, radius, z
      integer(int32) :: i

      if (curvature < dble(1.0e-5)) then
         print *, "Caution >> initLeaf >> curvature is too small < 1.0e-5"
         print *, "Then, ignored."
         return
      end if
      radius = 1.0d0/curvature
      leaf_L = maxval(obj%femdomain%mesh%nodcoord(:, 3)) - minval(obj%femdomain%mesh%nodcoord(:, 3))
      leaf_L = 0.50d0*leaf_L
      do i = 1, obj%femdomain%nn()
         z = obj%femdomain%mesh%nodcoord(i, 3)
         obj%femdomain%mesh%nodcoord(i, 2) = &
            obj%femdomain%mesh%nodcoord(i, 2) &
            - sqrt(radius*radius - leaf_L*leaf_L) &
            + sqrt(radius*radius - (z - leaf_L)*(z - leaf_L))
      end do

   end subroutine

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

      if (present(reset)) then
         if (reset .eqv. .true.) then
            call obj%femdomain%rotate(-obj%rot_x, -obj%rot_y, -obj%rot_z)
            obj%rot_x = 0.0d0
            obj%rot_y = 0.0d0
            obj%rot_z = 0.0d0
         end if
      end if

      origin1 = obj%getCoordinate("A")
      call obj%femdomain%rotate(x, y, z)
      obj%rot_x = obj%rot_x + input(default=0.0d0, option=x)
      obj%rot_y = obj%rot_y + input(default=0.0d0, option=y)
      obj%rot_z = obj%rot_z + input(default=0.0d0, option=z)
      origin2 = obj%getCoordinate("A")
      disp = origin1
      disp(:) = origin1(:) - origin2(:)
      call obj%femdomain%move(x=disp(1), y=disp(2), z=disp(3))

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

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

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

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

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

! ########################################
   subroutine connectleafleaf(obj, direct, leaf)
      class(leaf_), intent(inout) :: obj
      class(leaf_), intent(inout) :: leaf
      character(2), intent(in) :: direct
      real(real64), allocatable :: x1(:), x2(:), disp(:)

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

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

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

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

! ########################################
   subroutine connectLeafStem(obj, direct, Stem)
      class(leaf_), intent(inout) :: obj
      class(Stem_), intent(inout) :: stem
      character(2), intent(in) :: direct
      real(real64), allocatable :: x1(:), x2(:), disp(:)

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

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

! ########################################
   function getCoordinateleaf(obj, nodetype) result(ret)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) :: nodetype
      real(real64), allocatable :: ret(:)
      integer(int32) :: dimnum, n, i

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

      allocate (ret(dimnum))
      ret(:) = 0.0d0

      if (nodetype == "A" .or. nodetype == "a") then
         n = size(obj%I_planeNodeID)
         do i = 1, n
            ret(:) = ret(:) + obj%femdomain%mesh%nodcoord(obj%I_planeNodeID(i), :)
         end do
         ret(:) = 1.0d0/dble(n)*ret(:)

         !ret = obj%femdomain%mesh%nodcoord(obj%A_PointNodeID,:)
      end if
      if (nodetype == "B" .or. nodetype == "b") then
         !ret = obj%femdomain%mesh%nodcoord(obj%B_PointNodeID,:)
         n = size(obj%II_planeNodeID)
         do i = 1, n
            ret(:) = ret(:) + obj%femdomain%mesh%nodcoord(obj%II_planeNodeID(i), :)
         end do
         ret(:) = 1.0d0/dble(n)*ret(:)
      end if

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

! ########################################
   subroutine gmshleaf(obj, name)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%gmsh(Name=name)
      ! PPFD を出力
      call obj%femdomain%gmsh(Name=name//"_PPFD_", field=obj%PPFD)
      ! ソース量 を出力
      call obj%femdomain%gmsh(Name=name//"_SOURCE_", field=obj%source)
      ! 光合成速度 を出力
      call obj%femdomain%gmsh(Name=name//"_A_", field=obj%A)

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

! ########################################
   subroutine mshleaf(obj, name)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%msh(Name=name)
      ! PPFD を出力
      !call obj%femdomain%msh(Name=name//"_PPFD_",field=obj%PPFD)
      ! ソース量 を出力
      !call obj%femdomain%msh(Name=name//"_SOURCE_",field=obj%source)
      ! 光合成速度 を出力
      !call obj%femdomain%msh(Name=name//"_A_",field=obj%A)

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

! ########################################
   subroutine vtkleaf(obj, name, field_name)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) ::name
      character(*), optional, intent(in) ::field_name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%vtk(Name=name, field=field_name)
      ! PPFD を出力
      !call obj%femdomain%msh(Name=name//"_PPFD_",field=obj%PPFD)
      ! ソース量 を出力
      !call obj%femdomain%msh(Name=name//"_SOURCE_",field=obj%source)
      ! 光合成速度 を出力
      !call obj%femdomain%msh(Name=name//"_A_",field=obj%A)

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

! ########################################
   subroutine stlleaf(obj, name)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%stl(Name=name)
      ! PPFD を出力
      !call obj%femdomain%msh(Name=name//"_PPFD_",field=obj%PPFD)
      ! ソース量 を出力
      !call obj%femdomain%msh(Name=name//"_SOURCE_",field=obj%source)
      ! 光合成速度 を出力
      !call obj%femdomain%msh(Name=name//"_A_",field=obj%A)

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

! ########################################
   subroutine plyleaf(obj, name)
      class(leaf_), intent(inout) :: obj
      character(*), intent(in) ::name
      if (obj%femdomain%mesh%empty()) then
         return
      end if

      call obj%femdomain%ply(Name=name)
      ! PPFD を出力
      !call obj%femdomain%msh(Name=name//"_PPFD_",field=obj%PPFD)
      ! ソース量 を出力
      !call obj%femdomain%msh(Name=name//"_SOURCE_",field=obj%source)
      ! 光合成速度 を出力
      !call obj%femdomain%msh(Name=name//"_A_",field=obj%A)

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

! ########################################
   subroutine resizeleaf(obj, x, y, z, Length, Width)
      class(Leaf_), optional, intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z, Length, Width
      real(real64), allocatable :: origin1(:), origin2(:), disp(:), total_rot(:)

      if (present(Length) .or. present(Width)) then

         total_rot = zeros(3)
         total_rot(1) = obj%rot_x
         total_rot(2) = obj%rot_y
         total_rot(3) = obj%rot_z
         call obj%rotate(z=-total_rot(3))
         call obj%rotate(y=-total_rot(2))
         call obj%rotate(x=-total_rot(1))

         if (present(Length)) then
            call obj%femdomain%resize(z=Length)
         end if

         if (present(Width)) then
            call obj%femdomain%resize(x=Width)
         end if

         call obj%rotate(x=total_rot(1))
         call obj%rotate(y=total_rot(2))
         call obj%rotate(z=total_rot(3))

      else
         origin1 = obj%getCoordinate("A")
         call obj%femdomain%resize(x_len=x, y_len=y, z_len=z)
         origin2 = obj%getCoordinate("A")
         disp = origin1 - origin2
         call obj%move(x=disp(1), y=disp(2), z=disp(3))
      end if
   end subroutine
! ########################################

! ########################################
   function getLengthLeaf(obj) result(Length)
      class(Leaf_), optional, intent(inout) :: obj
      real(real64) :: Length
      real(real64), allocatable :: origin1(:), origin2(:), disp(:), total_rot(:)

      total_rot = zeros(3)
      total_rot(1) = obj%rot_x
      total_rot(2) = obj%rot_y
      total_rot(3) = obj%rot_z
      call obj%rotate(z=-total_rot(3))
      call obj%rotate(y=-total_rot(2))
      call obj%rotate(x=-total_rot(1))

      Length = obj%femdomain%zmax() - obj%femdomain%zmin()

      call obj%rotate(x=total_rot(1))
      call obj%rotate(y=total_rot(2))
      call obj%rotate(z=total_rot(3))

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

! ########################################
   function FullyExpandedLeaf(obj, threshold) result(ret_expanded)
      class(Leaf_), optional, intent(inout) :: obj
      real(real64), intent(in) :: threshold
      logical :: ret_expanded
      real(real64) :: length, full_length

      if (obj%getLength()/obj%final_length > threshold) then
         ret_expanded = .true.
      else
         ret_expanded = .false.
      end if

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

! ########################################
   subroutine rescaleleaf(obj, x, y, z)
      class(Leaf_), optional, intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z
      real(real64), allocatable :: origin1(:), origin2(:), disp(:)

      origin1 = obj%getCoordinate("A")
      call obj%femdomain%resize(x_rate=x, y_rate=y, z_rate=z)
      origin2 = obj%getCoordinate("A")
      disp = origin1 - origin2
      call obj%move(x=disp(1), y=disp(2), z=disp(3))
   end subroutine
! ########################################

! ########################################
!subroutine LayTracingLeaf(obj,maxPPFD,light,)
!    class(Leaf_),intent(inout) :: obj
!    class(Light_),intent(in) :: light
!    real(real64),intent(in) :: maxPPFD
!    integer(int32) :: i,j,n,m,node_id
!    real(real64) :: lx(3)
!    real(real64),allocatable :: Elem_x(:,:)
!    ! PPFDを計算する。
!    ! Photosynthetic photon flux density (PPFD)
!    ! micro-mol/m^2/s
!
!    ! 反射、屈折は無視、直線のみ
!
!    n=size(obj%FEMDomain%Mesh%ElemNod,2)
!    m=size(obj%FEMDomain%Mesh%NodCoord,2)
!
!    allocate(Elem_x(n,m) )
!    ! 要素ごと
!    do i=1, size(obj%FEMDomain%Mesh%ElemNod,1)
!        do j=1,size(obj%FEMDomain%Mesh%ElemNod,2)
!            node_id = obj%FEMDomain%Mesh%ElemNod(i,j)
!            Elem_x(j,:) = obj%FEMDomain%Mesh%NodCoord(node_id,:)
!        enddo
!        ! 要素座標 >> Elem_x(:,:)
!        ! 光源座標 >> lx(:)
!
!    enddo
!
!
!
!end subroutine
! ########################################

! ########################################
   subroutine photosynthesisLeaf(obj, dt, air)

      ! https://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/39102/1/67-013.pdf

      class(Leaf_), intent(inout) :: obj
      type(Air_), intent(in) :: air
      type(IO_) :: f
      real(real64), intent(in) :: dt
      ! Farquhar modelのパラメータ
      real(real64) :: A   ! CO2吸収速度
      real(real64) :: V_c ! カルボキシル化反応速度
      real(real64) :: V_o ! 酸素化反応速度

      real(real64) :: W_c! RuBPが飽和している場合のCO2吸収速度
      real(real64) :: W_j! RuBP供給が律速している場合のCO2吸収速度

      real(real64) :: V_cmax ! 最大カルボキシル化反応速度
      real(real64) :: V_omax ! 最大酸素化反応速度
      real(real64) :: O2 ! 酸素濃度
      real(real64) :: CO2 ! 二酸化炭素濃度
      real(real64) :: R_d ! なんだっけ

      real(real64) :: K_c ! CO2に対するミカエリス定数
      real(real64) :: K_o ! O2に対するミカエリス定数

      real(real64) :: J_ ! 電子伝達速度
      real(real64) :: I_ ! 光強度
      real(real64) :: phi ! I-J曲線の初期勾配
      real(real64) :: J_max !最大電子伝達速度
      real(real64) :: theta_r ! 曲線の凸度

      real(real64) :: pfd

      real(real64) :: Lambda, volume, gs, Ci, Ca

      integer(int32) :: i, element_id

      obj%temp = air%temp
      obj%CO2 = air%CO2
      obj%O2 = air%O2

      ! TT-model
      do i = 1, size(obj%source)
         ! 要素ごとに電子伝達速度を求める
         element_id = i
         pfd = obj%ppfd(element_id)

         ! I-J関係
         obj%J_ = 0.240d0*pfd/(sqrt(1.0d0 + (0.240d0*0.240d0)*pfd*pfd)/obj%J_max/obj%J_max)

         ! lambdaからV_omaxを推定
         obj%V_omax = obj%Lambda*(2.0d0*obj%V_cmax*obj%K_o)/(obj%K_c*O2)

         ! CO2固定速度の計算
         V_c = (obj%V_cmax*obj%CO2)/(obj%CO2 + obj%K_o*(1.0d0 + obj%O2/obj%K_o))
         V_o = (obj%V_omax*obj%O2)/(obj%O2 + obj%K_o*(1.0d0 + obj%CO2/obj%K_c))

         ! RuBPが飽和している場合のCO2吸収速度
         Ca = obj%CO2
         W_c = (obj%V_cmax*(Ca - obj%Lambda))/(Ca + obj%K_c*(1.0d0 + obj%O2/obj%K_o))

         ! RuBP供給が律速している場合のCO2吸収速度
         W_j = obj%J_*(Ca - obj%Lambda)/(4.0d0*Ca + 8.0d0*obj%Lambda) - obj%R_d

         if (W_j >= W_c) then
            A = W_c
         else
            A = W_j
         end if

         !A = gs(Ca - Ci)
         !A = gs*Ca - gs*Ci
         !A = gs*Ca - gs*Ci
         !gs_max = 0.6

         ! 要素体積を求める, m^3
         obj%A(element_id) = A
         volume = obj%femdomain%getVolume(elem=element_id)

         ! CO2拡散が考慮されていない，表皮がない．コンダクタンス∞
         ! コンダクタンスを光で．
         ! 水分の関数をあとでいれる

         !CO2固定量　mincro-mol/m-2/s
         ! ここ、体積あたりにする必要がある
         ! 一応、通常の葉の厚さを2mmとして、
         ! 1 micro-mol/m^2/sを、 1 micro-mol/ 0.0002m^3/s= 5000micro-mol/m^3/sとして計算
         ! また、ソース量はC6H12O6の質量gramとして換算する。
         ! CO2の分子量44.01g/mol
         ! C6H12O6の分子量180.16g/mol
         ! 6CO2 + 12H2O => C6H12O6 + 6H2O + 6O2
         ! よって、生成されるソース量は
         !               {CO2固定量,mol     }× {1/6 してグルコースmol}×グルコース分子量
         obj%source(i) = obj%source(i) + A*dt*5000.0d0*volume*1.0d0/6.0d0*180.160d0

      end do
!    ! For each elements, estimate photosynthesis by Farquhar model
!    do i=1,size(obj%source)
!
!        ! 光合成量の計算
!        ! Farquhar model
!        V_c = (V_cmax*CO2)/(CO2 + K_o * (1.0d0 + O2/K_o) )
!        V_o = (V_omax*O2 )/(O2 + K_o * (1.0d0 + CO2/K_c) )
!
!        Lambda = (V_omax*K_c*O2)/( 2.0d0 * V_cmax*K_o )
!
!        W_c = (V_cmax*(CO2 - Lambda))/(CO2 + K_c*(1.0d0 + O2/K_o)  )
!
!        J_ = (phi*I_ + J_max - &
!        sqrt( (phi*I_ + J_max)**(2.0d0) - 4.0d0*phi*I_*theta_r*J_max)&
!        /(2.0d0 * theta_r) )
!        W_j = J_ * (CO2 - Lambda)/(4.0d0 * CO2 + 8.0d0 * Lambda ) - R_d
!        ! CO2吸収速度
!        A = V_c + 0.50d0*V_o - R_d
!
!        if(W_j >= W_c )then
!            A = W_c
!        else
!            A = W_j
!        endif
!
!
!    enddo
!
!
   end subroutine

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

! ########################################
   function getPhotosynthesisSpeedPerVolumeLeaf(obj, dt, air) result(Speed_PV)

      ! https://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/39102/1/67-013.pdf

      class(Leaf_), intent(inout) :: obj
      type(Air_), intent(in) :: air
      type(IO_) :: f
      real(real64), intent(in) :: dt
      ! Farquhar modelのパラメータ
      real(real64) :: A   ! CO2吸収速度
      real(real64) :: V_c ! カルボキシル化反応速度
      real(real64) :: V_o ! 酸素化反応速度

      real(real64) :: W_c! RuBPが飽和している場合のCO2吸収速度
      real(real64) :: W_j! RuBP供給が律速している場合のCO2吸収速度

      real(real64) :: V_cmax ! 最大カルボキシル化反応速度
      real(real64) :: V_omax ! 最大酸素化反応速度
      real(real64) :: O2 ! 酸素濃度
      real(real64) :: CO2 ! 二酸化炭素濃度
      real(real64) :: R_d ! 呼吸速度

      real(real64) :: K_c ! CO2に対するミカエリス定数
      real(real64) :: K_o ! O2に対するミカエリス定数

      real(real64) :: J_ ! 電子伝達速度
      real(real64) :: I_ ! 光強度
      real(real64) :: phi ! I-J曲線の初期勾配
      real(real64) :: J_max !最大電子伝達速度
      real(real64) :: theta_r ! 曲線の凸度

      real(real64) :: pfd

      real(real64) :: Lambda, volume
      real(real128) :: buf
      real(real64), allocatable :: Speed_PV(:)

      integer(int32) :: i, element_id

      obj%temp = air%temp
      obj%CO2 = air%CO2
      obj%O2 = air%O2

      Speed_PV = zeros(obj%femdomain%ne())
      ! TT-model
      do i = 1, size(Speed_PV)
         ! 要素ごとに電子伝達速度を求める
         element_id = i
         pfd = obj%ppfd(element_id)
         obj%J_ = 0.240d0*pfd/(sqrt(1.0d0 + (0.240d0*0.240d0)*pfd*pfd)/obj%J_max/obj%J_max)

         ! lambdaからV_omaxを推定
         !obj%V_omax = obj%Lambda* 2.0d0 * obj%V_cmax*obj%K_o/(obj%K_c*O2)
         buf = obj%Lambda*2.0d0*obj%V_cmax*obj%K_o
         obj%V_omax = dble(buf/obj%K_c/obj%O2)

         ! CO2固定速度の計算
         V_c = (obj%V_cmax*obj%CO2)/(obj%CO2 + obj%K_o*(1.0d0 + obj%O2/obj%K_o))
         V_o = (obj%V_omax*obj%O2)/(obj%O2 + obj%K_o*(1.0d0 + obj%CO2/obj%K_c))

         ! RuBPが飽和している場合のCO2吸収速度
         W_c = (obj%V_cmax*(obj%CO2 - obj%Lambda))/(obj%CO2 + obj%K_c*(1.0d0 + obj%O2/obj%K_o))

         ! RuBP供給が律速している場合のCO2吸収速度
         W_j = obj%J_*(obj%CO2 - obj%Lambda)/(4.0d0*obj%CO2 + 8.0d0*obj%Lambda) - obj%R_d

         if (W_j >= W_c) then
            A = W_c
         else
            A = W_j
         end if
         ! 要素体積を求める, m^3
         obj%A(element_id) = A
         volume = obj%femdomain%getVolume(elem=element_id)

         !CO2固定量　mincro-mol/m^2/s
         ! ここ、体積あたりにする必要がある
         ! 一応、通常の葉の厚さを0.2mmとして、
         ! 1 micro-mol/m^2/sを、 1 micro-mol/ 0.002m^3/s= 500micro-mol/m^3/sとして計算
         ! また、ソース量はC6H12O6の質量gramとして換算する。
         ! CO2の分子量44.01g/mol
         ! C6H12O6の分子量180.16g/mol
         ! 6CO2 + 12H2O => C6H12O6 + 6H2O + 6O2
         ! よって、生成されるソース量(micro-gram)は
         !               {CO2固定量,mol     }× {1/6 してグルコースmol}×グルコース分子量
         Speed_PV(i) = A*dt*5000.0d0*1.0d0/6.0d0*180.160d0

      end do
   end function

! ####################################################################
   subroutine adjustLeaf(obj, width)
      class(Leaf_), intent(inout) :: obj
      real(real64), intent(in) :: width(:, :)

   end subroutine

   function getVolumeLeaf(obj) result(ret)
      class(Leaf_), intent(in) :: obj
      real(real64) :: ret
      integer(int32) :: i, j

      ret = 0.0d0

      if (obj%femdomain%mesh%empty()) then
         return
      end if

      do i = 1, obj%femdomain%ne()
         ret = ret + obj%femdomain%getVolume(elem=i)
      end do

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

! ########################################
   function getBiomassLeaf(obj) result(ret)
      class(Leaf_), intent(in) :: obj
      real(real64) :: ret
      integer(int32) :: i, j

      ret = 0.0d0

      if (obj%femdomain%mesh%empty()) then
         return
      end if

      do i = 1, obj%femdomain%ne()
         ret = ret + obj%femdomain%getVolume(elem=i)*obj%drydensity(i)
      end do

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

! ########################################
   function emptyLeaf(obj) result(leaf_is_empty)
      class(Leaf_), intent(in) :: obj
      logical :: leaf_is_empty

      leaf_is_empty = obj%femdomain%empty()

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

   function getLeafAreaLeaf(obj) result(LeafArea)
      class(Leaf_), intent(in) :: obj
      real(real64) :: LeafArea

      LeafArea = obj%getVolume()/obj%thickness

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

! ################################################################
   pure function getRadiusLeaf(obj) result(radius)
      class(Leaf_), intent(in) :: obj
      real(real64), allocatable :: Points(:, :)
      real(real64) :: radius, center(3)

      Points = obj%femdomain%mesh%nodcoord

      ! search Intersect leaf
      center(1) = sum(Points(:, 1))/dble(size(Points, 1))
      center(2) = sum(Points(:, 2))/dble(size(Points, 1))
      center(3) = sum(Points(:, 3))/dble(size(Points, 1))

      Points(:, 1) = Points(:, 1) - center(1)
      Points(:, 2) = Points(:, 2) - center(2)

      radius = maxval(Points(:, 1)*Points(:, 1) + Points(:, 2)*Points(:, 2))
      radius = sqrt(radius)

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

! ################################################################
   pure function getCenterLeaf(obj) result(Center)
      class(Leaf_), intent(in) :: obj
      real(real64), allocatable :: Points(:, :)
      real(real64) :: center(3)

      Points = obj%femdomain%mesh%nodcoord

      ! search Intersect leaf
      center(1) = sum(Points(:, 1))/dble(size(Points, 1))
      center(2) = sum(Points(:, 2))/dble(size(Points, 1))
      center(3) = sum(Points(:, 3))/dble(size(Points, 1))

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

   subroutine syncleaf(obj, from, mpid)
      class(Leaf_), intent(inout) :: obj
      integer(int32), intent(in) :: from
      type(MPI_), intent(inout) :: mpid

      call obj%FEMDomain%sync(from=from, mpid=mpid)
      call mpid%bcast(from=from, val=obj%LeafSurfaceNode2D)
      call mpid%bcast(from=from, val=obj%ShapeFactor) !
      call mpid%bcast(from=from, val=obj%Thickness) !
      call mpid%bcast(from=from, val=obj%length) !
      call mpid%bcast(from=from, val=obj%width) !
      call mpid%BcastMPIRealVecFixedSize(from=from, val=obj%center) !(3)
      call mpid%bcast(from=from, val=obj%MaxThickness) !
      call mpid%bcast(from=from, val=obj%Maxlength) !
      call mpid%bcast(from=from, val=obj%Maxwidth) !
      call mpid%BcastMPIRealVecFixedSize(from=from, val=obj%center_bottom) !(3)
      call mpid%BcastMPIRealVecFixedSize(from=from, val=obj%center_top) !(3)
      call mpid%BcastMPIRealVecFixedSize(from=from, val=obj%outer_normal_bottom) !(3)
      call mpid%BcastMPIRealVecFixedSize(from=from, val=obj%outer_normal_top) !(3)
      call mpid%bcast(from=from, val=obj%source)
      call mpid%bcast(from=from, val=obj%ppfd)
      call mpid%bcast(from=from, val=obj%A)
      call mpid%bcast(from=from, val=obj%Division)
      !type(leaf_),pointer ::  pleaf
      !type(Peti_),pointer ::  pPeti
      call mpid%bcast(from=from, val=obj%rot_x) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%rot_y) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%rot_z) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%disp_x) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%disp_y) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%disp_z) ! = 0.0d0
      call mpid%bcast(from=from, val=obj%shaperatio) ! = 0.30d0
      call mpid%bcast(from=from, val=obj%minwidth) !
      call mpid%bcast(from=from, val=obj%minlength) !
      call mpid%bcast(from=from, val=obj%MinThickness) !

      call mpid%bcast(from=from, val=obj%I_planeNodeID) !(:)
      call mpid%bcast(from=from, val=obj%I_planeElementID) !(:)
      call mpid%bcast(from=from, val=obj%II_planeNodeID) !(:)
      call mpid%bcast(from=from, val=obj%II_planeElementID) !(:)
      call mpid%bcast(from=from, val=obj%A_PointNodeID)!
      call mpid%bcast(from=from, val=obj%B_PointNodeID)!
      call mpid%bcast(from=from, val=obj%A_PointElementID)!
      call mpid%bcast(from=from, val=obj%B_PointElementID)!
      call mpid%bcast(from=from, val=obj%xnum)! = 10
      call mpid%bcast(from=from, val=obj%ynum)! = 10
      call mpid%bcast(from=from, val=obj%znum)! = 10

      ! phisiological parameters

      call mpid%bcast(from=from, val=obj%V_cmax) ! = 100.0d0 ! 最大カルボキシル化反応速度, mincro-mol/m-2/s
      call mpid%bcast(from=from, val=obj%V_omax) ! = 100.0d0 ! 最大酸素化反応速度, mincro-mol/m-2/s, lambdaから推定
      call mpid%bcast(from=from, val=obj%O2) ! = 380.0d0! 酸素濃度, ppm
      call mpid%bcast(from=from, val=obj%CO2) !=202000.0d0! 二酸化炭素濃度, ppm
      call mpid%bcast(from=from, val=obj%R_d) !=1.0d0 ! 暗呼吸速度, mincro-mol/m-2/s

      call mpid%bcast(from=from, val=obj%K_c) !=272.380d0 ! CO2に対するミカエリス定数
      call mpid%bcast(from=from, val=obj%K_o) !=165820.0d0 ! O2に対するミカエリス定数

      call mpid%bcast(from=from, val=obj%J_) !=0.0d0 ! 電子伝達速度
      call mpid%bcast(from=from, val=obj%I_) !=0.0d0 ! 光強度
      call mpid%bcast(from=from, val=obj%phi) !=0.0d0 ! I-J曲線の初期勾配
      call mpid%bcast(from=from, val=obj%J_max) !=180.0d0 !最大電子伝達速度,mincro-mol/m-2/s
      call mpid%bcast(from=from, val=obj%theta_r) !=0.0d0 ! 曲線の凸度

      call mpid%bcast(from=from, val=obj%maxPPFD) !=1.0d0 ! micro-mol/m^2/s

      call mpid%bcast(from=from, val=obj%Lambda) != 37.430d0 ! 暗呼吸速度を無視した時のCO2補償点ppm
      call mpid%bcast(from=from, val=obj%temp) !=303.0d0 ! temp

      ! physical parameter
      call mpid%bcast(from=from, val=obj%DryDensity)
      call mpid%bcast(from=from, val=obj%WaterContent)

      ! For deformation analysis
      call mpid%bcast(from=from, val=obj%YoungModulus)

      call mpid%bcast(from=from, val=obj%PoissonRatio)
      call mpid%bcast(from=from, val=obj%CarbonDiffusionCoefficient)
      call mpid%bcast(from=from, val=obj%Density)
      call mpid%bcast(from=from, val=obj%Stress)
      call mpid%bcast(from=from, val=obj%Displacement)

      call mpid%bcast(from=from, val=obj%BoundaryTractionForce)
      call mpid%bcast(from=from, val=obj%BoundaryDisplacement)

   end subroutine

   subroutine syncLeafVector(obj, from, mpid)
      type(Leaf_), allocatable, intent(inout) :: obj(:)
      integer(int32), intent(in) :: from
      type(MPI_), intent(inout) :: mpid
      integer(int32) :: vec_size, i

      vec_size = 0
      if (mpid%myrank == from) then
         if (.not. allocated(obj)) then
            vec_size = -1
         else
            vec_size = size(obj)
         end if
      end if
      call mpid%bcast(from=from, val=vec_size)
      if (vec_size < 1) then
         return
      end if

      if (from /= mpid%myrank) then
         if (allocated(obj)) then
            deallocate (obj)
         end if
         allocate (obj(vec_size))
      end if

      do i = 1, vec_size
         call obj(i)%sync(from=from, mpid=mpid)
      end do

   end subroutine
! ##################################################################
   function getNormalVectorLeaf(obj, ElementID) result(ret)
      class(Leaf_), intent(inout) :: obj
      integer(int32), intent(in) :: ElementID
      real(real64), allocatable :: ret(:), x_A(:), x_B(:), x_C(:), x_D(:), x(:), &
                                   n_AC(:), n_CB(:), n_BD(:), n_DA(:)
      real(real64) :: min_norm
      integer(int32) :: num_effective_norms, nd

      nd = obj%femdomain%nd()
      allocate (ret(nd))
      x_A = obj%femdomain%centerPosition(ElementID=obj%A_PointElementID)
      x_B = obj%femdomain%centerPosition(ElementID=obj%B_PointElementID)
      x_C = obj%femdomain%centerPosition(ElementID=obj%C_PointElementID)
      x_D = obj%femdomain%centerPosition(ElementID=obj%D_PointElementID)
      x = obj%femdomain%centerPosition(ElementID=ElementID)

      n_AC = cross_product(x_C - x_A, x - x_C)
      n_CB = cross_product(x_B - x_C, x - x_B)
      n_BD = cross_product(x_D - x_B, x - x_D)
      n_DA = cross_product(x_A - x_D, x - x_A)

      num_effective_norms = 0
      min_norm = dble(1.0e-6)
      if (norm(n_AC) < min_norm) then
         n_AC = zeros(nd)
      else
         num_effective_norms = num_effective_norms + 1
         n_AC = n_AC(:)/norm(n_AC)
      end if

      if (norm(n_CB) < min_norm) then
         n_CB = zeros(nd)
      else
         num_effective_norms = num_effective_norms + 1
         n_CB = n_CB(:)/norm(n_CB)
      end if

      if (norm(n_BD) < min_norm) then
         n_BD = zeros(nd)
      else
         num_effective_norms = num_effective_norms + 1
         n_BD = n_BD(:)/norm(n_BD)
      end if

      if (norm(n_DA) < dble(1.0e-6)) then
         n_DA = zeros(nd)
      else
         num_effective_norms = num_effective_norms + 1
         n_DA = n_DA(:)/norm(n_DA)
      end if

      if (num_effective_norms == 0) then
         print *, "ERROR :: getNormalVectorLeaf >> no valid normal vector is found"
         stop
      end if
      ret = (n_AC + n_CB + n_BD + n_DA)/dble(num_effective_norms)
      ret = ret/norm(ret)*(-1.0d0)
      ! get outer nomal of the element

      !
      !           D %%%%%%%%%%%%%%%%%%%%%%%%%%%  B
      !         %%                        %   %
      !        %%                    %      %%
      !      %%                 %          %%
      !     %%            %              %%
      !     %%      %                  %%
      !     %%                       %%
      !   A   %%                  %%
      !      <I> %%%%%%%%%%%%%%%% C

   end function
! #################################################################
   subroutine change_length_or_width_Leaf(this, dt)
      class(Leaf_), intent(inout)::this
      real(real64), intent(in) :: dt
      real(real64) :: Length, Width

      if (this%already_grown) then
         ! ignore growth for this
         return
      end if

      if (this%femdomain%empty()) then
         return
      end if

      ! logistic curve
      ! automatic growth
      this%my_time = this%my_time + dt
      ! growth curve: logistic function
      Length = this%final_length &
               /(1.0d0 + &
                 (this%final_length/this%initial_length - 1.0d0) &
                 *exp(-this%length_growth_ratio*this%my_time))

      Width = this%final_Width &
              /(1.0d0 + &
                (this%final_Width/this%initial_Width - 1.0d0) &
                *exp(-this%Width_growth_ratio*this%my_time))

      call this%resize(Length=Length, Width=Width)
      return

   end subroutine

! ###############################################################
   subroutine removeLeaf(this)
      class(Leaf_), intent(inout) :: this

      call this%FEMDomain%remove()
      if (allocated(this%LeafSurfaceNode2D)) deallocate (this%LeafSurfaceNode2D)
      this%ShapeFactor = 0.0d0
      this%Thickness = 0.0d0
      this%length = 0.0d0
      this%width = 0.0d0
      this%center = 0.0d0
      this%MaxThickness = 0.0d0
      this%Maxlength = 0.0d0
      this%Maxwidth = 0.0d0
      this%center_bottom = 0.0d0
      this%center_top = 0.0d0
      this%outer_normal_bottom = 0.0d0
      this%outer_normal_top = 0.0d0
      if (allocated(this%source)) deallocate (this%source)

      this%Division = 0
      if (associated(this%pleaf)) nullify (this%pleaf)
      if (associated(this%pPeti)) nullify (this%pPeti)

      this%rot_x = 0.0d0
      this%rot_y = 0.0d0
      this%rot_z = 0.0d0
      this%disp_x = 0.0d0
      this%disp_y = 0.0d0
      this%disp_z = 0.0d0
      this%shaperatio = 0.30d0
      this%minwidth = 0.0d0
      this%minlength = 0.0d0
      this%MinThickness = 0.0d0

      ! id in multi-leaf
      this%LeafID = -1
      this%already_grown = .false.

      if (allocated(this%I_planeNodeID)) deallocate (this%I_planeNodeID)
      if (allocated(this%I_planeElementID)) deallocate (this%I_planeElementID)
      if (allocated(this%II_planeNodeID)) deallocate (this%II_planeNodeID)
      if (allocated(this%II_planeElementID)) deallocate (this%II_planeElementID)
      this%A_PointNodeID = 0
      this%B_PointNodeID = 0
      this%C_PointNodeID = 0
      this%D_PointNodeID = 0
      this%A_PointElementID = 0
      this%B_PointElementID = 0
      this%C_PointElementID = 0
      this%D_PointElementID = 0
      this%xnum = 10
      this%ynum = 10
      this%znum = 10

      ! phisiological parameters

      this%V_cmax = 100.0d0 ! 最大カルボキシル化反応速度, mincro-mol/m-2/s
      this%V_omax = 100.0d0 ! 最大酸素化反応速度, mincro-mol/m-2/s, lambdaから推定
      this%O2 = 380.0d0! 酸素濃度, ppm
      this%CO2 = 202000.0d0! 二酸化炭素濃度, ppm
      this%R_d = 1.0d0 ! 暗呼吸速度, mincro-mol/m-2/s

      this%K_c = 272.380d0 ! CO2に対するミカエリス定数
      this%K_o = 165820.0d0 ! O2に対するミカエリス定数

      this%J_ = 0.0d0 ! 電子伝達速度
      this%I_ = 0.0d0 ! 光強度
      this%phi = 0.0d0 ! I-J曲線の初期勾配
      this%J_max = 180.0d0 !最大電子伝達速度,mincro-mol/m-2/s
      this%theta_r = 0.0d0 ! 曲線の凸度

      this%maxPPFD = 1000.0d0 ! micro-mol/m^2/s

      this%Lambda = 37.430d0 ! 暗呼吸速度を無視した時のCO2補償点ppm
      this%temp = 303.0d0 ! temp

      ! physical parameter
      if (allocated(this%DryDensity)) deallocate (this%DryDensity)
      if (allocated(this%WaterContent)) deallocate (this%WaterContent)

      ! For deformation analysis
      if (allocated(this%YoungModulus)) deallocate (this%YoungModulus)
      if (allocated(this%PoissonRatio)) deallocate (this%PoissonRatio)
      if (allocated(this%Density)) deallocate (this%Density)
      if (allocated(this%CarbonDiffusionCoefficient)) deallocate (this%CarbonDiffusionCoefficient)
      if (allocated(this%Stress)) deallocate (this%Stress)
      if (allocated(this%Displacement)) deallocate (this%Displacement)

      if (allocated(this%BoundaryTractionForce)) deallocate (this%BoundaryTractionForce)
      if (allocated(this%BoundaryDisplacement)) deallocate (this%BoundaryDisplacement)

      ! growth parameters
      this%my_time = 0.0d0
      this%initial_width = 0.0010d0 ! 1.0 mm
      this%initial_length = 0.0010d0 ! 1.0 mm
      this%final_width = 0.060d0   !  60.0 mm
      this%final_length = 0.100d0   ! 100.0 mm
      this%width_growth_ratio = 1.0d0/4.0d0   !
      this%length_growth_ratio = 1.0d0/4.0d0   !

   end subroutine

   function nnLeaf(this) result(ret)
      class(Leaf_), intent(in) :: this
      integer(int32) :: ret

      ret = this%femdomain%nn()
   end function

   function neLeaf(this) result(ret)
      class(Leaf_), intent(in) :: this
      integer(int32) :: ret

      ret = this%femdomain%ne()
   end function

! ############################################################
   function append_leaf_object_vector(arg1, arg2) result(ret)
      type(Leaf_), allocatable, intent(in) :: arg1(:), arg2(:)
      type(Leaf_), allocatable :: ret(:)

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

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

end module
