四阶龙格库塔法程序—— FORTRAN语言编写

代码:

SUBROUTINE runge_kutta()
!关于Runge-Kutta方法，该方法是用来解形如y'=f(t,y)的常微分方程的
!经典的4阶R-K方法的公式如下：
! Yn+1 = Yn + h/6 * (K1+2K2+2K3+K4)
!其中
! K1=f(Tn,Yn)
! K2=f(Tn+h/2,Yn+h/2*K1)
! K3=f(Tn+h/2,Yn+h/2*K2)
! K4=f(Tn+h,Yn+h*K3)
!
type(node_fish),pointer :: p
!
integer :: i,j
!
!-----STEP 1 ：保存Yn
do i=1,Num_lon
do j=1,Num_lat
Yn%nutrients(i,j) = patches(i,j)%core%nutrients
Yn%planktons(i,j) = patches(i,j)%core%planktons
end do
end do
!
if (anchovy%Num_node>Lim_fish) call abort('ERROR(dt01):too many fishes')
!
p=>anchovy%head%next
do i=1,anchovy%Num_node
Yn%fish(i)=p%value%weight
p=>p%next
end do
!
!-----STEP 2 ：依次计算K1-K4
!-------STEP 2.1 ：求K1，时间是t，Yn不变
K1=integral_function(rt_current,Yn)
!-------STEP 2.2 ：求K2，时间更新到t+h/2，Yn更新到Yn+(h/2)*K1
K2=integral_function(rt_current+rt_step/2.0,Yn+(rt_step/2.0)*K1)
!-------STEP 2.3 ：求K3，时间更新到t+h/2，Yn更新到Yn+(h/2)*K2
K3=integral_function(rt_current+rt_step/2.0,Yn+(rt_step/2.0)*K2)
!-------STEP 2.4 ：求K4，时间更新到t+h，Yn更新到Yn+h*K3
K4=integral_function(rt_current+rt_step,Yn+rt_step*K3)
!
!-----STEP 3 ：由Yn和K1-K4计算Yn+1
Yn_1 = Yn + (rt_step/6.0)*(K1+2.0*K2+2.0*K3+K4)
!
!-----STEP 4 ：赋值
do i=1,Num_lon
do j=1,Num_lat
if (Yn_1%nutrients(i,j)%NO3 < min_nutrient%NO3 ) Yn_1%nutrients(i,j)%NO3 = min_nutrient%NO3
if (Yn_1%nutrients(i,j)%NH4 < min_nutrient%NH4 ) Yn_1%nutrients(i,j)%NH4 = min_nutrient%NH4
if (Yn_1%nutrients(i,j)%PON < min_nutrient%PON ) Yn_1%nutrients(i,j)%PON = min_nutrient%PON
if (Yn_1%nutrients(i,j)%DON < min_nutrient%DON ) Yn_1%nutrients(i,j)%DON = min_nutrient%DON
if (Yn_1%nutrients(i,j)%SiOH4 < min_nutrient%SiOH4) Yn_1%nutrients(i,j)%SiOH4 = min_nutrient%SiOH4
if (Yn_1%nutrients(i,j)%Opal < min_nutrient%Opal ) Yn_1%nutrients(i,j)%Opal = min_nutrient%Opal
if (Yn_1%nutrients(i,j)%HPO4 < min_nutrient%HPO4 ) Yn_1%nutrients(i,j)%HPO4 = min_nutrient%HPO4
if (Yn_1%nutrients(i,j)%POP < min_nutrient%POP ) Yn_1%nutrients(i,j)%POP = min_nutrient%POP
if (Yn_1%nutrients(i,j)%DOP < min_nutrient%DOP ) Yn_1%nutrients(i,j)%DOP = min_nutrient%DOP
patches(i,j)%core%nutrients = Yn_1%nutrients(i,j)
if (Yn_1%planktons(i,j)%phyto_small < min_plankton%phyto_small ) Yn_1%planktons(i,j)%phyto_small = min_plankton%phyto_small
if (Yn_1%planktons(i,j)%phyto_large < min_plankton%phyto_large ) Yn_1%planktons(i,j)%phyto_large = min_plankton%phyto_large
if (Yn_1%planktons(i,j)%zoo_small < min_plankton%zoo_small ) Yn_1%planktons(i,j)%zoo_small = min_plankton%zoo_small
if (Yn_1%planktons(i,j)%zoo_large < min_plankton%zoo_large
) Yn_1%planktons(i,j)%zoo_large

= min_plankton%zoo_large
if (Yn_1%planktons(i,j)%zoo_predatory < min_plankton%zoo_predatory ) Yn_1%planktons(i,j)%zoo_predatory = min_plankton%zoo_predatory
patches(i,j)%core%planktons = Yn_1%planktons(i,j)
end do
end do
!
p=>anchovy%head%next
do i=1,anchovy%Num_node
p%value%weight = Yn_1%fish(i)
p%value%growth = Yn_1%fish(i) - Yn%fish(i)
call update_half_gene(p%value,rt_current+rt_step)
call update_fish_state(p%value,rt_current+rt_step)
p=>p%next
end do
!
END SUBROUTINE runge_kutta