计算过程中数值不稳定性 有偿

NAN00

计算图灵斑样貌过程中，使用如下代码进行多次迭代，发现在控制相同输入、参数情况下，有时会得到不同的结果，已经排查不是并行导致问题，希望老师们帮忙看看
subroutine myosin_reaction_diffusion2(Nc,nb_cell,num_nb,Dam,Dcr,am_a,cr_r,ka_on,ka_off,kr_off,delta,am_conc,cr_conc)
    implicit none
    integer, intent(in) :: Nc,nb_cell(20,Nc),num_nb(Nc)
    double precision, intent(in) :: Dam,Dcr,am_a,cr_r,ka_on,ka_off,kr_off,delta
    double precision, intent(inout) :: am_conc(Nc), cr_conc(Nc)

    integer :: i,j,i1
    double precision :: temp1_am,temp2_am,temp3_am,temp1_cr,temp2_cr
    double precision :: am_conc_old(Nc), cr_conc_old(Nc),delta01

    am_conc_old = am_conc; cr_conc_old = cr_conc; delta01 = delta * 0.5d0
!$omp parallel private(i,j,i1,temp1_am,temp2_am,temp3_am,temp1_cr,temp2_cr)
!$omp do  
    do i = 1,Nc
        ! temp1_am = 0.d0; temp1_cr = 0.d0
        ! do j = 1,num_nb(i)
        !     i1 = nb_cell(j,i)

        !     temp1_am = temp1_am + ( am_conc_old(i1) - am_conc_old(i) )
        !     temp1_cr = temp1_cr + ( cr_conc_old(i1) - cr_conc_old(i) )
        ! end do
        ! temp1_am = Dam * temp1_am; temp1_cr = Dcr * temp1_cr
        temp1_am = Dam * sum(am_conc_old(nb_cell(1:num_nb(i),i)) - am_conc_old(i))
        temp1_cr = Dcr * sum(cr_conc_old(nb_cell(1:num_nb(i),i)) - cr_conc_old(i))

        temp2_am = ka_on * (am_conc_old(i)**2) * cr_conc_old(i)
        temp2_cr = -kr_off * (am_conc_old(i)**2) * cr_conc_old(i)

        temp3_am = -ka_off * am_conc_old(i)

        am_conc(i) = am_conc_old(i) + (temp1_am + temp2_am + temp3_am + am_a) * delta01

        cr_conc(i) = cr_conc_old(i) + (temp1_cr + temp2_cr + cr_r) * delta01
    end do
!$omp end do
!$omp end parallel

end subroutine

chiangtp

1. please check (num_nb(i)>=1 .AND. num_nb(i)<=20)

   do i = 1, Nc                                                      
     IF( num_nb(i)<1 .OR. num_nb(i)>20 ) WRITE(*,*) num_nb(i)

2. "在控制相同输入、参数情况下，有时会得到不同的结果" ---> (???) "数值不稳定"