press_release/code/update_cg2d.F

C $Header: /u/gcmpack/MITgcm_contrib/ksnow/press_release/code/update_cg2d.F,v 1.3 2017/02/04 18:55:11 dgoldberg Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_PRESSURE_RELEASE_CODE
# include "SHELFICE_OPTIONS.h"
#endif

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 */
#ifdef ALLOW_PRESSURE_RELEASE_CODE
# include "SHELFICE.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:
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
#ifdef ALLOW_PRESSURE_RELEASE_CODE
      _RL     curPrelVisc, avgGrd
      _RL depth, eff_depth, depth_fac, pinit
      _RL drag_fac
#endif
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)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
            IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN        

             
             drag_fac = _recip_hFacW(i,j,k,bi,bj)**2*recip_drF(k)**2*
     &        pReleaseDamp * viscArNr(k) * 
     &        1./(1+exp(-10./cg2dminColumnEps*
     &         (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))

             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
     &               + faceArea*dU_psFacX(i,j,k,bi,bj)
     &               / (1 + drag_fac*deltaTmom)
     &                         *recip_dxC(i,j,bi,bj)
            ELSE
#endif
             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
     &               + faceArea*dU_psFacX(i,j,k,bi,bj)
     &                         *recip_dxC(i,j,bi,bj)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
            ENDIF
#endif

            faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
     &                *_hFacS(i,j,k,bi,bj)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
            IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        

             drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
     &        pReleaseDamp * viscArNr(k)* 
     &        1./(1+exp(-10./cg2dminColumnEps*
     &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))

            
             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
     &               + faceArea*dV_psFacY(i,j,k,bi,bj)
     &               / (1 + drag_fac*deltaTmom)
     &                         *recip_dyC(i,j,bi,bj)
            ELSE
#endif
             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
     &               + faceArea*dV_psFacY(i,j,k,bi,bj)
     &                         *recip_dyC(i,j,bi,bj)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
            ENDIF
#endif
           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

#ifdef ALLOW_PRESSURE_RELEASE_CODE
            IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
             drag_fac = _recip_hFacW(i,j,k,bi,bj)**2*recip_drF(k)**2*
     &        pReleaseDamp * viscArNr(k)* 
     &        1./(1+exp(-10./cg2dminColumnEps*
     &         (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2))) 
            ELSE
             drag_fac = 0. _d 0
            ENDIF
#endif
            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)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
     &       / (1 + drag_fac*deltaTmom)        
#endif
#ifdef ALLOW_PRESSURE_RELEASE_CODE
            IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN        
             drag_fac = _recip_hFacS(i,j,k,bi,bj)**2*recip_drF(k)**2*
     &        pReleaseDamp * viscArNr(k) * 
     &        1./(1+exp(-10./cg2dminColumnEps*
     &         (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
            ELSE
             drag_fac = 0. _d 0
            ENDIF
#endif
            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)
#ifdef ALLOW_PRESSURE_RELEASE_CODE
     &       / (1 + drag_fac*deltaTmom)        
#endif
           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
#ifdef ALLOW_PRESSURE_RELEASE_CODE
        DO j=1,sNy+1
         DO i=1,sNx+1
          IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN 
           IF (maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj).eq.1.) THEN

            avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i-1,j,bi,bj))            
            curPrelVisc = max((0.55+avgGrd*0.45),0. _d 0)*pReleaseVisc
            
            depth = depthColW(i,j,bi,bj)
            eff_depth = cg2dminColumnEps * 
     &       (1.-2./PI * COS (PI*depth / (2.*cg2dminColumnEps)))
            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj) + 
     &                        curPrelVisc*(eff_depth-depth)
     &                        *_dyG(i,j,bi,bj)*
     &                        recip_dxC(i,j,bi,bj) * cg2dNorm * 
     &                        implicSurfPress*implicDiv2DFlow /
     &                        deltaTmom
     
           ENDIF
          ENDIF
          IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN 
           IF (maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj).eq.1.) THEN

            avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i,j-1,bi,bj))            
            curPrelVisc = max((0.55+avgGrd*0.45),0. _d 0)*pReleaseVisc
            
            depth = depthColS(i,j,bi,bj)
            eff_depth = cg2dminColumnEps * 
     &       (1.-2./PI * COS (PI*depth / (2.*cg2dminColumnEps)))
            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj) + 
     &                        curPrelVisc*(eff_depth-depth)
     &                        *_dxG(i,j,bi,bj)*
     &                        recip_dyC(i,j,bi,bj) * cg2dNorm * 
     &                        implicSurfPress*implicDiv2DFlow /
     &                        deltaTmom
           ENDIF
          ENDIF
         ENDDO
        ENDDO
#endif
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	dgoldberg	1.4	C $Header: /u/gcmpack/MITgcm_contrib/ksnow/press_release/code/update_cg2d.F,v 1.3 2017/02/04 18:55:11 dgoldberg Exp $
2	ksnow	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6			#ifdef ALLOW_PRESSURE_RELEASE_CODE
7			# include "SHELFICE_OPTIONS.h"
8			#endif
9
10			CBOP
11			C !ROUTINE: UPDATE_CG2D
12			C !INTERFACE:
13			SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid )
14			C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| SUBROUTINE UPDATE_CG2D
17			C \| o Update 2d conjugate gradient solver operators
18			C \| account for Free-Surf effect on total column thickness
19			C ==========================================================
20			C \| This routine is based on INI_CG2D, and simplified. It is
21			C \| used when the non-linear free surface mode is activated
22			C \| or when bottom depth is part of the control vector.
23			C ==========================================================
24			C \ev
25
26			C !USES:
27			IMPLICIT NONE
28			C === Global variables ===
29			#include "SIZE.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32			#include "GRID.h"
33			#include "SURFACE.h"
34			#include "CG2D.h"
35	ksnow	1.2	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
36			# include "DYNVARS.h"
37			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
38	ksnow	1.1	#ifdef ALLOW_PRESSURE_RELEASE_CODE
39			# include "SHELFICE.h"
40			#endif
41
42			C !INPUT/OUTPUT PARAMETERS:
43			C === Routine arguments ===
44			C myTime :: Current time in simulation
45			C myIter :: Current iteration number in simulation
46			C myThid :: Thread number for this instance of the routine.
47			_RL myTime
48			INTEGER myIter
49			INTEGER myThid
50
51			C !LOCAL VARIABLES:
52			C-- Note : compared to "INI_CG2D", no needs to compute again
53			C the solver normalisation factor or the solver tolerance
54			C === Local variables ===
55			C bi,bj :: tile indices
56			C i,j,k :: Loop counters
57			C faceArea :: Temporary used to hold cell face areas.
58			INTEGER bi, bj
59			INTEGER i, j, k, ks
60			_RL faceArea
61			_RL pW_tmp, pS_tmp
62	ksnow	1.2	LOGICAL updatePreCond
63	ksnow	1.1	#ifdef ALLOW_PRESSURE_RELEASE_CODE
64			_RL curPrelVisc, avgGrd
65	ksnow	1.2	_RL depth, eff_depth, depth_fac, pinit
66			_RL drag_fac
67	ksnow	1.1	#endif
68			CEOP
69
70			C-- Decide when to update cg2d Preconditioner :
71			IF ( cg2dPreCondFreq.EQ.0 ) THEN
72			updatePreCond = .FALSE.
73			ELSE
74			updatePreCond = ( myIter.EQ.nIter0 )
75			IF ( MOD(myIter,cg2dPreCondFreq).EQ.0 ) updatePreCond=.TRUE.
76			ENDIF
77
78			C-- Initialise laplace operator
79			C aW2d: integral in Z Ax/dX
80			C aS2d: integral in Z Ay/dY
81			DO bj=myByLo(myThid),myByHi(myThid)
82			DO bi=myBxLo(myThid),myBxHi(myThid)
83			DO j=1-OLy,sNy+OLy
84			DO i=1-OLx,sNx+OLx
85			aW2d(i,j,bi,bj) = 0. _d 0
86			aS2d(i,j,bi,bj) = 0. _d 0
87			#ifdef ALLOW_AUTODIFF
88			aC2d(i,j,bi,bj) = 0. _d 0
89			#endif /* ALLOW_AUTODIFF */
90			ENDDO
91			ENDDO
92	ksnow	1.2	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
93			IF ( selectImplicitDrag.EQ.2 ) THEN
94			DO k=1,Nr
95			DO j=1,sNy+1
96			DO i=1,sNx+1
97			faceArea = _dyG(i,j,bi,bj)drF(k)deepFacC(k)*rhoFacC(k)
98			& *_hFacW(i,j,k,bi,bj)
99			#ifdef ALLOW_PRESSURE_RELEASE_CODE
100			IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN
101
102
103			drag_fac = _recip_hFacW(i,j,k,bi,bj)*2recip_drF(k)*2
104	dgoldberg	1.4	& pReleaseDamp * viscArNr(k) *
105	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
106			& (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
107
108			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
109			& + faceArea*dU_psFacX(i,j,k,bi,bj)
110			& / (1 + drag_fac*deltaTmom)
111	dgoldberg	1.4	& *recip_dxC(i,j,bi,bj)
112	ksnow	1.2	ELSE
113			#endif
114			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
115			& + faceArea*dU_psFacX(i,j,k,bi,bj)
116			& *recip_dxC(i,j,bi,bj)
117			#ifdef ALLOW_PRESSURE_RELEASE_CODE
118			ENDIF
119			#endif
120
121			faceArea = _dxG(i,j,bi,bj)drF(k)deepFacC(k)*rhoFacC(k)
122			& *_hFacS(i,j,k,bi,bj)
123			#ifdef ALLOW_PRESSURE_RELEASE_CODE
124			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
125
126			drag_fac = _recip_hFacS(i,j,k,bi,bj)*2recip_drF(k)*2
127	dgoldberg	1.4	& pReleaseDamp * viscArNr(k)*
128	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
129			& (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
130
131
132			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
133			& + faceArea*dV_psFacY(i,j,k,bi,bj)
134			& / (1 + drag_fac*deltaTmom)
135	dgoldberg	1.4	& *recip_dyC(i,j,bi,bj)
136	ksnow	1.2	ELSE
137			#endif
138			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
139			& + faceArea*dV_psFacY(i,j,k,bi,bj)
140			& *recip_dyC(i,j,bi,bj)
141			#ifdef ALLOW_PRESSURE_RELEASE_CODE
142			ENDIF
143			#endif
144			ENDDO
145			ENDDO
146			ENDDO
147			ELSE
148			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
149			DO k=1,Nr
150			DO j=1,sNy+1
151			DO i=1,sNx+1
152	ksnow	1.1	C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
153	ksnow	1.2
154			#ifdef ALLOW_PRESSURE_RELEASE_CODE
155			IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN
156			drag_fac = _recip_hFacW(i,j,k,bi,bj)*2recip_drF(k)*2
157	dgoldberg	1.4	& pReleaseDamp * viscArNr(k)*
158	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
159			& (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
160			ELSE
161			drag_fac = 0. _d 0
162			ENDIF
163	ksnow	1.1	#endif
164	ksnow	1.2	faceArea = _dyG(i,j,bi,bj)*drF(k)
165			& *_hFacW(i,j,k,bi,bj)
166
167
168			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
169			& + faceArea*recip_dxC(i,j,bi,bj)
170			#ifdef ALLOW_PRESSURE_RELEASE_CODE
171			& / (1 + drag_fac*deltaTmom)
172			#endif
173			#ifdef ALLOW_PRESSURE_RELEASE_CODE
174			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
175			drag_fac = _recip_hFacS(i,j,k,bi,bj)*2recip_drF(k)*2
176	dgoldberg	1.4	& pReleaseDamp * viscArNr(k) *
177	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
178			& (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
179			ELSE
180			drag_fac = 0. _d 0
181			ENDIF
182			#endif
183			faceArea = _dxG(i,j,bi,bj)*drF(k)
184			& *_hFacS(i,j,k,bi,bj)
185
186			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
187			& + faceArea*recip_dyC(i,j,bi,bj)
188			#ifdef ALLOW_PRESSURE_RELEASE_CODE
189			& / (1 + drag_fac*deltaTmom)
190			#endif
191			ENDDO
192	ksnow	1.1	ENDDO
193			ENDDO
194	ksnow	1.2	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
195			ENDIF
196			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
197	ksnow	1.1	DO j=1,sNy+1
198			DO i=1,sNx+1
199			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
200			& implicSurfPressimplicDiv2DFlow
201			#ifdef ALLOW_OBCS
202			& maskInC(i,j,bi,bj)maskInC(i-1,j,bi,bj)
203			#endif
204			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
205			& implicSurfPressimplicDiv2DFlow
206			#ifdef ALLOW_OBCS
207			& maskInC(i,j,bi,bj)maskInC(i,j-1,bi,bj)
208			#endif
209			ENDDO
210			ENDDO
211			#ifdef ALLOW_PRESSURE_RELEASE_CODE
212			DO j=1,sNy+1
213			DO i=1,sNx+1
214			IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN
215			IF (maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj).eq.1.) THEN
216
217			avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i-1,j,bi,bj))
218			curPrelVisc = max((0.55+avgGrd0.45),0. _d 0)pReleaseVisc
219
220			depth = depthColW(i,j,bi,bj)
221			eff_depth = cg2dminColumnEps *
222			& (1.-2./PI * COS (PIdepth / (2.cg2dminColumnEps)))
223			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj) +
224			& curPrelVisc*(eff_depth-depth)
225			& _dyG(i,j,bi,bj)
226			& recip_dxC(i,j,bi,bj) * cg2dNorm *
227			& implicSurfPress*implicDiv2DFlow /
228			& deltaTmom
229
230			ENDIF
231			ENDIF
232			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
233			IF (maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj).eq.1.) THEN
234
235			avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i,j-1,bi,bj))
236			curPrelVisc = max((0.55+avgGrd0.45),0. _d 0)pReleaseVisc
237
238			depth = depthColS(i,j,bi,bj)
239			eff_depth = cg2dminColumnEps *
240			& (1.-2./PI * COS (PIdepth / (2.cg2dminColumnEps)))
241			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj) +
242			& curPrelVisc*(eff_depth-depth)
243			& _dxG(i,j,bi,bj)
244			& recip_dyC(i,j,bi,bj) * cg2dNorm *
245			& implicSurfPress*implicDiv2DFlow /
246			& deltaTmom
247			ENDIF
248			ENDIF
249			ENDDO
250			ENDDO
251			#endif
252			C-- compute matrix main diagonal :
253			IF ( deepAtmosphere ) THEN
254			DO j=1,sNy
255			DO i=1,sNx
256			ks = kSurfC(i,j,bi,bj)
257			aC2d(i,j,bi,bj) = -(
258			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
259			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
260			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)*deepFac2F(ks)
261	ksnow	1.2	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
262	ksnow	1.1	& )
263			ENDDO
264			ENDDO
265			ELSE
266			DO j=1,sNy
267			DO i=1,sNx
268			aC2d(i,j,bi,bj) = -(
269			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
270			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
271			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)
272	ksnow	1.2	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
273	ksnow	1.1	& )
274			ENDDO
275			ENDDO
276			ENDIF
277			C- end bi,bj loops
278			ENDDO
279			ENDDO
280
281			IF ( updatePreCond ) THEN
282			C-- Update overlap regions
283			CALL EXCH_XY_RS(aC2d, myThid)
284
285			C-- Initialise preconditioner
286			DO bj=myByLo(myThid),myByHi(myThid)
287			DO bi=myBxLo(myThid),myBxHi(myThid)
288			DO j=1,sNy+1
289			DO i=1,sNx+1
290			IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
291			pC(i,j,bi,bj) = 1. _d 0
292			ELSE
293			pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj)
294			ENDIF
295			pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
296			IF ( pW_tmp .EQ. 0. ) THEN
297			pW(i,j,bi,bj) = 0.
298			ELSE
299			pW(i,j,bi,bj) =
300			& -aW2d(i,j,bi,bj)/((cg2dpcOffDFac pW_tmp)*2 )
301			ENDIF
302			pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
303			IF ( pS_tmp .EQ. 0. ) THEN
304			pS(i,j,bi,bj) = 0.
305			ELSE
306			pS(i,j,bi,bj) =
307			& -aS2d(i,j,bi,bj)/((cg2dpcOffDFac pS_tmp)*2 )
308			ENDIF
309			ENDDO
310			ENDDO
311			ENDDO
312			ENDDO
313			C- if update Preconditioner : end
314			ENDIF
315
316			RETURN
317			END