pkg/generic_advdiff/gad_u3c4_impl_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_u3c4_impl_r.F,v 1.3 2004/12/04 00:22:25 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C     !ROUTINE: GAD_FLUXLIMIT_IMPL_R
C     !INTERFACE:
      SUBROUTINE GAD_U3C4_IMPL_R( 
     I           bi,bj,k, iMin,iMax,jMin,jMax, 
     I           advectionScheme, deltaTarg, rTrans,
     O           a5d, b5d, c5d, d5d, e5d,
     I           myThid )

C     !DESCRIPTION:

C     Compute matrix element to solve vertical advection
C     \begin{enumerate}
C     \item implicitly using 3rd order upwind, or 
C     \item 4th order Centered advection schemes.
C     \end{enumerate}
C     Also, the contribution of vertical transport at interface k
C     is added to matrix lines k and k-1 

C     !USES:
      IMPLICIT NONE

C     == Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GAD.h"

C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj           :: tile indices
C k               :: vertical level
C iMin,iMax       :: computation domain
C jMin,jMax       :: computation domain
C advectionScheme :: advection scheme to use
C deltaTarg       :: time step
C rTrans          :: vertical volume transport
C tFld            :: tracer field
C a5d             :: 2nd  lower diag of pentadiagonal matrix 
C b5d             :: 1rst lower diag of pentadiagonal matrix 
C c5d             :: main diag       of pentadiagonal matrix 
C d5d             :: 1rst upper diag of pentadiagonal matrix 
C e5d             :: 2nd  upper diag of pentadiagonal matrix 
C myThid          :: thread number
      INTEGER bi,bj,k
      INTEGER iMin,iMax,jMin,jMax
      INTEGER advectionScheme
      _RL deltaTarg(Nr)
      _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL a5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL b5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL c5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL d5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL e5d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER myThid

C == Local Variables ==
C i,j             :: loop indices
C kp1             :: =min( k+1 , Nr )
C km2             :: =max( k-2 , 1 )
C rCenter         :: centered contribution
C rUpwind         :: upwind   contribution
      LOGICAL flagC4
      INTEGER i,j,kp1,km2
      _RL rCenter, rUpwind
      _RL rC4km, rC4kp, rU1k, rU3km, rU3kp
      _RL mskM, mskP, maskM2, maskP1
CEOP

C--   process interior interface only:
      IF ( k.GT.1 .AND. k.LE.Nr ) THEN

      km2=MAX(1,k-2)
      kp1=MIN(Nr,k+1)
      maskP1 = 1. _d 0 
      maskM2 = 1. _d 0 
      IF ( k.LE.2 ) maskM2 = 0. _d 0
      IF ( k.GE.Nr) maskP1 = 0. _d 0
      flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH 
     &         .AND. k.GT.2 .AND. k.LT.Nr

C--    Add centered & upwind contributions
       DO j=jMin,jMax
         DO i=iMin,iMax
           rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkFac
           mskM   = maskC(i,j,km2,bi,bj)*maskM2
           mskP   = maskC(i,j,kp1,bi,bj)*maskP1
           rC4km  = oneSixth*rCenter*mskM
           rC4kp  = oneSixth*rCenter*mskP
           IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
            rUpwind= 0. _d 0
            rU3km  = 0. _d 0
            rU3kp  = 0. _d 0
           ELSE
            rU1k   = oneSixth*abs(rCenter)
            rUpwind= rU1k+rU1k
            rU3km  = rU1k*mskM
            rU3kp  = rU1k*mskP
           ENDIF
           a5d(i,j,k)   = a5d(i,j,k)
     &                  - (rC4km - rU3km)
     &                   *deltaTarg(k)
     &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
           b5d(i,j,k)   = b5d(i,j,k)
     &                  + (rCenter + rC4km - rUpwind - rU3km)
     &                   *deltaTarg(k)
     &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
           c5d(i,j,k)   = c5d(i,j,k)
     &                  + (rCenter + rC4kp + rUpwind + rU3kp)
     &                   *deltaTarg(k)
     &                    *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
           d5d(i,j,k)   = d5d(i,j,k)
     &                  - (rC4kp + rU3kp)
     &                   *deltaTarg(k)
     &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
           b5d(i,j,k-1) = b5d(i,j,k-1)
     &                  + (rC4km - rU3km)
     &                   *deltaTarg(k-1)
     &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
           c5d(i,j,k-1) = c5d(i,j,k-1)
     &                  - (rCenter + rC4km - rUpwind - rU3km)
     &                   *deltaTarg(k-1)
     &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
           d5d(i,j,k-1) = d5d(i,j,k-1)
     &                  - (rCenter + rC4kp + rUpwind + rU3kp)
     &                   *deltaTarg(k-1)
     &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
           e5d(i,j,k-1) = e5d(i,j,k-1)
     &                  + (rC4kp + rU3kp)
     &                   *deltaTarg(k-1)
     &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
         ENDDO
       ENDDO

C--   process interior interface only: end
      ENDIF

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_u3c4_impl_r.F,v 1.3 2004/12/04 00:22:25 jmc Exp $
2	C $Name: $
3
4	#include "GAD_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: GAD_FLUXLIMIT_IMPL_R
8	C !INTERFACE:
9	SUBROUTINE GAD_U3C4_IMPL_R(
10	I bi,bj,k, iMin,iMax,jMin,jMax,
11	I advectionScheme, deltaTarg, rTrans,
12	O a5d, b5d, c5d, d5d, e5d,
13	I myThid )
14
15	C !DESCRIPTION:
16
17	C Compute matrix element to solve vertical advection
18	C \begin{enumerate}
19	C \item implicitly using 3rd order upwind, or
20	C \item 4th order Centered advection schemes.
21	C \end{enumerate}
22	C Also, the contribution of vertical transport at interface k
23	C is added to matrix lines k and k-1
24
25	C !USES:
26	IMPLICIT NONE
27
28	C == Global variables ===
29	#include "SIZE.h"
30	#include "GRID.h"
31	#include "EEPARAMS.h"
32	#include "PARAMS.h"
33	#include "GAD.h"
34
35	C !INPUT/OUTPUT PARAMETERS:
36	C == Routine Arguments ==
37	C bi,bj :: tile indices
38	C k :: vertical level
39	C iMin,iMax :: computation domain
40	C jMin,jMax :: computation domain
41	C advectionScheme :: advection scheme to use
42	C deltaTarg :: time step
43	C rTrans :: vertical volume transport
44	C tFld :: tracer field
45	C a5d :: 2nd lower diag of pentadiagonal matrix
46	C b5d :: 1rst lower diag of pentadiagonal matrix
47	C c5d :: main diag of pentadiagonal matrix
48	C d5d :: 1rst upper diag of pentadiagonal matrix
49	C e5d :: 2nd upper diag of pentadiagonal matrix
50	C myThid :: thread number
51	INTEGER bi,bj,k
52	INTEGER iMin,iMax,jMin,jMax
53	INTEGER advectionScheme
54	_RL deltaTarg(Nr)
55	_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56	_RL a5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57	_RL b5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58	_RL c5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59	_RL d5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
60	_RL e5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
61	INTEGER myThid
62
63	C == Local Variables ==
64	C i,j :: loop indices
65	C kp1 :: =min( k+1 , Nr )
66	C km2 :: =max( k-2 , 1 )
67	C rCenter :: centered contribution
68	C rUpwind :: upwind contribution
69	LOGICAL flagC4
70	INTEGER i,j,kp1,km2
71	_RL rCenter, rUpwind
72	_RL rC4km, rC4kp, rU1k, rU3km, rU3kp
73	_RL mskM, mskP, maskM2, maskP1
74	CEOP
75
76	C-- process interior interface only:
77	IF ( k.GT.1 .AND. k.LE.Nr ) THEN
78
79	km2=MAX(1,k-2)
80	kp1=MIN(Nr,k+1)
81	maskP1 = 1. _d 0
82	maskM2 = 1. _d 0
83	IF ( k.LE.2 ) maskM2 = 0. _d 0
84	IF ( k.GE.Nr) maskP1 = 0. _d 0
85	flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH
86	& .AND. k.GT.2 .AND. k.LT.Nr
87
88	C-- Add centered & upwind contributions
89	DO j=jMin,jMax
90	DO i=iMin,iMax
91	rCenter= 0.5 _d 0 rTrans(i,j)recip_rA(i,j,bi,bj)*rkFac
92	mskM = maskC(i,j,km2,bi,bj)*maskM2
93	mskP = maskC(i,j,kp1,bi,bj)*maskP1
94	rC4km = oneSixthrCentermskM
95	rC4kp = oneSixthrCentermskP
96	IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
97	rUpwind= 0. _d 0
98	rU3km = 0. _d 0
99	rU3kp = 0. _d 0
100	ELSE
101	rU1k = oneSixth*abs(rCenter)
102	rUpwind= rU1k+rU1k
103	rU3km = rU1k*mskM
104	rU3kp = rU1k*mskP
105	ENDIF
106	a5d(i,j,k) = a5d(i,j,k)
107	& - (rC4km - rU3km)
108	& *deltaTarg(k)
109	& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
110	b5d(i,j,k) = b5d(i,j,k)
111	& + (rCenter + rC4km - rUpwind - rU3km)
112	& *deltaTarg(k)
113	& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
114	c5d(i,j,k) = c5d(i,j,k)
115	& + (rCenter + rC4kp + rUpwind + rU3kp)
116	& *deltaTarg(k)
117	& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
118	d5d(i,j,k) = d5d(i,j,k)
119	& - (rC4kp + rU3kp)
120	& *deltaTarg(k)
121	& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
122	b5d(i,j,k-1) = b5d(i,j,k-1)
123	& + (rC4km - rU3km)
124	& *deltaTarg(k-1)
125	& recip_hFacC(i,j,k-1,bi,bj)recip_drF(k-1)
126	c5d(i,j,k-1) = c5d(i,j,k-1)
127	& - (rCenter + rC4km - rUpwind - rU3km)
128	& *deltaTarg(k-1)
129	& recip_hFacC(i,j,k-1,bi,bj)recip_drF(k-1)
130	d5d(i,j,k-1) = d5d(i,j,k-1)
131	& - (rCenter + rC4kp + rUpwind + rU3kp)
132	& *deltaTarg(k-1)
133	& recip_hFacC(i,j,k-1,bi,bj)recip_drF(k-1)
134	e5d(i,j,k-1) = e5d(i,j,k-1)
135	& + (rC4kp + rU3kp)
136	& *deltaTarg(k-1)
137	& recip_hFacC(i,j,k-1,bi,bj)recip_drF(k-1)
138	ENDDO
139	ENDDO
140
141	C-- process interior interface only: end
142	ENDIF
143
144	RETURN
145	END