press_release/code/update_cg2d.F

C $Header: /u/gcmpack/MITgcm_contrib/ksnow/press_release/code/update_cg2d.F,v 1.2 2017/01/30 16:32:17 ksnow 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*
     &        12000. * 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)
            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*
     &        12000. * 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)
            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*
     &        12000. * 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*
     &        12000. * 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.3	C $Header: /u/gcmpack/MITgcm_contrib/ksnow/press_release/code/update_cg2d.F,v 1.2 2017/01/30 16:32:17 ksnow 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.3	& 12000. * 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			ELSE
112			#endif
113			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
114			& + faceArea*dU_psFacX(i,j,k,bi,bj)
115			& *recip_dxC(i,j,bi,bj)
116			#ifdef ALLOW_PRESSURE_RELEASE_CODE
117			ENDIF
118			#endif
119
120			faceArea = _dxG(i,j,bi,bj)drF(k)deepFacC(k)*rhoFacC(k)
121			& *_hFacS(i,j,k,bi,bj)
122			#ifdef ALLOW_PRESSURE_RELEASE_CODE
123			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
124
125			drag_fac = _recip_hFacS(i,j,k,bi,bj)*2recip_drF(k)*2
126	dgoldberg	1.3	& 12000. * viscArNr(k)*
127	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
128			& (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
129
130
131			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
132			& + faceArea*dV_psFacY(i,j,k,bi,bj)
133			& / (1 + drag_fac*deltaTmom)
134			ELSE
135			#endif
136			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
137			& + faceArea*dV_psFacY(i,j,k,bi,bj)
138			& *recip_dyC(i,j,bi,bj)
139			#ifdef ALLOW_PRESSURE_RELEASE_CODE
140			ENDIF
141			#endif
142			ENDDO
143			ENDDO
144			ENDDO
145			ELSE
146			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
147			DO k=1,Nr
148			DO j=1,sNy+1
149			DO i=1,sNx+1
150	ksnow	1.1	C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
151	ksnow	1.2
152			#ifdef ALLOW_PRESSURE_RELEASE_CODE
153			IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN
154			drag_fac = _recip_hFacW(i,j,k,bi,bj)*2recip_drF(k)*2
155	dgoldberg	1.3	& 12000. * viscArNr(k)*
156	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
157			& (-1.)*(depthColW(i,j,bi,bj)-cg2dminColumnEps/2)))
158			ELSE
159			drag_fac = 0. _d 0
160			ENDIF
161	ksnow	1.1	#endif
162	ksnow	1.2	faceArea = _dyG(i,j,bi,bj)*drF(k)
163			& *_hFacW(i,j,k,bi,bj)
164
165
166			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
167			& + faceArea*recip_dxC(i,j,bi,bj)
168			#ifdef ALLOW_PRESSURE_RELEASE_CODE
169			& / (1 + drag_fac*deltaTmom)
170			#endif
171			#ifdef ALLOW_PRESSURE_RELEASE_CODE
172			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
173			drag_fac = _recip_hFacS(i,j,k,bi,bj)*2recip_drF(k)*2
174	dgoldberg	1.3	& 12000. * viscArNr(k) *
175	ksnow	1.2	& 1./(1+exp(-10./cg2dminColumnEps*
176			& (-1.)*(depthColS(i,j,bi,bj)-cg2dminColumnEps/2)))
177			ELSE
178			drag_fac = 0. _d 0
179			ENDIF
180			#endif
181			faceArea = _dxG(i,j,bi,bj)*drF(k)
182			& *_hFacS(i,j,k,bi,bj)
183
184			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
185			& + faceArea*recip_dyC(i,j,bi,bj)
186			#ifdef ALLOW_PRESSURE_RELEASE_CODE
187			& / (1 + drag_fac*deltaTmom)
188			#endif
189			ENDDO
190	ksnow	1.1	ENDDO
191			ENDDO
192	ksnow	1.2	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
193			ENDIF
194			#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
195	ksnow	1.1	DO j=1,sNy+1
196			DO i=1,sNx+1
197			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
198			& implicSurfPressimplicDiv2DFlow
199			#ifdef ALLOW_OBCS
200			& maskInC(i,j,bi,bj)maskInC(i-1,j,bi,bj)
201			#endif
202			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
203			& implicSurfPressimplicDiv2DFlow
204			#ifdef ALLOW_OBCS
205			& maskInC(i,j,bi,bj)maskInC(i,j-1,bi,bj)
206			#endif
207			ENDDO
208			ENDDO
209			#ifdef ALLOW_PRESSURE_RELEASE_CODE
210			DO j=1,sNy+1
211			DO i=1,sNx+1
212			IF (depthColW(i,j,bi,bj).lt.cg2dminColumnEps) THEN
213			IF (maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj).eq.1.) THEN
214
215			avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i-1,j,bi,bj))
216			curPrelVisc = max((0.55+avgGrd0.45),0. _d 0)pReleaseVisc
217
218			depth = depthColW(i,j,bi,bj)
219			eff_depth = cg2dminColumnEps *
220			& (1.-2./PI * COS (PIdepth / (2.cg2dminColumnEps)))
221			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj) +
222			& curPrelVisc*(eff_depth-depth)
223			& _dyG(i,j,bi,bj)
224			& recip_dxC(i,j,bi,bj) * cg2dNorm *
225			& implicSurfPress*implicDiv2DFlow /
226			& deltaTmom
227
228			ENDIF
229			ENDIF
230			IF (depthColS(i,j,bi,bj).lt.cg2dminColumnEps) THEN
231			IF (maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj).eq.1.) THEN
232
233			avgGrd = min (grdFactor(i,j,bi,bj),grdFactor(i,j-1,bi,bj))
234			curPrelVisc = max((0.55+avgGrd0.45),0. _d 0)pReleaseVisc
235
236			depth = depthColS(i,j,bi,bj)
237			eff_depth = cg2dminColumnEps *
238			& (1.-2./PI * COS (PIdepth / (2.cg2dminColumnEps)))
239			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj) +
240			& curPrelVisc*(eff_depth-depth)
241			& _dxG(i,j,bi,bj)
242			& recip_dyC(i,j,bi,bj) * cg2dNorm *
243			& implicSurfPress*implicDiv2DFlow /
244			& deltaTmom
245			ENDIF
246			ENDIF
247			ENDDO
248			ENDDO
249			#endif
250			C-- compute matrix main diagonal :
251			IF ( deepAtmosphere ) THEN
252			DO j=1,sNy
253			DO i=1,sNx
254			ks = kSurfC(i,j,bi,bj)
255			aC2d(i,j,bi,bj) = -(
256			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
257			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
258			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)*deepFac2F(ks)
259	ksnow	1.2	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
260	ksnow	1.1	& )
261			ENDDO
262			ENDDO
263			ELSE
264			DO j=1,sNy
265			DO i=1,sNx
266			aC2d(i,j,bi,bj) = -(
267			& aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
268			& +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
269			& +freeSurfFaccg2dNormrecip_Bo(i,j,bi,bj)
270	ksnow	1.2	& *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
271	ksnow	1.1	& )
272			ENDDO
273			ENDDO
274			ENDIF
275			C- end bi,bj loops
276			ENDDO
277			ENDDO
278
279			IF ( updatePreCond ) THEN
280			C-- Update overlap regions
281			CALL EXCH_XY_RS(aC2d, myThid)
282
283			C-- Initialise preconditioner
284			DO bj=myByLo(myThid),myByHi(myThid)
285			DO bi=myBxLo(myThid),myBxHi(myThid)
286			DO j=1,sNy+1
287			DO i=1,sNx+1
288			IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
289			pC(i,j,bi,bj) = 1. _d 0
290			ELSE
291			pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj)
292			ENDIF
293			pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
294			IF ( pW_tmp .EQ. 0. ) THEN
295			pW(i,j,bi,bj) = 0.
296			ELSE
297			pW(i,j,bi,bj) =
298			& -aW2d(i,j,bi,bj)/((cg2dpcOffDFac pW_tmp)*2 )
299			ENDIF
300			pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
301			IF ( pS_tmp .EQ. 0. ) THEN
302			pS(i,j,bi,bj) = 0.
303			ELSE
304			pS(i,j,bi,bj) =
305			& -aS2d(i,j,bi,bj)/((cg2dpcOffDFac pS_tmp)*2 )
306			ENDIF
307			ENDDO
308			ENDDO
309			ENDDO
310			ENDDO
311			C- if update Preconditioner : end
312			ENDIF
313
314			RETURN
315			END