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)} }}}$ (258)

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\}.$ (259)

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)} }}.$ (260)

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-12-28