Second-order Extended Modified Patankar scheme


INTERFACE:

    subroutine emp_2(dt,numc,nlev,cc,right_hand_side)
DESCRIPTION:

Here, the second-order Extended Modified Patankar scheme for biogeochemical models is coded, with two evaluations of the right-hand side per time step:


$\displaystyle \vec{c}^{(1)}$ $\displaystyle =$ $\displaystyle \vec{c}^n + \Delta t \: \vec{f}(t^n ,\vec{c}^n )\prod\limits_{j \in J^n } {\frac{{c_j^{(1)} }}{{c_j^n }}}$  
$\displaystyle \vec{c}^{n + 1}$ $\displaystyle =$ $\displaystyle \vec{c}^n + \frac{{\Delta t}}{2}\left( {\vec{f}(t^n ,\vec{c}^n ) ...
...{(1)} )} \right)\prod\limits_{k \in K^n } {\frac{{c_k^{n + 1} }}{{c_k^{(1)} }}}$ (256)

where


$\displaystyle J^n$ $\displaystyle =$ $\displaystyle \left\{ {i:f_i (t^n ,\vec{c}^n ) < 0, i \in \{ 1,...,I\} } \right\}$  
$\displaystyle K^n$ $\displaystyle =$ $\displaystyle \left\{ {i:f_i (t^n ,\vec{c}^n ) + f_i (t^{n+1} ,\vec{c}^{(1)} ) < 0, i \in \{ 1,...,I\} } \right\}.$ (257)

The first step is identical to a step with the first-order EMP scheme. The second step mathmatically identical to a step with the first-order scheme if we rewrite it as


$\displaystyle {\vec{c}^{n + 1} = \vec{c}^n + \Delta t \: \vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)})\prod\limits_{k \in
K^n} {\frac{c_k^{n + 1} }{c_k^n }}}$
    with $\displaystyle \vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)}) = \frac{1}{2}\left...
...1},\vec{c}^{(1)})} \right)\prod\limits_{k \in K^n} {\frac{c_k^n }{c_k^{(1)} }}.$ (258)

Therefore, this scheme can be implemented as two consecutive steps with the first-order scheme, the second using $ \vec{h}(t^n,t^{n + 1},\vec{c}^n,\vec{c}^{(1)})$. The non-linear problem of each consecutive step is solved in auxiliary subroutine findp_bisection.

For more details, see Bruggeman et al. (2006).


USES:

    IMPLICIT NONE
INPUT PARAMETERS:
    REALTYPE, intent(in)                :: dt
    integer, intent(in)                 :: numc,nlev
   !INPUT/OUTPUT PARAMETER:
    REALTYPE, intent(inout)              :: cc(1:numc,0:nlev)
 
    interface
       subroutine right_hand_side(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): Jorn Bruggeman
LOCAL VARIABLES:
   logical  :: first
   integer  :: i,ci
   REALTYPE :: pi, rhs(1:numc,0:nlev), cc_med(1:numc,0:nlev), rhs_med(1:numc,0:nlev)

Karsten Bolding 2012-01-24