module inoput_gw3 use gnuplot implicit none contains !=================================== subroutine input_gw3(restart, max_angle, LDL, max_length, & Apical_dominance_p, seedling_day, harvesting_day, RUE, max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate, limit_leaf_day, & peti_length, max_pod_number, source_convection_rate, source_limitation, & seed_pod_rate, threshold) real(8), intent(out)::max_angle, max_length, RUE, LDL, max_leaf_angle, & max_leaf_length, max_leaf_width, leaf_shape_rate, peti_length, & source_convection_rate, source_limitation, seed_pod_rate, threshold integer, intent(out)::restart, Apical_dominance_p, limit_leaf_day, max_pod_number integer, allocatable, intent(out)::seedling_day(:), harvesting_day(:) character*50 infile print *, "growth simulator 1.0" print *, "input file name" read (*, *) infile open (10, file=infile) allocate (seedling_day(3), harvesting_day(3)) read (10, *) max_angle, max_leaf_angle read (10, *) max_length read (10, *) RUE, limit_leaf_day read (10, *) seedling_day(1), seedling_day(2), seedling_day(3) read (10, *) harvesting_day(1), harvesting_day(2), harvesting_day(3) read (10, *) LDL read (10, *) Apical_dominance_p read (10, *) max_leaf_length, max_leaf_width, peti_length, leaf_shape_rate read (10, *) max_pod_number, source_convection_rate, source_limitation, & seed_pod_rate read (10, *) restart close (10) end subroutine input_gw3 !=================================== subroutine input_daylength(day, daylength) integer, allocatable, intent(out)::day(:, :) real(8), allocatable, intent(out)::daylength(:) integer i, j, n open (10, file="daylength.txt") read (10, *) n allocate (day(n, 3), daylength(n)) do i = 1, n read (10, *) day(i, 1:3), daylength(i)!yy,mm,dd,daylength end do close (10) end subroutine input_daylength !=================================== subroutine gnuplot_gw3(d_connect, d_node, step, day, day_ID) integer, intent(in)::d_connect(:, :), step, day(:, :), day_ID real(8), allocatable::nod_coord(:, :) real(8), intent(in)::d_node(:, :) integer, allocatable::elem_nod(:, :) integer i, j, n n = size(d_node, 1) allocate (nod_coord(n, 3), elem_nod(n, 2)) do i = 1, n nod_coord(i, 1:3) = d_node(i, 1:3) elem_nod(i, 1) = i elem_nod(i, 2) = d_connect(i, 1) end do call gnuplot_out(elem_nod, nod_coord, step, day, day_ID) end subroutine gnuplot_gw3 !=================================== subroutine gnuplot_leaf_gw3(leaf_coord, peti_coord, d_node, step) real(8), intent(in) ::d_node(:, :), leaf_coord(:, :, :, :), peti_coord(:, :, :) integer, intent(in)::step integer i, j character filename1*17 character filename2*17 !======gnuplot動画出力用コマンド=========== !データ作成用プログラム !連番ファイルの作成----------------------- i = step write (filename1, '("gnu_pet", i6.6, ".txt")') i ! ここでファイル名を生成している !write (filename2, '("gnu_pet", i6.6, ".png")') i open (40, file=filename1) !---------------------------------------- !png出力用スクリプトファイル----------------- !open(50,file='png_script.gp',position='append') !write(50,*)'set output ','"',filename2,'"' !write(50,*)'replot ','"',filename1,'" with vectors' !============================================= do i = 1, size(peti_coord, 1) if (d_node(i, 5) == 0.0d0) then write (40, *) d_node(1, 1:3) write (40, *) d_node(1, 1:3) write (40, *) ' ' write (40, *) ' ' cycle end if write (40, *) peti_coord(i, 1, 1:3) write (40, *) d_node(i, 1:3) write (40, *) ' ' write (40, *) ' ' write (40, *) peti_coord(i, 2, 1:3) write (40, *) peti_coord(i, 1, 1:3) write (40, *) ' ' write (40, *) ' ' write (40, *) peti_coord(i, 3, 1:3) write (40, *) peti_coord(i, 1, 1:3) write (40, *) ' ' write (40, *) ' ' write (40, *) peti_coord(i, 4, 1:3) write (40, *) peti_coord(i, 1, 1:3) write (40, *) ' ' write (40, *) ' ' end do close (40) !close(50) i = step write (filename1, '("gnu_lef", i6.6, ".txt")') i ! ここでファイル名を生成している !write (filename2, '("gnu_pet", i6.6, ".png")') i open (40, file=filename1) !---------------------------------------- !png出力用スクリプトファイル----------------- !open(50,file='png_script.gp',position='append') !write(50,*)'set output ','"',filename2,'"' !write(50,*)'replot ','"',filename1,'" with lines' !============================================= do i = 2, size(leaf_coord, 1) if (d_node(i, 5) == 0.0d0) then write (40, *) d_node(1, 1:3) write (40, *) d_node(1, 1:3) write (40, *) ' ' write (40, *) ' ' cycle end if do j = 1, size(leaf_coord, 2) !plot leaf write (40, *) leaf_coord(i, j, 1, 1:3) write (40, *) leaf_coord(i, j, 2, 1:3) write (40, *) leaf_coord(i, j, 3, 1:3) write (40, *) leaf_coord(i, j, 4, 1:3) write (40, *) leaf_coord(i, j, 1, 1:3) write (40, *) ' ' write (40, *) ' ' end do end do close (40) !close(50) end subroutine gnuplot_leaf_gw3 !=================================== subroutine output_gw3(d_connect, d_node) integer, intent(in)::d_connect(:, :) real(8), intent(in)::d_node(:, :) integer i, n open (10, file="output_gw.d") n = size(d_connect, 1) do i = 1, n write (10, *) d_connect(i, :) end do do i = 1, n write (10, *) d_node(i, :) end do close (10) end subroutine output_gw3 !=================================== subroutine output_factor(d_node, leaf_coord, occupied_area, step) real(8), intent(in)::d_node(:, :), leaf_coord(:, :, :, :) real(8), intent(out)::occupied_area real(8), allocatable::a(:), b(:) real(8) total_factor, r, r_tr, pi, x1, x2, l1, l2, LA integer, intent(in)::step integer i, j, k allocate (a(3), b(3)) pi = 3.1415926535d0 if (step == 1) then open (10, file="resource.txt") write (10, *) "days, resource/area, LAI, total resource, area(m^2), leaf area(m^2)" else open (10, file="resource.txt", position='append') end if r = 0.0d0 !compute occupied areaby circle do i = 1, size(leaf_coord, 1) !node ID do j = 1, size(leaf_coord, 2)! Leaf ID do k = 1, size(leaf_coord, 3)! edge ID x1 = leaf_coord(i, j, k, 1) x2 = leaf_coord(i, j, k, 2) r_tr = dsqrt(x1*x1 + x2*x2) if (r < r_tr) then r = r_tr else cycle end if end do end do end do occupied_area = r*r*pi*0.01*0.01 LA = 0.0d0 !compute total leaf area do i = 1, size(leaf_coord, 1) !node ID do j = 1, size(leaf_coord, 2)! Leaf ID a(:) = leaf_coord(i, j, 3, :) - leaf_coord(i, j, 1, :) b(:) = leaf_coord(i, j, 4, :) - leaf_coord(i, j, 2, :) l1 = dsqrt(dot_product(a, a)) l2 = dsqrt(dot_product(b, b)) LA = LA + 0.50d0*l1*l2*0.01*0.01 end do end do total_factor = 0.0d0 do i = 1, size(d_node, 1) total_factor = total_factor + d_node(i, 5) end do write (10, *) step, int(total_factor/occupied_area), LA/occupied_area, & total_factor, occupied_area, LA close (10) end subroutine output_factor !=================================== subroutine output_harvest(d_pod, structural_TDW, d_node, seed_pod_rate, & occupied_area) real(8), intent(in)::d_pod(:, :), structural_TDW, d_node(:, :), & seed_pod_rate, occupied_area real(8) TDW, Harvest_TDW, seed_weight, yield, soyseed_per_dw, pod_weight integer i, j TDW = structural_TDW Harvest_TDW = 0.0d0 do i = 1, size(d_pod, 1) do j = 1, size(d_pod, 2) Harvest_TDW = Harvest_TDW + d_pod(i, j) end do end do soyseed_per_dw = 2.0d0 seed_weight = seed_pod_rate/(1.0d0 + soyseed_per_dw*seed_pod_rate) & *Harvest_TDW pod_weight = 1.0d0/(1.0d0 + soyseed_per_dw*seed_pod_rate) & *Harvest_TDW TDW = TDW + seed_weight + pod_weight do i = 1, size(d_node, 1) TDW = TDW + d_node(i, 5) end do open (10, file="harvest.txt") write (10, *) "seed weight (g) :", seed_weight write (10, *) "Yield (kg/10a) :", seed_weight/occupied_area write (10, *) "Total Dry Weight(g) :", TDW write (10, *) "TDW yield (kg/10a) :", TDW/occupied_area write (10, *) "Harvest Index :", seed_weight/TDW close (10) end subroutine output_harvest !=================================== end module inoput_gw3