Фортран подпрограмма не возвращает значения

Question

У меня есть эта программа -

      PROGRAM main
          real*8 :: P1,t,m
          t=0.d1
          P1=0.d0
          open(1,file="popr.txt")
          do while(t.lt.11.d3)
              t=t+2.d1
              CALL Mod_PoPs(1,0.d0,t,2.d3,10.d3,2.d5,2.d3,3,P1) ! most of the values passed here are junk or used for something not shown in the part of the subroutine here
              !write (1,*) t,P1
              write (*,*) t,P1
          enddo
          close(1)
      END PROGRAM

Это всего лишь тестовый код, поэтому я не очень беспокоюсь о надлежащих стандартах (я все время путаюсь с ними, поэтому я только исправляю все на последних этапах). Однако я также попытался поместить подпрограмму в реальный код, соответствующий стандартам 2008 года, но она все еще не работала. Соответствующая часть подпрограммы Mod_PoPs такова -

      SUBROUTINE Mod_PoPs(opt,Pops0,t,t_st,t_end,per,inte,n,Pops)
          integer, intent(in) :: opt,n
          real*8, intent(in) :: per,Pops0,t_st,t_end,t,inte
          real*8, intent(out) :: Pops
          if (opt.eq.1) then
              if(t.eq.t_st) then
                  Pops=per+Pops0
                  write (*,*) "change Pops here: ", Pops
              else
                  Pops=Pops0
              endif
          endif
      END

А вот соответствующая часть вывода -

 1980.00000000000       0.000000000000000E+000
 change Pops here:    200000.000000000     
 2000.00000000000       0.000000000000000E+000

Как видите, значение P1 должно измените на 200000, когда t = 2000, и когда вызывается подпрограмма, она явно переходит в правильное состояние при правильном t, но каким-то образом это изменение не отражается в выходных данных. Я не смог понять это. Я пытался использовать gfortran и ifort, но ничего не получалось. Это имеет еще меньше смысла, потому что в подпрограмме есть другие случаи (когда opt! = 1), и во всех них Pops успешно изменяется. Я не включил другие случаи для ясности, но вот пример другого случая, который прекрасно работает -

          if (opt.eq.9) then
              if (t.ge.t_st.and.t.le.t_e) then
                  Pops=Pops0+per
              else
                  Pops=Pops0
              endif
          endif

Я подумал, что это может быть проблемой с плавающим сравнением, особенно равенством. Но есть и другие случаи, в которых есть равенство, но они работают нормально. Я опубликую всю подпрограмму, если потребуется. Я также пытался не использовать равенство, как это -

      if (opt.eq.1) then
          if(abs(t-t_st).lt.2.d1) then
              Pops=per+Pops0
              write (*,*) "change Pops here: ", Pops
          else
              Pops=Pops0
          endif
      endif

Все еще не работает.

Я смотрел на другие вопросы с похожими названиями, но ничего не помогло. Я подумал, что это может быть крайний случай, определяющий c для моей системы и не воспроизводимый, поэтому я попытался запустить тот же код на моем сервере, но тот же результат и там. То, что я хочу знать, это

1. , если эта ошибка даже воспроизводима
2. как мне ее исправить