problem with subroutines input/output - application: ODE solver

[1ramirez1]

Dear Forum,

I tried to use a Runge-Kutta method to solve a system of three ODEs (dNdt, dPdt, dZdt).
The code is splitted in the main program, the ODE subroutine and a solver subroutine.
I am very new to Fortran (90), therefore my problem could be very basic and resulting from misunderstanding of basic fortran concepts, even though I read quite a few Introductions explaining the exchange of values between the main program and subroutines.
My problem is that the calculated values of the state variables (N,P,Z) in the ODE and the solver subroutines are not returned to the main program and hence the state variables stay constant over the time loop.
Attached is the basic code outline:

Thank you very much in advance for your help!

PROGRAM NPZ
  
  USE solver 

  (Declarations...)

  (Initial Cond.)

  ! combine all state variables:
  y(1) = N
  y(2) = P
  y(3) = Z


   DO k = 1,tend,dt                   
            
 
      CALL rk4(y,dydt,yout)     ! Call the solver subroutine
    
     N = y(1)
     P = y(2)
     Z = y(3)
     
  END DO

END PROGRAM




MODULE solver
       CONTAINS


      
      SUBROUTINE ode(y,f)      
      
      (Declarations)


      ! De-combine state variables
      N = y(1)
      P = y(2)
      Z = y(3)

      dNdt = ...
      dPdt = ...
      dZdt = ...

      f(1) = dNdt
      f(2) = dPdt 
      f(3) = dZdt

      dydt = f


      RETURN
    END SUBROUTINE ode

     





    SUBROUTINE rk4(y,dydt,yout)  

      (Declarations)
 
      ! De-combine state variables
      N = y1(1)
      P = y1(2)
      Z = y1(3)


    
      CALL ode(y,dydt)      
      DO i=1,imax
         s1(i) = dydt(i)   
      END DO


      DO i=1,imax
         ytemp1(i) = y(i) + (dt/2.)*s1(i)
      END DO


      CALL ode(ytemp1,dydt)
      DO i=1,imax
         s2(i) = dydt(i)   
      END DO


      DO i=1,imax
         ytemp2(i) = y(i)+(dt/2.)*s2(i)
      END DO

      CALL ode(ytemp2,dydt)
      DO i=1,imax
         s3(i) = dydt(i)
      END DO



      DO i=1,imax
         ytemp3(i) = y(i)+dt*s3(i)
      END DO


      CALL ode(ytemp3,dydt)
      DO i=1,imax
         s4(i) = dydt(i)
      END DO


       DO i=1,imax
        yout(i) = y(i) + (dt/6.)*(s1(i) + 2*s2(i) + 2*s3(i) + s4(i)) 
        y(i) = yout(i)
        dydt(i) = s1(i)   
       END DO
       
 

      RETURN
    END SUBROUTINE rk4




END MODULE solver

 

[1ramirez1]

Another problem is that a namelist that contains all the parameters is not global available to the subroutines.

I define and read the namelist in the main program as follows:

REAL :: alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero
INTEGER :: dt, tend

NAMELIST /NPZNAME/ alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, dt, tend

OPEN(UNIT=41,FILE='npz.nml',STATUS='OLD')
READ(41,NML=NPZNAME)
! Don't close here !

What statement do I have to use in the subroutines that they can access those parameters?
Thanks!