General ODE solver (Source File: ode_solvers.F90)

INTERFACE:

    subroutine ode_solver(solver,numc,nlev,dt,cc,right_hand_side_rhs,right_hand_side_ppdd)

DESCRIPTION:

Here, 10 different numerical solvers for the right hand sides of the biogeochemical models are implemented for computing the ordinary differential equations (ODEs) which are calculated as the second step of the operational split method for the complete biogeochemical models. The remaining ODE is

$$\partial_t c_i = P_i(\vec{c}) -D_i(\vec{c}), \;\; i = 1,\ldots,I,$$ (247)

with $c_i$ denoting the concentrations of state variables. The right hand side denotes the reaction terms, which are composed of contributions $d_{i,j}(\vec{c})$, which represent reactive fluxes from $c_i$ to $c_j$, and in turn, $p_{i,j}(\vec{c})$ are reactive fluxes from $c_j$ received by $c_i$, see equation (265).

These methods are:

1. First-order explicit (not unconditionally positive)
2. Second order explicit Runge-Kutta (not unconditionally positive)
3. Fourth-order explicit Runge-Kutta (not unconditionally positive)
4. First-order Patankar (not conservative)
5. Second-order Patankar-Runge-Kutta (not conservative)
6. Fourth-order Patankar-Runge-Kutta (does not work, not conservative)
7. First-order Modified Patankar (conservative and positive)
8. Second-order Modified Patankar-Runge-Kutta (conservative and positive)
9. Fourth-order Modified Patankar-Runge-Kutta (does not work, conservative and positive)
10. First-order Extended Modified Patankar (stoichiometrically conservative and positive)
11. Second-order Extended Modified Patankar-Runge-Kutta (stoichiometrically conservative and positive)

The schemes 1 - 5 and 7 - 8 have been described in detail by Burchard et al. (2003). Later, Bruggeman et al. (2006) could show that the Modified Patankar schemes 7 - 8 are only conservative for one limiting nutrient and therefore they developed the Extended Modified Patankar (EMP) schemes 10 and 11 which are also stoichiometrically conservative. Patankar and Modified Patankar schemes of fourth order have not yet been developed, such that choices 6 and 9 do not work yet.

The call to ode_solver() requires a little explanation. The first argument solver is an integer and specifies which solver to use. The arguments numc and nlev gives the dimensions of the data structure cc i.e. cc(1:numc,0:nlev). dt is simply the time step. The last argument is the most complicated. right_hand_side is a subroutine with a fixed argument list. The subroutine evaluates the right hand side of the ODE and may be called more than once during one time-step - for higher order schemes. For an example of a correct right_hand_side have a look at e.g. do_bio_npzd()

USES:

    IMPLICIT NONE

INPUT PARAMETERS:

    integer, intent(in)                 :: solver,nlev,numc
    REALTYPE, intent(in)                :: dt
   !INPUT/OUTPUT PARAMETER:
    REALTYPE, intent(inout)             :: cc(1:numc,0:nlev)
 
    interface
       subroutine right_hand_side_ppdd(first,numc,nlev,cc,pp,dd)
          logical, intent(in)                  :: first
          integer, intent(in)                  :: numc,nlev
          REALTYPE, intent(in)                 :: cc(1:numc,0:nlev)
 
          REALTYPE, intent(out)                :: pp(1:numc,1:numc,0:nlev)
          REALTYPE, intent(out)                :: dd(1:numc,1:numc,0:nlev)
       end
    end interface
 
    interface
       subroutine right_hand_side_rhs(first,numc,nlev,cc,rhs)
          logical, intent(in)                  :: first
          integer, intent(in)                  :: numc,nlev
          REALTYPE, intent(in)                 :: cc(1:numc,0:nlev)
          REALTYPE, intent(out)                :: rhs(1:numc,0:nlev)
       end
    end interface

REVISION HISTORY:

    Original author(s): Hans Burchard, Karsten Bolding

LOCAL VARIABLES: