model/src/adams_bashforth3.F

6b924d96b1 Jean*0001 #include "CPP_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: ADAMS_BASHFORTH3
                0005 C     !INTERFACE:
                0006       SUBROUTINE ADAMS_BASHFORTH3(
94e3f14b29 Jean*0007      I                     bi, bj, kArg, kSize,
6b924d96b1 Jean*0008      U                     gTracer, gTrNm,
38e8b1b15b Jean*0009      O                     AB_gTr,
43065bee0a Jean*0010      I                     startAB, myIter, myThid )
6b924d96b1 Jean*0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
43065bee0a Jean*0013 C     | S/R ADAMS_BASHFORTH3
                0014 C     | o Extrapolate forward in time using third order
                0015 C     |   Adams-Bashforth method.
                0016 C     *==========================================================*
                0017 C     | Either apply to tendency (kArg>0) at level k=kArg,
                0018 C     |     or apply to state variable (kArg=0) for all levels
6b924d96b1 Jean*0019 C     *==========================================================*
                0020 C     \ev
43065bee0a Jean*0021 C Extrapolate forward in time using 2 A.B. parameters (alpha,beta),
                0022 C either tendency gX :
6b924d96b1 Jean*0023 C \begin{equation*}
                0024 C gX^{n+1/2} = (1 + \alpha + \beta) gX^{n}
                0025 C              - (\alpha + 2 \beta) gX^{n-1}
                0026 C                          + \beta  gX^{n-2}
                0027 C \end{equation*}
43065bee0a Jean*0028 C     or state variable X :
                0029 C \begin{equation*}
                0030 C  X^{n+1/2} = (1 + \alpha + \beta) X^{n}
                0031 C              - (\alpha + 2 \beta) X^{n-1}
                0032 C                          + \beta  X^{n-2}
                0033 C \end{equation*}
6b924d96b1 Jean*0034 C with:
                0035 C (alpha,beta)=(1/2,5/12) : AB-3, stable until CFL = 0.724
                0036 C     (note: beta=0.281105 give the Max stability: 0.78616)
                0037 C (alpha,beta)=(1/2+abEps,0) : return to previous quasi AB-2.
38e8b1b15b Jean*0038 C (alpha,beta)=(0,0)         : 1rst.Order forward time stepping
6b924d96b1 Jean*0039
                0040 C     !USES:
                0041       IMPLICIT NONE
                0042 C     == Global variables ===
                0043 #include "SIZE.h"
                0044 #include "EEPARAMS.h"
                0045 #include "PARAMS.h"
                0046
                0047 C     !INPUT/OUTPUT PARAMETERS:
                0048 C     == Routine Arguments ==
43065bee0a Jean*0049 C     bi,bj   :: Tile indices
                0050 C     kArg    :: if >0: apply AB on tendency at level k=kArg
                0051 C             :: if =0: apply AB on state variable and process all levels
fed646132e Jean*0052 C     kSize   :: 3rd dimension of gTracer
38e8b1b15b Jean*0053 C     gTracer ::  in: Tendency/State at current time
23a7f3050f Jean*0054 C             :: out(kArg >0): Extrapolated Tendency at current time
                0055 C     gTrNm   ::  in: Tendency/State at previous time
                0056 C             :: out(kArg >0): Save tendency at current time
                0057 C             :: out(kArg =0): Extrapolated State at current time
38e8b1b15b Jean*0058 C     AB_gTr  :: Adams-Bashforth tendency increment
117ee419f5 Jean*0059 C     startAB :: number of previous time level available to start/restart AB
6b924d96b1 Jean*0060 C     myIter  :: Current time step number
38e8b1b15b Jean*0061 C     myThid  :: my Thread Id. number
94e3f14b29 Jean*0062       INTEGER bi, bj, kArg, kSize
23a7f3050f Jean*0063       _RL  gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize)
fed646132e Jean*0064       _RL  gTrNm  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,2)
38e8b1b15b Jean*0065       _RL  AB_gTr (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117ee419f5 Jean*0066       INTEGER startAB
6b924d96b1 Jean*0067       INTEGER myIter, myThid
                0068
                0069 #ifdef ALLOW_ADAMSBASHFORTH_3
                0070 C     !LOCAL VARIABLES:
                0071 C     == Local variables ==
fed646132e Jean*0072 C     i, j        :: Loop counters
                0073 C     k, kl       :: level indices
6b924d96b1 Jean*0074 C     m1,m2       :: indices for the 2 previous time-step Tendency
                0075 C     ab1,ab2,ab3 :: Adams bashforth extrapolation weights.
fed646132e Jean*0076       INTEGER i,j, k, kl, m1,m2
6b924d96b1 Jean*0077       _RL ab0, ab1, ab2
                0078 CEOP
                0079
972c0130ec Jean*0080       m1 = 1 + MOD(myIter+1,2)
                0081       m2 = 1 + MOD( myIter ,2)
6b924d96b1 Jean*0082
                0083 C     Adams-Bashforth timestepping weights
117ee419f5 Jean*0084       IF ( myIter.EQ.nIter0 .AND. startAB.EQ.0 ) THEN
38e8b1b15b Jean*0085        ab0 = 0. _d 0
6b924d96b1 Jean*0086        ab1 = 0. _d 0
                0087        ab2 = 0. _d 0
117ee419f5 Jean*0088       ELSEIF ( (myIter.EQ.nIter0   .AND. startAB.EQ.1)
                0089      &    .OR. (myIter.EQ.1+nIter0 .AND. startAB.EQ.0) ) THEN
38e8b1b15b Jean*0090        ab0 =  alph_AB
6b924d96b1 Jean*0091        ab1 = -alph_AB
                0092        ab2 = 0. _d 0
                0093       ELSE
38e8b1b15b Jean*0094        ab0 =  alph_AB + beta_AB
6b924d96b1 Jean*0095        ab1 = -alph_AB - 2.*beta_AB
                0096        ab2 =  beta_AB
                0097       ENDIF
                0098
                0099 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0100
43065bee0a Jean*0101       IF ( kArg.EQ.0 ) THEN
                0102 C-    Extrapolate forward in time the state variable, with AB weights:
94e3f14b29 Jean*0103         DO k=1,kSize
                0104          DO j=1-OLy,sNy+OLy
                0105           DO i=1-OLx,sNx+OLx
23a7f3050f Jean*0106            AB_gTr(i,j) = ab0*gTracer(i,j,k)
38e8b1b15b Jean*0107      &                 + ab1*gTrNm(i,j,k,bi,bj,m1)
                0108      &                 + ab2*gTrNm(i,j,k,bi,bj,m2)
23a7f3050f Jean*0109            gTrNm(i,j,k,bi,bj,m2) = gTracer(i,j,k) + AB_gTr(i,j)
43065bee0a Jean*0110           ENDDO
                0111          ENDDO
                0112         ENDDO
                0113       ELSE
                0114 C-    Extrapolate forward in time the tendency, with AB weights:
fed646132e Jean*0115         kl = kArg
                0116         k = MIN( kArg, kSize )
94e3f14b29 Jean*0117         DO j=1-OLy,sNy+OLy
                0118          DO i=1-OLx,sNx+OLx
23a7f3050f Jean*0119            AB_gTr(i,j) = ab0*gTracer(i,j,k)
fed646132e Jean*0120      &                 + ab1*gTrNm(i,j,kl,bi,bj,m1)
                0121      &                 + ab2*gTrNm(i,j,kl,bi,bj,m2)
                0122            gTrNm(i,j,kl,bi,bj,m2) = gTracer(i,j,k)
23a7f3050f Jean*0123            gTracer(i,j,k) = gTracer(i,j,k) + AB_gTr(i,j)
43065bee0a Jean*0124          ENDDO
                0125         ENDDO
                0126 C---
                0127       ENDIF
6b924d96b1 Jean*0128
                0129 #endif /* ALLOW_ADAMSBASHFORTH_3 */
                0130
                0131       RETURN
                0132       END