model/src/update_cg2d.F

C $Header: /u/gcmpack/MITgcm/model/src/update_cg2d.F,v 1.10 2012/06/17 13:13:08 jmc 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     | used when the non-linear free surface mode is activated
C     | or when bottom depth is part of the control vector.
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_SOLVE4_PS_AND_DRAG
# include "DYNVARS.h"
#endif /* ALLOW_SOLVE4_PS_AND_DRAG */

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:
C-- Note : compared to "INI_CG2D", no needs to compute again
C   the solver normalisation factor or the solver tolerance
C     === Local variables ===
C     bi,bj  :: tile indices
C     i,j,k  :: Loop counters
C     faceArea :: Temporary used to hold cell face areas.
      INTEGER bi, bj
      INTEGER i, j, k, ks
      _RL     faceArea
      _RL     pW_tmp, pS_tmp
      LOGICAL updatePreCond
CEOP

C--   Decide when to update cg2d Preconditioner :
      IF ( cg2dPreCondFreq.EQ.0 ) THEN
        updatePreCond = .FALSE.
      ELSE
        updatePreCond = ( myIter.EQ.nIter0 )
        IF ( MOD(myIter,cg2dPreCondFreq).EQ.0 ) updatePreCond=.TRUE.
      ENDIF

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-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          aW2d(i,j,bi,bj) = 0. _d 0
          aS2d(i,j,bi,bj) = 0. _d 0
#ifdef ALLOW_AUTODIFF
          aC2d(i,j,bi,bj) = 0. _d 0
#endif /* ALLOW_AUTODIFF */
         ENDDO
        ENDDO
#ifdef ALLOW_SOLVE4_PS_AND_DRAG
        IF ( selectImplicitDrag.EQ.2 ) THEN
         DO k=1,Nr
          DO j=1,sNy+1
           DO i=1,sNx+1
            faceArea = _dyG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
     &                *_hFacW(i,j,k,bi,bj)
            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
     &               + faceArea*dU_psFacX(i,j,k,bi,bj)
     &                         *recip_dxC(i,j,bi,bj)
            faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
     &                *_hFacS(i,j,k,bi,bj)
            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
     &               + faceArea*dV_psFacY(i,j,k,bi,bj)
     &                         *recip_dyC(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ELSE
#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
         DO k=1,Nr
          DO j=1,sNy+1
           DO i=1,sNx+1
C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
            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_SOLVE4_PS_AND_DRAG
        ENDIF
#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
        DO j=1,sNy+1
         DO i=1,sNx+1
          aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
     &                     *implicSurfPress*implicDiv2DFlow
#ifdef ALLOW_OBCS
     &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
#endif
          aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
     &                     *implicSurfPress*implicDiv2DFlow
#ifdef ALLOW_OBCS
     &                  *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
#endif
         ENDDO
        ENDDO
C--   compute matrix main diagonal :
        IF ( deepAtmosphere ) THEN
         DO j=1,sNy
          DO i=1,sNx
           ks = kSurfC(i,j,bi,bj)
           aC2d(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)*deepFac2F(ks)
     &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
     &                        )
          ENDDO
         ENDDO
        ELSE
         DO j=1,sNy
          DO i=1,sNx
           aC2d(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
        ENDIF
C-    end bi,bj loops
       ENDDO
      ENDDO

      IF ( updatePreCond ) THEN
C--   Update overlap regions
      CALL EXCH_XY_RS(aC2d, 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 ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
            pC(i,j,bi,bj) = 1. _d 0
          ELSE
           pC(i,j,bi,bj) =  1. _d 0 / aC2d(i,j,bi,bj)
          ENDIF
          pW_tmp = aC2d(i,j,bi,bj)+aC2d(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 = aC2d(i,j,bi,bj)+aC2d(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
C-    if update Preconditioner : end
      ENDIF

      RETURN
      END
1	jmc	1.11	C $Header: /u/gcmpack/MITgcm/model/src/update_cg2d.F,v 1.10 2012/06/17 13:13:08 jmc 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	jmc	1.7	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	cnh	1.2	C ==========================================================
17	jmc	1.7	C \| This routine is based on INI_CG2D, and simplified. It is
18	jmc	1.9	C \| used when the non-linear free surface mode is activated
19			C \| or when bottom depth is part of the control vector.
20	cnh	1.2	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	jmc	1.11	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
33			# include "DYNVARS.h"
34			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
35	jmc	1.1
36	cnh	1.2	C !INPUT/OUTPUT PARAMETERS:
37	jmc	1.1	C === Routine arguments ===
38	jmc	1.8	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	jmc	1.1	_RL myTime
42			INTEGER myIter
43			INTEGER myThid
44	cnh	1.2
45			C !LOCAL VARIABLES:
46	jmc	1.1	C-- Note : compared to "INI_CG2D", no needs to compute again
47	jmc	1.9	C the solver normalisation factor or the solver tolerance
48	jmc	1.1	C === Local variables ===
49	jmc	1.7	C bi,bj :: tile indices
50	jmc	1.8	C i,j,k :: Loop counters
51	jmc	1.7	C faceArea :: Temporary used to hold cell face areas.
52	jmc	1.1	INTEGER bi, bj
53	jmc	1.8	INTEGER i, j, k, ks
54	jmc	1.1	_RL faceArea
55	jmc	1.3	_RL pW_tmp, pS_tmp
56	jmc	1.6	LOGICAL updatePreCond
57	cnh	1.2	CEOP
58	jmc	1.1
59	jmc	1.6	C-- Decide when to update cg2d Preconditioner :
60			IF ( cg2dPreCondFreq.EQ.0 ) THEN
61			updatePreCond = .FALSE.
62			ELSE
63			updatePreCond = ( myIter.EQ.nIter0 )
64			IF ( MOD(myIter,cg2dPreCondFreq).EQ.0 ) updatePreCond=.TRUE.
65			ENDIF
66
67	jmc	1.1	C-- Initialise laplace operator
68			C aW2d: integral in Z Ax/dX
69			C aS2d: integral in Z Ay/dY
70			DO bj=myByLo(myThid),myByHi(myThid)
71			DO bi=myBxLo(myThid),myBxHi(myThid)
72	jmc	1.10	DO j=1-OLy,sNy+OLy
73			DO i=1-OLx,sNx+OLx
74	jmc	1.8	aW2d(i,j,bi,bj) = 0. _d 0
75			aS2d(i,j,bi,bj) = 0. _d 0
76	jmc	1.10	#ifdef ALLOW_AUTODIFF
77			aC2d(i,j,bi,bj) = 0. _d 0
78			#endif /* ALLOW_AUTODIFF */
79	jmc	1.1	ENDDO
80			ENDDO
81	jmc	1.11	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
82			IF ( selectImplicitDrag.EQ.2 ) THEN
83			DO k=1,Nr
84			DO j=1,sNy+1
85			DO i=1,sNx+1
86			faceArea = _dyG(i,j,bi,bj)drF(k)deepFacC(k)*rhoFacC(k)
87			& *_hFacW(i,j,k,bi,bj)
88			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
89			& + faceArea*dU_psFacX(i,j,k,bi,bj)
90			& *recip_dxC(i,j,bi,bj)
91			faceArea = _dxG(i,j,bi,bj)drF(k)deepFacC(k)*rhoFacC(k)
92			& *_hFacS(i,j,k,bi,bj)
93			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
94			& + faceArea*dV_psFacY(i,j,k,bi,bj)
95			& *recip_dyC(i,j,bi,bj)
96			ENDDO
97			ENDDO
98			ENDDO
99			ELSE
100			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
101			DO k=1,Nr
102			DO j=1,sNy+1
103			DO i=1,sNx+1
104	jmc	1.7	C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
105	jmc	1.11	faceArea = _dyG(i,j,bi,bj)*drF(k)
106			& *_hFacW(i,j,k,bi,bj)
107			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
108			& + faceArea*recip_dxC(i,j,bi,bj)
109			faceArea = _dxG(i,j,bi,bj)*drF(k)
110			& *_hFacS(i,j,k,bi,bj)
111			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
112			& + faceArea*recip_dyC(i,j,bi,bj)
113			ENDDO
114	jmc	1.1	ENDDO
115			ENDDO
116	jmc	1.11	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
117			ENDIF
118			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
119	jmc	1.8	DO j=1,sNy+1
120			DO i=1,sNx+1
121			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
122			& implicSurfPressimplicDiv2DFlow
123	jmc	1.1	#ifdef ALLOW_OBCS
124	jmc	1.8	& maskInC(i,j,bi,bj)maskInC(i-1,j,bi,bj)
125	jmc	1.1	#endif
126	jmc	1.8	aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
127	jmc	1.1	& implicSurfPressimplicDiv2DFlow
128	jmc	1.8	#ifdef ALLOW_OBCS
129			& maskInC(i,j,bi,bj)maskInC(i,j-1,bi,bj)
130			#endif
131	jmc	1.1	ENDDO
132			ENDDO
133	jmc	1.7	C-- compute matrix main diagonal :
134			IF ( deepAtmosphere ) THEN
135	jmc	1.8	DO j=1,sNy
136			DO i=1,sNx
137	jmc	1.9	ks = kSurfC(i,j,bi,bj)
138	jmc	1.8	aC2d(i,j,bi,bj) = -(
139			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
140			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
141			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)*deepFac2F(ks)
142	jmc	1.11	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
143	jmc	1.7	& )
144			ENDDO
145	jmc	1.3	ENDDO
146	jmc	1.7	ELSE
147	jmc	1.8	DO j=1,sNy
148			DO i=1,sNx
149			aC2d(i,j,bi,bj) = -(
150			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
151			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
152			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)
153	jmc	1.11	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
154	jmc	1.7	& )
155			ENDDO
156	jmc	1.6	ENDDO
157			ENDIF
158			C- end bi,bj loops
159	jmc	1.1	ENDDO
160			ENDDO
161	jmc	1.7
162	jmc	1.6	IF ( updatePreCond ) THEN
163	jmc	1.1	C-- Update overlap regions
164	jmc	1.7	CALL EXCH_XY_RS(aC2d, myThid)
165	jmc	1.1
166			C-- Initialise preconditioner
167			DO bj=myByLo(myThid),myByHi(myThid)
168			DO bi=myBxLo(myThid),myBxHi(myThid)
169	jmc	1.8	DO j=1,sNy+1
170			DO i=1,sNx+1
171			IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
172			pC(i,j,bi,bj) = 1. _d 0
173	jmc	1.1	ELSE
174	jmc	1.8	pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj)
175	jmc	1.1	ENDIF
176	jmc	1.8	pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
177	jmc	1.3	IF ( pW_tmp .EQ. 0. ) THEN
178	jmc	1.8	pW(i,j,bi,bj) = 0.
179	jmc	1.1	ELSE
180	jmc	1.8	pW(i,j,bi,bj) =
181			& -aW2d(i,j,bi,bj)/((cg2dpcOffDFac pW_tmp)*2 )
182	jmc	1.1	ENDIF
183	jmc	1.8	pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
184	jmc	1.3	IF ( pS_tmp .EQ. 0. ) THEN
185	jmc	1.8	pS(i,j,bi,bj) = 0.
186	jmc	1.1	ELSE
187	jmc	1.8	pS(i,j,bi,bj) =
188			& -aS2d(i,j,bi,bj)/((cg2dpcOffDFac pS_tmp)*2 )
189	jmc	1.1	ENDIF
190			ENDDO
191			ENDDO
192			ENDDO
193			ENDDO
194	jmc	1.6	C- if update Preconditioner : end
195			ENDIF
196	jmc	1.1
197			RETURN
198			END