module growthmod implicit none contains !=================================== subroutine initialize_gw3(d_connect, d_node, peti_coord, leaf_coord, max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate, peti_length, max_length) integer, allocatable, intent(out)::d_connect(:, :) real(8), allocatable, intent(out)::d_node(:, :), peti_coord(:, :, :), leaf_coord(:, :, :, :) real(8), intent(in)::max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate, peti_length, max_length real(8) r1, r2, x_coord, & y_coord, z_coord, direction, pi, angle, x_coord_1, y_coord_1, z_coord_1, & peti_RUE, mid_length, width, length real(8), allocatable::a(:), b(:), z(:), nvec(:), a_unit(:), b_unit(:) integer int_num, real_value integer i, j, n, m, o, parent, itr integer :: seedsize integer, allocatable :: seed(:) real :: rnd integer :: c !時間を入れる pi = 3.14159265350d0 int_num = 3 real_value = 5 allocate (d_connect(2, int_num), d_node(2, real_value)) allocate (peti_coord(2, 4, 3), leaf_coord(2, 3, 4, 3)) allocate (a(3), b(3), z(3), nvec(3), a_unit(3), b_unit(3)) leaf_coord(:, :, :, :) = 0.0d0 peti_coord(:, :, :) = 0.0d0 !initisalize d_connect d_connect(1, 1) = 1 d_connect(1, 2) = -1 !Apical No. d_connect(1, 3) = 0 d_connect(2, 1) = 1 d_connect(2, 2) = 0 !Apical No. d_connect(2, 3) = 0 !initialize d_node d_node(1, 1) = 0.0d0 d_node(1, 2) = 0.0d0 d_node(1, 3) = 0.0d0 d_node(1, 4) = 1.0d0 d_node(1, 5) = 0.0d0 d_node(2, 1) = 0.0d0 d_node(2, 2) = 0.0d0 d_node(2, 3) = max_length d_node(2, 4) = 1.0d0 d_node(2, 5) = 0.0d0 !initialize peti_coord(:,:,:)============ !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) length = ((rnd - 0.50d0) + 1.0d0)*peti_length call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) direction = rnd*360.0d0 r1 = cos(pi*(angle/180.0d0))*length r2 = sin(pi*(angle/180.0d0))*length x_coord = d_node(2, 1) y_coord = d_node(2, 2) z_coord = d_node(2, 3) x_coord = x_coord + r2*cos(pi*(direction/180.0d0)) y_coord = y_coord + r2*sin(pi*(direction/180.0d0)) z_coord = z_coord + r1 peti_coord(1, :, :) = 0.0d0 peti_coord(2, 1, 1) = x_coord peti_coord(2, 1, 2) = y_coord peti_coord(2, 1, 3) = z_coord !second-fourth do i = 2, 4 !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 2) then peti_RUE = 0.1 else peti_RUE = 0.05 end if length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length*peti_RUE call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 2) then direction = rnd*360.0d0 elseif (i == 3) then direction = direction + 120 elseif (i == 4) then direction = direction - 240 end if r1 = cos(pi*(angle/180.0d0))*length r2 = sin(pi*(angle/180.0d0))*length x_coord_1 = x_coord y_coord_1 = y_coord z_coord_1 = z_coord x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0)) y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0)) z_coord_1 = z_coord_1 + r1 peti_coord(2, i, 1) = x_coord_1 peti_coord(2, i, 2) = y_coord_1 peti_coord(2, i, 3) = z_coord_1 end do !leaf generate do i = 1, 3 !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length width = ((rnd - 0.50d0) + 1.0d0)*max_leaf_width call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 1) then direction = rnd*360.0d0 elseif (i == 2) then direction = direction + 120 elseif (i == 3) then direction = direction - 240 end if r1 = cos(pi*(angle/180.0d0))*length r2 = sin(pi*(angle/180.0d0))*length x_coord_1 = peti_coord(2, i, 1) y_coord_1 = peti_coord(2, i, 2) z_coord_1 = peti_coord(2, i, 3) x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0)) y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0)) z_coord_1 = z_coord_1 + r1 leaf_coord(2, i, 3, 1) = x_coord_1 leaf_coord(2, i, 3, 2) = y_coord_1 leaf_coord(2, i, 3, 3) = z_coord_1 leaf_coord(2, i, 1, 1) = peti_coord(2, i, 1) leaf_coord(2, i, 1, 2) = peti_coord(2, i, 2) leaf_coord(2, i, 1, 3) = peti_coord(2, i, 3) ! a(1) = leaf_coord(2, i, 3, 1) - peti_coord(2, i, 1) a(2) = leaf_coord(2, i, 3, 2) - peti_coord(2, i, 2) a(3) = leaf_coord(2, i, 3, 3) - peti_coord(2, i, 3) z(1) = 0.0d0 z(2) = 0.0d0 z(3) = 1.0d0 a_unit(:) = a(:)/dsqrt(dot_product(a, a)) b(:) = z(:) - dot_product(a, z)*a_unit(:) b_unit(:) = b(:)/dsqrt(dot_product(b, b)) nvec(1) = a_unit(2)*b_unit(3) - a_unit(3)*b_unit(2) nvec(2) = a_unit(3)*b_unit(1) - a_unit(1)*b_unit(3) nvec(3) = a_unit(1)*b_unit(2) - a_unit(2)*b_unit(1) nvec(:) = nvec(:)/dsqrt(dot_product(nvec, nvec)) mid_length = 1.0d0/(1.0d0 + leaf_shape_rate)*length j = 2 x_coord_1 = peti_coord(2, i, 1) y_coord_1 = peti_coord(2, i, 2) z_coord_1 = peti_coord(2, i, 3) x_coord_1 = x_coord_1 + mid_length*a_unit(1) y_coord_1 = y_coord_1 + mid_length*a_unit(2) z_coord_1 = z_coord_1 + mid_length*a_unit(3) leaf_coord(2, i, j, 1) = x_coord_1 - 0.50d0*width*nvec(1) leaf_coord(2, i, j, 2) = y_coord_1 - 0.50d0*width*nvec(2) leaf_coord(2, i, j, 3) = z_coord_1 - 0.50d0*width*nvec(3) j = 4 x_coord_1 = peti_coord(2, i, 1) y_coord_1 = peti_coord(2, i, 2) z_coord_1 = peti_coord(2, i, 3) x_coord_1 = x_coord_1 + mid_length*a_unit(1) y_coord_1 = y_coord_1 + mid_length*a_unit(2) z_coord_1 = z_coord_1 + mid_length*a_unit(3) leaf_coord(2, i, j, 1) = x_coord_1 + 0.50d0*width*nvec(1) leaf_coord(2, i, j, 2) = y_coord_1 + 0.50d0*width*nvec(2) leaf_coord(2, i, j, 3) = z_coord_1 + 0.50d0*width*nvec(3) end do end subroutine initialize_gw3 !=================================== subroutine growth_gw3(node, d_connect, d_node, max_angle, max_length, & Apical_dominance_p, max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate, & peti_coord, leaf_coord, peti_length) integer, allocatable::d_connect_es(:, :) integer, allocatable, intent(inout)::d_connect(:, :) real(8), allocatable::d_node_es(:, :), peti_coord_es(:, :, :), & leaf_coord_es(:, :, :, :) real(8), allocatable, intent(inout)::d_node(:, :), peti_coord(:, :, :), & leaf_coord(:, :, :, :) real(8), intent(in)::max_angle, max_length, peti_length real(8), intent(in)::max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate real(8) angle, length, direction, r1, r2, pi, x_coord, y_coord, z_coord integer, intent(in)::node, Apical_dominance_p real(8) x_coord_1, y_coord_1, z_coord_1, peti_RUE, mid_length, width real(8), allocatable::a(:), b(:), z(:), nvec(:), a_unit(:), b_unit(:) integer i, j, n, m, o, parent, itr integer :: seedsize integer, allocatable :: seed(:) real :: rnd integer :: c !時間を入れる pi = 3.14159265350d0 allocate (a(3), b(3), z(3), nvec(3), a_unit(3), b_unit(3)) n = size(d_connect, 1) m = size(d_connect, 2) o = size(d_node, 2) allocate (d_connect_es(n + 1, m), d_node_es(n + 1, o)) d_connect_es(:, :) = 0 allocate (peti_coord_es(n, 4, 3), leaf_coord_es(n, 3, 4, 3)) !copy do i = 1, n d_connect_es(i, :) = d_connect(i, :) d_node_es(i, :) = d_node(i, :) end do !generate new node d_connect_es(n + 1, 1) = node if (d_connect_es(node, 2) == -1) then !basis return elseif (d_connect_es(node, 2) == 0) then !aptical node => Get aptical dominance d_connect_es(n + 1, 2) = 0 d_connect_es(node, 2) = 1 parent = d_connect_es(node, 1) itr = 1 do itr = itr + 1 if (d_connect_es(parent, 2) == -1) then !base exit else if (d_connect_es(parent, 2) + 1 >= itr) then d_connect_es(parent, 2) = d_connect_es(parent, 2) + 1 parent = d_connect_es(parent, 1) else cycle end if end if end do elseif (d_connect_es(node, 2) >= Apical_dominance_p) then !generate new aptical d_connect_es(n + 1, 2) = 0 parent = d_connect_es(node, 1) itr = 1 do itr = itr + 1 if (d_connect_es(parent, 2) == -1) then !base exit else if (d_connect_es(parent, 2) + 1 == itr) then d_connect_es(parent, 2) = d_connect_es(parent, 2) + 1 parent = d_connect_es(parent, 1) else exit end if end if end do else !not generate return end if d_connect_es(n + 1, 3) = 0 !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) length = ((rnd - 0.50d0) + 1.0d0)*max_length call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) direction = rnd*360.0d0 r1 = abs(cos(pi*(angle/180.0d0)))*length r2 = abs(sin(pi*(angle/180.0d0)))*length x_coord = d_node(node, 1) y_coord = d_node(node, 2) z_coord = d_node(node, 3) x_coord = x_coord + r2*cos(pi*(direction/180.0d0)) y_coord = y_coord + r2*sin(pi*(direction/180.0d0)) z_coord = z_coord + r1 d_node_es(n + 1, 1) = x_coord d_node_es(n + 1, 2) = y_coord d_node_es(n + 1, 3) = z_coord d_node_es(n + 1, 4) = 1.0d0 d_node_es(n + 1, 5) = 0.0d0 deallocate (d_connect, d_node) n = size(d_connect_es, 1) m = size(d_connect_es, 2) o = size(d_node_es, 2) allocate (d_connect(n, m), d_node(n, o)) d_connect(:, :) = d_connect_es(:, :) d_node(:, :) = d_node_es(:, :) deallocate (d_connect_es, d_node_es) !==================================== !generate peti and leaf !======================== !expand peti_coord,leaf_coord !copy n = size(d_connect, 1) leaf_coord_es(:, :, :, :) = leaf_coord(:, :, :, :) peti_coord_es(:, :, :) = peti_coord(:, :, :) deallocate (peti_coord, leaf_coord) allocate (peti_coord(n, 4, 3), leaf_coord(n, 3, 4, 3)) do i = 1, size(peti_coord, 1) peti_coord(i, :, :) = peti_coord_es(i, :, :) leaf_coord(i, :, :, :) = leaf_coord_es(i, :, :, :) end do deallocate (peti_coord_es, leaf_coord_es) n = size(d_connect, 1) !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) length = ((rnd - 0.50d0) + 1.0d0)*peti_length call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) direction = rnd*360.0d0 r1 = cos(pi*(angle/180.0d0))*length r2 = sin(pi*(angle/180.0d0))*length x_coord = d_node(n, 1) y_coord = d_node(n, 2) z_coord = d_node(n, 3) x_coord = x_coord + r2*cos(pi*(direction/180.0d0)) y_coord = y_coord + r2*sin(pi*(direction/180.0d0)) z_coord = z_coord + r1 peti_coord(1, :, :) = 0.0d0 peti_coord(n, 1, 1) = x_coord peti_coord(n, 1, 2) = y_coord peti_coord(n, 1, 3) = z_coord !second-fourth do i = 2, 4 !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 2) then peti_RUE = 0.10d0 else peti_RUE = 0.050d0 end if length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length*peti_RUE call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 2) then direction = rnd*360.0d0 elseif (i == 3) then direction = direction + 120 elseif (i == 4) then direction = direction - 240 end if r1 = abs(cos(pi*(angle/180.0d0)))*length r2 = abs(sin(pi*(angle/180.0d0)))*length x_coord_1 = x_coord y_coord_1 = y_coord z_coord_1 = z_coord x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0)) y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0)) z_coord_1 = z_coord_1 + r1 peti_coord(n, i, 1) = x_coord_1 peti_coord(n, i, 2) = y_coord_1 peti_coord(n, i, 3) = z_coord_1 end do !leaf generate do i = 1, 3 !random number generation call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) angle = ((rnd - 0.50d0) + 1.0d0)*max_leaf_angle call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) length = ((rnd - 0.50d0) + 1.0d0)*max_leaf_length width = ((rnd - 0.50d0) + 1.0d0)*max_leaf_width call system_clock(count=c) !時間を取得 call random_seed(size=seedsize) allocate (seed(seedsize)) call random_seed(get=seed) seed = c !時間を全部に代入 call random_number(rnd) !rndに乱数をセット deallocate (seed) if (i == 1) then direction = rnd*360.0d0 elseif (i == 2) then direction = direction + 120 elseif (i == 3) then direction = direction - 240 end if r1 = abs(cos(pi*(angle/180.0d0)))*length r2 = abs(sin(pi*(angle/180.0d0)))*length x_coord_1 = peti_coord(n, i, 1) y_coord_1 = peti_coord(n, i, 2) z_coord_1 = peti_coord(n, i, 3) x_coord_1 = x_coord_1 + r2*cos(pi*(direction/180.0d0)) y_coord_1 = y_coord_1 + r2*sin(pi*(direction/180.0d0)) z_coord_1 = z_coord_1 + r1 leaf_coord(n, i, 3, 1) = x_coord_1 leaf_coord(n, i, 3, 2) = y_coord_1 leaf_coord(n, i, 3, 3) = z_coord_1 leaf_coord(n, i, 1, 1) = peti_coord(n, i, 1) leaf_coord(n, i, 1, 2) = peti_coord(n, i, 2) leaf_coord(n, i, 1, 3) = peti_coord(n, i, 3) ! a(1) = leaf_coord(n, i, 3, 1) - peti_coord(n, i, 1) a(2) = leaf_coord(n, i, 3, 2) - peti_coord(n, i, 2) a(3) = leaf_coord(n, i, 3, 3) - peti_coord(n, i, 3) z(1) = 0.0d0 z(2) = 0.0d0 z(3) = 1.0d0 a_unit(:) = a(:)/dsqrt(dot_product(a, a)) b(:) = z(:) - dot_product(a, z)*a_unit(:) b_unit(:) = b(:)/dsqrt(dot_product(b, b)) nvec(1) = a_unit(2)*b_unit(3) - a_unit(3)*b_unit(2) nvec(2) = a_unit(3)*b_unit(1) - a_unit(1)*b_unit(3) nvec(3) = a_unit(1)*b_unit(2) - a_unit(2)*b_unit(1) nvec(:) = nvec(:)/dsqrt(dot_product(nvec, nvec)) mid_length = 1.0d0/(1.0d0 + leaf_shape_rate)*length j = 2 x_coord_1 = peti_coord(n, i, 1) y_coord_1 = peti_coord(n, i, 2) z_coord_1 = peti_coord(n, i, 3) x_coord_1 = x_coord_1 + mid_length*a_unit(1) y_coord_1 = y_coord_1 + mid_length*a_unit(2) z_coord_1 = z_coord_1 + mid_length*a_unit(3) leaf_coord(n, i, j, 1) = x_coord_1 - 0.50d0*width*nvec(1) leaf_coord(n, i, j, 2) = y_coord_1 - 0.50d0*width*nvec(2) leaf_coord(n, i, j, 3) = z_coord_1 - 0.50d0*width*nvec(3) j = 4 x_coord_1 = peti_coord(n, i, 1) y_coord_1 = peti_coord(n, i, 2) z_coord_1 = peti_coord(n, i, 3) x_coord_1 = x_coord_1 + mid_length*a_unit(1) y_coord_1 = y_coord_1 + mid_length*a_unit(2) z_coord_1 = z_coord_1 + mid_length*a_unit(3) leaf_coord(n, i, j, 1) = x_coord_1 + 0.50d0*width*nvec(1) leaf_coord(n, i, j, 2) = y_coord_1 + 0.50d0*width*nvec(2) leaf_coord(n, i, j, 3) = z_coord_1 + 0.50d0*width*nvec(3) end do end subroutine growth_gw3 !=================================== subroutine update_factor_gw3(d_connect, d_node, RUE, num_node, leaf_coord & , day, daylength, limit_leaf_day, day_ID) real(8), intent(inout)::d_node(:, :) real(8), intent(in)::RUE, leaf_coord(:, :, :, :), daylength(:) real(8), allocatable::photosynthesis(:, :), radiation(:, :) real(8) factor, light_area, total_area, old_RUEmeter integer, intent(in)::d_connect(:, :), num_node, day_ID, day(:, :), limit_leaf_day integer i, j, month, node_ID, leaf_ID allocate (photosynthesis(num_node, 3), radiation(12, 2)) photosynthesis(:, :) = 0.0d0 open (60, file="radiation.txt") do i = 1, 12 read (60, *) radiation(i, 1:2) end do close (60) !生成してから1日経ってから光合成開始 !日照面積判定 do i = 1, num_node do j = 1, size(leaf_coord, 2) !compute light area node_ID = i leaf_ID = j call detect_light_area(leaf_coord, light_area, total_area, node_ID, leaf_ID) photosynthesis(i, 1) = light_area photosynthesis(i, 2) = total_area !光合成量＝日射時間(hr)×面積当たり時間当たり日射量(MJ/m^2/hr)×日照面積(m^2)×日射量あたり光合成(g/MJ) month = day(day_ID, 2) if (radiation(month, 1) == 0.0d0) then photosynthesis(i, 3) = 0.0d0 end if if (d_connect(i, 3) == -1) then old_RUEmeter = 0.0d0 else old_RUEmeter = dble(d_connect(i, 3))/dble(limit_leaf_day) if (old_RUEmeter < 0) stop"debugb" end if photosynthesis(i, 3) = daylength(day_ID)*radiation(month, 2)/radiation(month, 1) & *light_area*RUE*old_RUEmeter end do end do do i = 1, num_node !compute factor factor = d_node(i, 5) factor = factor + photosynthesis(i, 3) d_node(i, 5) = factor end do end subroutine update_factor_gw3 !=================================== subroutine detect_start(start, finish, seedling_day, harvesting_day, day) integer, intent(out)::start, finish integer, intent(in)::seedling_day(:), harvesting_day(:), day(:, :) integer i, j, n, yy, mm, dd, diff do i = 1, size(day, 1) yy = day(i, 1) mm = day(i, 2) dd = day(i, 3) diff = abs(seedling_day(1) - yy) & + abs(seedling_day(2) - mm) & + abs(seedling_day(3) - dd) if (diff == 0) then start = i exit end if end do do i = 1, size(day, 1) yy = day(i, 1) mm = day(i, 2) dd = day(i, 3) diff = abs(harvesting_day(1) - yy) & + abs(harvesting_day(2) - mm) & + abs(harvesting_day(3) - dd) if (diff == 0) then finish = i exit end if end do end subroutine detect_start !=================================== subroutine count_br_per_node(node_number, br_per_node, d_connect) integer, intent(in)::node_number, d_connect(:, :) integer, intent(out)::br_per_node integer i br_per_node = 0 do i = 1, size(d_connect, 1) if (node_number == d_connect(i, 1)) then br_per_node = br_per_node + 1 else cycle end if end do if (br_per_node >= 3) stop"debug" end subroutine count_br_per_node !=================================== subroutine detect_light_area(leaf_coord, light_area, total_area, node_ID, leaf_ID) real(8), intent(in)::leaf_coord(:, :, :, :) real(8), intent(out)::light_area, total_area real(8), allocatable::mid_x(:), shape_function(:), evaluate_point(:, :), eva_x(:), rvec(:) & , x1(:), x2(:), y1(:), y2(:) real(8) gzi_1, gzi_2, eva_r, mid_r, both_r, l1, l2, l3 integer, intent(in)::node_ID, leaf_ID integer i, j, k, l, m, shade_or_not, count_shade allocate (mid_x(3), shape_function(4), evaluate_point(4, 2), eva_x(3), rvec(2), x1(2), x2(2)) allocate (y1(3), y2(3)) !2D linear interpolation,4 points evaluate_point(1, 1) = -0.50d0 evaluate_point(1, 2) = -0.50d0 evaluate_point(2, 1) = -0.50d0 evaluate_point(2, 2) = 0.50d0 evaluate_point(3, 1) = 0.50d0 evaluate_point(3, 2) = 0.50d0 evaluate_point(4, 1) = 0.50d0 evaluate_point(4, 2) = -0.50d0 y1(1:3) = leaf_coord(node_ID, leaf_ID, 1, 1:3) - leaf_coord(node_ID, leaf_ID, 3, 1:3) y2(1:3) = leaf_coord(node_ID, leaf_ID, 2, 1:3) - leaf_coord(node_ID, leaf_ID, 4, 1:3) l1 = dsqrt(dot_product(y1, y1)) l2 = dsqrt(dot_product(y2, y2)) total_area = 0.50d0*l1*l2/100.0d0/100.0d0 count_shade = 0 do k = 1, 4 !evaluate points gzi_1 = evaluate_point(k, 1) gzi_2 = evaluate_point(k, 2) shape_function(1) = 0.250d0*(1.0d0 - gzi_1)*(1.0d0 - gzi_1) shape_function(2) = 0.250d0*(1.0d0 - gzi_1)*(1.0d0 + gzi_1) shape_function(3) = 0.250d0*(1.0d0 + gzi_1)*(1.0d0 + gzi_1) shape_function(4) = 0.250d0*(1.0d0 + gzi_1)*(1.0d0 - gzi_1) eva_x(:) = 0.0d0 do l = 1, 4 eva_x(:) = eva_x(:) + leaf_coord(node_ID, leaf_ID, l, :) & *shape_function(l) end do rvec(:) = leaf_coord(node_ID, leaf_ID, k, 1:2) eva_r = dsqrt(dot_product(rvec - eva_x(1:2), rvec - eva_x(1:2))) shade_or_not = 0 do i = 1, size(leaf_coord, 1) do j = 1, size(leaf_coord, 2) mid_x(:) = 0.0d0 do l = 1, 4 mid_x(:) = mid_x(:) + 0.250d0*leaf_coord(i, j, l, :) end do !if z coordinate of the midpoint of the leaf > evaluate point ! => possible node if (mid_x(3) < eva_x(3)) then cycle end if x1(1:2) = leaf_coord(i, j, 1, 1:2) - leaf_coord(i, j, 3, 1:2) x2(1:2) = leaf_coord(i, j, 2, 1:2) - leaf_coord(i, j, 4, 1:2) mid_r = dsqrt(dot_product(x1, x1)) + dsqrt(dot_product(x2, x2)) mid_r = 0.50d0*mid_r both_r = dsqrt(dot_product(mid_x(1:2) - eva_x(1:2), mid_x(1:2) - eva_x(1:2))) if (both_r > eva_r + mid_r) then cycle else shade_or_not = 1 exit end if end do if (shade_or_not == 1) then exit end if end do if (shade_or_not == 1) then cycle else count_shade = count_shade + 1 end if end do light_area = total_area/4.0d0*dble(count_shade) end subroutine detect_light_area !=================================== subroutine get_old(d_connect, limit_leaf_day) integer, intent(inout)::d_connect(:, :) integer, intent(in)::limit_leaf_day integer i do i = 1, size(d_connect, 1) if (d_connect(i, 3) == 0) then !first day d_connect(i, 3) = limit_leaf_day elseif (d_connect(i, 3) == 1) then d_connect(i, 3) = -1 elseif (d_connect(i, 3) == -1) then d_connect(i, 3) = -1 elseif (limit_leaf_day >= d_connect(i, 3) .and. d_connect(i, 3) >= 2) then d_connect(i, 3) = d_connect(i, 3) - 1 else stop"invalid d_connect,3" end if end do end subroutine get_old !=================================== subroutine generate_pod(d_node, d_pod, source_convection_rate, source_limitation) real(8), intent(inout)::d_node(:, :), d_pod(:, :) real(8), intent(in)::source_convection_rate, source_limitation real(8) total_source, extra_source integer i, j, pn do i = 1, size(d_node, 1) total_source = d_node(i, 5) pn = 0 do j = 1, size(d_pod, 2) if (d_pod(i, j) == source_limitation) then cycle elseif (d_pod(i, j) < source_limitation) then pn = j exit else stop"invalid source/pod" end if end do if (pn == 0) then !sink was occupied cycle end if if (total_source >= source_convection_rate) then d_node(i, 5) = total_source - source_convection_rate d_pod(i, pn) = d_pod(i, pn) + source_convection_rate elseif (source_convection_rate > total_source .and. & total_source > 0.0d0) then d_pod(i, pn) = d_pod(i, pn) + total_source d_node(i, 5) = 0.0d0 elseif (total_source == 0.0d0) then cycle else stop"invalid source < 0" end if extra_source = 0.0d0 do j = 1, size(d_pod, 2) if (extra_source /= 0.0d0) then d_pod(i, j) = d_pod(i, j) + extra_source end if if (d_pod(i, j) > source_limitation) then extra_source = d_pod(i, j) - source_limitation d_pod(i, j) = source_limitation end if end do if (extra_source /= 0.0d0) then d_node(i, 5) = d_node(i, 5) + extra_source extra_source = 0.0d0 end if end do end subroutine generate_pod !=================================== end module growthmod