!
!sns.f90
!


program main
  implicit none
  
  integer(kind=4),parameter:: neq = 6  !number of equations
  integer(kind=4),parameter:: nmax = 365  !number of step

  !parameters for runge.f90
  real(kind=8),parameter:: eps = 1.0d-10
  external ex2
  integer(kind=4) ierr,i,n,iseed
  real(kind=8),dimension(neq,12):: work
  real(kind=8),dimension(neq):: y0,yn 
  real(kind=8) x0,xn,dx !time

  real(8) up(2)

  !intial conditions
  n=0
  x0=0.0d0           !time
  dx=1.0d0/365.0d0   !timestep, dx=1day

  y0(1) = 100.0d0   !yaruki 1
  y0(2) = 100.0d0   !yaruki 2

  y0(3) = 10.0d0    !aki 1
  y0(4) = 5.0d0     !aki 2
  y0(5) = 0.0d0     !interaction 1
  y0(6) = 0.0d0     !interaction 2
  
  up(1) = 7.0d0    !recover b  1
  up(2) = 7.0d0    !recover b  2

  open(10, file='sns.dat')
  write(10, '(7f10.4)') 365.0d0*x0, y0(1), y0(2), y0(3), y0(4), y0(5), y0(6)
  
  do n=1, nmax
     xn = x0 + dx
     
     call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work)	!runge.f90
     x0 = xn  !time evolution
     
     do i = 1,neq
        y0(i) = yn(i)  !update y0
     end do    
     
     write(10, '(7f10.4)') 365.0d0*x0, y0(1), y0(2), y0(3), y0(4), y0(5), y0(6)   

     if( mod( int(x0*365.0d0), 20)==0 ) then       !how often upload (each 20days)
        y0(5) = -up(1)
        if( y0(1)>100.0d0 ) then
        y0(5) = y0(5) -up(1) + up(1)*  tanh( y0(1)/100.0d0 -1.0d0 )
        end if
     end if

     if( mod( int(x0*365.0d0), 30)==0 ) then      !each 30days
        y0(6) = -up(2)
        if( y0(2)>100.0d0 ) then
        y0(6) = y0(6) -up(2) + up(2)*  tanh( y0(6)/100.0d0 -1.0d0 )
        end if
     end if

    !if yaruki is less than 0, no more sns.
     do i=1, 2
        if(y0(i)<1) then
           goto 001
        end if
     end do
     
  end do
  
001 close(10)  
end program main


!function
real(kind=8) function ex2(i,x,y)
  implicit none  
  
  integer(kind=4),parameter:: neq=6
  integer(kind=4) i
  real(kind=8) x
  real(kind=8),dimension(neq):: y
  
  !equations
  if(i==1) then  !yaruki
     ex2 = -y(3)* y(1) -y(5)* (y(2)/ y(1)) * 100.0d0
  else if (i==2) then  !yaruki
     ex2 = -y(4)* y(2) -y(6)* (y(1)/ y(2)) * 100.0d0
  else if (i==3) then  !aki
     ex2 = 10.0d0* tanh( -1.0d0 + 100.0d0/y(1) )
  else if (i==4) then  !aki
     ex2 = 5.0d0* tanh( -1.0d0 + 100.0d0/y(2) )
  else if (i==5) then   !interaction
     ex2 = 100.0d0* ( 100.0d0/y(1) )* tanh( 1.0d0 - y(1)/100.0d0 )
  else                  !interaction
     ex2 = 100.0d0* ( 100.0d0/y(2) )* tanh( 1.0d0 - y(2)/100.0d0 )
  end if
  
  return
end function ex2


