model/src/update_cg2d.F

C $Header: /u/gcmpack/MITgcm/model/src/update_cg2d.F,v 1.4 2003/10/09 04:19:18 edhill Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: UPDATE_CG2D
C     !INTERFACE:
      SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE UPDATE_CG2D                                    
C     | o Update 2d conjugate gradient solver operators           
C     |   account for Free-Surf effect on total column thickness  
C     *==========================================================*
C     | This routine is based on INI_CG2D, and simplified. It is 
C     | only active when the non-linear free surface mode of 
C     | equations is active.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "CG2D.h"
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime - Current time in simulation
C     myIter - Current iteration number in simulation
C     myThid - Thread number for this instance of the routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
#ifdef NONLIN_FRSURF
C-- Note : compared to "INI_CG2D", no needs to compute again
C   the solver norn=malisation factor of the solver tolerance
C     === Local variables ===
C     bi,bj  - Loop counters
C     I,J,K
C     faceArea - Temporary used to hold cell face areas.
      INTEGER bi, bj
      INTEGER I, J, K
      _RL     faceArea
      _RL     pW_tmp, pS_tmp
CEOP

C--   Initialise laplace operator
C     aW2d: integral in Z Ax/dX
C     aS2d: integral in Z Ay/dY
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1,sNy+1
         DO I=1,sNx+1
          aW2d(I,J,bi,bj) = 0. _d 0
          aS2d(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
        DO K=1,Nr
         DO J=1,sNy+1
          DO I=1,sNx+1
           faceArea = _dyG(I,J,bi,bj)*drF(K)
     &               *_hFacW(I,J,K,bi,bj)
           aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) + 
     &             faceArea*recip_dxC(I,J,bi,bj)
           faceArea = _dxG(I,J,bi,bj)*drF(K)
     &               *_hFacS(I,J,K,bi,bj)
           aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) + 
     &             faceArea*recip_dyC(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
         DO I=1,sNx
          IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
          IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
          IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
          IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
         ENDDO
         DO J=1,sNy
          IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
          IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
          IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
          IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
         ENDDO
        ENDIF
#endif
        DO J=1,sNy+1
         DO I=1,sNx+1
          aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm
     &                     *implicSurfPress*implicDiv2DFlow
          aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm
     &                     *implicSurfPress*implicDiv2DFlow
         ENDDO
        ENDDO
C--   Start to compute preconditioner matrix (use cg2d_q as temporary array)
        DO J=1,sNy
         DO I=1,sNx
          cg2d_q(I,J,bi,bj) = -(
     &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
     &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
     &    +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)* 
     &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
     &    )
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      
C--   Update overlap regions
      _EXCH_XY_R8(cg2d_q, myThid)

C--   Initialise preconditioner
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1,sNy+1
         DO I=1,sNx+1
          IF ( cg2d_q(I,J,bi,bj) .EQ. 0. ) THEN
            pC(I,J,bi,bj) = 1. _d 0
          ELSE
           pC(I,J,bi,bj) =  1. _d 0 / cg2d_q(I,J,bi,bj)
          ENDIF
          pW_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I-1,J,bi,bj)
          IF ( pW_tmp .EQ. 0. ) THEN
           pW(I,J,bi,bj) = 0.
          ELSE
           pW(I,J,bi,bj) = 
     &     -aW2d(I,J,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 )
          ENDIF
          pS_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I,J-1,bi,bj)
          IF ( pS_tmp .EQ. 0. ) THEN
           pS(I,J,bi,bj) = 0.
          ELSE
           pS(I,J,bi,bj) =
     &     -aS2d(I,J,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 )
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#endif /* NONLIN_FRSURF */

      RETURN
      END
1	jmc	1.5	C $Header: /u/gcmpack/MITgcm/model/src/update_cg2d.F,v 1.4 2003/10/09 04:19:18 edhill Exp $
2	jmc	1.1	C $Name: $
3
4	edhill	1.4	#include "PACKAGES_CONFIG.h"
5	jmc	1.1	#include "CPP_OPTIONS.h"
6	edhill	1.4
7	cnh	1.2	CBOP
8			C !ROUTINE: UPDATE_CG2D
9			C !INTERFACE:
10			SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid )
11			C !DESCRIPTION: \bv
12			C ==========================================================
13			C \| SUBROUTINE UPDATE_CG2D
14			C \| o Update 2d conjugate gradient solver operators
15			C \| account for Free-Surf effect on total column thickness
16			C ==========================================================
17			C \| This routine is based on INI_CG2D, and simplified. It is
18			C \| only active when the non-linear free surface mode of
19			C \| equations is active.
20			C ==========================================================
21			C \ev
22	jmc	1.1
23	cnh	1.2	C !USES:
24	jmc	1.1	IMPLICIT NONE
25			C === Global variables ===
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "SURFACE.h"
31			#include "CG2D.h"
32			#ifdef ALLOW_OBCS
33			#include "OBCS.h"
34			#endif
35
36	cnh	1.2	C !INPUT/OUTPUT PARAMETERS:
37	jmc	1.1	C === Routine arguments ===
38			C myTime - Current time in simulation
39			C myIter - Current iteration number in simulation
40			C myThid - Thread number for this instance of the routine.
41			_RL myTime
42			INTEGER myIter
43			INTEGER myThid
44	cnh	1.2
45			C !LOCAL VARIABLES:
46	jmc	1.1	#ifdef NONLIN_FRSURF
47			C-- Note : compared to "INI_CG2D", no needs to compute again
48			C the solver norn=malisation factor of the solver tolerance
49			C === Local variables ===
50			C bi,bj - Loop counters
51			C I,J,K
52			C faceArea - Temporary used to hold cell face areas.
53			INTEGER bi, bj
54			INTEGER I, J, K
55			_RL faceArea
56	jmc	1.3	_RL pW_tmp, pS_tmp
57	cnh	1.2	CEOP
58	jmc	1.1
59			C-- Initialise laplace operator
60			C aW2d: integral in Z Ax/dX
61			C aS2d: integral in Z Ay/dY
62			DO bj=myByLo(myThid),myByHi(myThid)
63			DO bi=myBxLo(myThid),myBxHi(myThid)
64	jmc	1.3	DO J=1,sNy+1
65			DO I=1,sNx+1
66	jmc	1.1	aW2d(I,J,bi,bj) = 0. _d 0
67			aS2d(I,J,bi,bj) = 0. _d 0
68			ENDDO
69			ENDDO
70			DO K=1,Nr
71	jmc	1.3	DO J=1,sNy+1
72			DO I=1,sNx+1
73	jmc	1.1	faceArea = _dyG(I,J,bi,bj)*drF(K)
74			& *_hFacW(I,J,K,bi,bj)
75			aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
76			& faceArea*recip_dxC(I,J,bi,bj)
77			faceArea = _dxG(I,J,bi,bj)*drF(K)
78			& *_hFacS(I,J,K,bi,bj)
79			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
80			& faceArea*recip_dyC(I,J,bi,bj)
81			ENDDO
82			ENDDO
83			ENDDO
84			#ifdef ALLOW_OBCS
85			IF (useOBCS) THEN
86			DO I=1,sNx
87			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
88			IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
89			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
90			IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
91			ENDDO
92			DO J=1,sNy
93			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
94			IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
95			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
96			IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
97			ENDDO
98			ENDIF
99			#endif
100	jmc	1.3	DO J=1,sNy+1
101			DO I=1,sNx+1
102	jmc	1.1	aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm
103			& implicSurfPressimplicDiv2DFlow
104			aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm
105			& implicSurfPressimplicDiv2DFlow
106			ENDDO
107			ENDDO
108	jmc	1.3	C-- Start to compute preconditioner matrix (use cg2d_q as temporary array)
109			DO J=1,sNy
110			DO I=1,sNx
111			cg2d_q(I,J,bi,bj) = -(
112			& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
113			& +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
114			& +freeSurfFaccg2dNormrecip_Bo(I,J,bi,bj)*
115	jmc	1.5	& rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
116	jmc	1.3	& )
117			ENDDO
118			ENDDO
119	jmc	1.1	ENDDO
120			ENDDO
121
122			C-- Update overlap regions
123	jmc	1.3	_EXCH_XY_R8(cg2d_q, myThid)
124	jmc	1.1
125			C-- Initialise preconditioner
126			DO bj=myByLo(myThid),myByHi(myThid)
127			DO bi=myBxLo(myThid),myBxHi(myThid)
128	jmc	1.3	DO J=1,sNy+1
129			DO I=1,sNx+1
130			IF ( cg2d_q(I,J,bi,bj) .EQ. 0. ) THEN
131	jmc	1.1	pC(I,J,bi,bj) = 1. _d 0
132			ELSE
133	jmc	1.3	pC(I,J,bi,bj) = 1. _d 0 / cg2d_q(I,J,bi,bj)
134	jmc	1.1	ENDIF
135	jmc	1.3	pW_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I-1,J,bi,bj)
136			IF ( pW_tmp .EQ. 0. ) THEN
137	jmc	1.1	pW(I,J,bi,bj) = 0.
138			ELSE
139			pW(I,J,bi,bj) =
140	jmc	1.3	& -aW2d(I,J,bi,bj)/((cg2dpcOffDFac pW_tmp)*2 )
141	jmc	1.1	ENDIF
142	jmc	1.3	pS_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I,J-1,bi,bj)
143			IF ( pS_tmp .EQ. 0. ) THEN
144	jmc	1.1	pS(I,J,bi,bj) = 0.
145			ELSE
146			pS(I,J,bi,bj) =
147	jmc	1.3	& -aS2d(I,J,bi,bj)/((cg2dpcOffDFac pS_tmp)*2 )
148	jmc	1.1	ENDIF
149			ENDDO
150			ENDDO
151			ENDDO
152			ENDDO
153
154			#endif /* NONLIN_FRSURF */
155
156			RETURN
157			END