pkg/generic_advdiff/gad_u3c4_impl_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_u3c4_impl_r.F,v 1.6 2005/10/20 17:03:08 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

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