39 |
|
|
40 |
C !LOCAL VARIABLES: |
C !LOCAL VARIABLES: |
41 |
C === Local variables === |
C === Local variables === |
42 |
C xG, yG - Global coordinate location. |
C bi,bj :: tile indices |
43 |
C zG |
C I,J,K :: Loop counters |
|
C iG, jG - Global coordinate index |
|
|
C bi,bj - Loop counters |
|
44 |
C faceArea - Temporary used to hold cell face areas. |
C faceArea - Temporary used to hold cell face areas. |
|
C I,J,K |
|
45 |
C myNorm - Work variable used in calculating normalisation factor |
C myNorm - Work variable used in calculating normalisation factor |
46 |
C sumArea - Work variable used to compute the total Domain Area |
C sumArea - Work variable used to compute the total Domain Area |
47 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
48 |
INTEGER bi, bj |
INTEGER bi, bj |
49 |
INTEGER I, J, K |
INTEGER I, J, K, ks |
50 |
_RL faceArea |
_RL faceArea |
51 |
_RS myNorm |
_RS myNorm |
52 |
_RL aC, aCw, aCs |
_RS aC, aCw, aCs |
53 |
CEOP |
CEOP |
54 |
|
|
55 |
C-- Initialize arrays in common blocs (CG2D.h) ; not really necessary |
C-- Initialize arrays in common blocs (CG2D.h) ; not really necessary |
56 |
C but safer when EXCH do not fill all the overlap regions. |
C but safer when EXCH do not fill all the overlap regions. |
57 |
DO bj=myByLo(myThid),myByHi(myThid) |
DO bj=myByLo(myThid),myByHi(myThid) |
58 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
DO bi=myBxLo(myThid),myBxHi(myThid) |
60 |
DO I=1-OLx,sNx+OLx |
DO I=1-OLx,sNx+OLx |
61 |
aW2d(I,J,bi,bj) = 0. _d 0 |
aW2d(I,J,bi,bj) = 0. _d 0 |
62 |
aS2d(I,J,bi,bj) = 0. _d 0 |
aS2d(I,J,bi,bj) = 0. _d 0 |
63 |
|
aC2d(I,J,bi,bj) = 0. _d 0 |
64 |
pW(I,J,bi,bj) = 0. _d 0 |
pW(I,J,bi,bj) = 0. _d 0 |
65 |
pS(I,J,bi,bj) = 0. _d 0 |
pS(I,J,bi,bj) = 0. _d 0 |
66 |
pC(I,J,bi,bj) = 0. _d 0 |
pC(I,J,bi,bj) = 0. _d 0 |
91 |
DO K=1,Nr |
DO K=1,Nr |
92 |
DO J=1,sNy |
DO J=1,sNy |
93 |
DO I=1,sNx |
DO I=1,sNx |
94 |
|
C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect |
95 |
faceArea = _dyG(I,J,bi,bj)*drF(K) |
faceArea = _dyG(I,J,bi,bj)*drF(K) |
96 |
& *_hFacW(I,J,K,bi,bj) |
& *_hFacW(I,J,K,bi,bj) |
97 |
aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) + |
aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) |
98 |
& implicSurfPress*implicDiv2DFlow |
& + implicSurfPress*implicDiv2DFlow |
99 |
& *faceArea*recip_dxC(I,J,bi,bj) |
& *faceArea*recip_dxC(I,J,bi,bj) |
100 |
faceArea = _dxG(I,J,bi,bj)*drF(K) |
faceArea = _dxG(I,J,bi,bj)*drF(K) |
101 |
& *_hFacS(I,J,K,bi,bj) |
& *_hFacS(I,J,K,bi,bj) |
102 |
aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) + |
aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) |
103 |
& implicSurfPress*implicDiv2DFlow |
& + implicSurfPress*implicDiv2DFlow |
104 |
& *faceArea*recip_dyC(I,J,bi,bj) |
& *faceArea*recip_dyC(I,J,bi,bj) |
105 |
ENDDO |
ENDDO |
106 |
ENDDO |
ENDDO |
107 |
ENDDO |
ENDDO |
145 |
ENDDO |
ENDDO |
146 |
ENDDO |
ENDDO |
147 |
ENDDO |
ENDDO |
148 |
|
|
149 |
C-- Update overlap regions |
C-- Update overlap regions |
150 |
CcnhDebugStarts |
CcnhDebugStarts |
151 |
C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid ) |
C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid ) |
188 |
CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) |
CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) |
189 |
_END_MASTER( myThid ) |
_END_MASTER( myThid ) |
190 |
CcnhDebugEnds |
CcnhDebugEnds |
191 |
|
|
192 |
C-- Initialise preconditioner |
C-- Initialise preconditioner |
193 |
C Note. 20th May 1998 |
C Note. 20th May 1998 |
194 |
C I made a weird discovery! In the model paper we argue |
C I made a weird discovery! In the model paper we argue |
207 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
DO bi=myBxLo(myThid),myBxHi(myThid) |
208 |
DO J=1,sNy |
DO J=1,sNy |
209 |
DO I=1,sNx |
DO I=1,sNx |
210 |
|
ks = ksurfC(I,J,bi,bj) |
211 |
pC(I,J,bi,bj) = 1. _d 0 |
pC(I,J,bi,bj) = 1. _d 0 |
212 |
aC = -( |
aC = -( |
213 |
& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj) |
& aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj) |
214 |
& +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj) |
& +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj) |
215 |
& +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)* |
& +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks) |
216 |
& rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf |
& *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf |
217 |
& ) |
& ) |
218 |
aCs = -( |
aCs = -( |
219 |
& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj) |
& aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj) |
220 |
& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj) |
& +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj) |
221 |
& +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)* |
& +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*deepFac2F(ks) |
222 |
& rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf |
& *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf |
223 |
& ) |
& ) |
224 |
aCw = -( |
aCw = -( |
225 |
& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj) |
& aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj) |
226 |
& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj) |
& +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj) |
227 |
& +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)* |
& +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*deepFac2F(ks) |
228 |
& rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf |
& *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf |
229 |
& ) |
& ) |
230 |
IF ( aC .EQ. 0. ) THEN |
IF ( aC .EQ. 0. ) THEN |
231 |
pC(I,J,bi,bj) = 1. _d 0 |
pC(I,J,bi,bj) = 1. _d 0 |
235 |
IF ( aC + aCw .EQ. 0. ) THEN |
IF ( aC + aCw .EQ. 0. ) THEN |
236 |
pW(I,J,bi,bj) = 0. |
pW(I,J,bi,bj) = 0. |
237 |
ELSE |
ELSE |
238 |
pW(I,J,bi,bj) = |
pW(I,J,bi,bj) = |
239 |
& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 ) |
& -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 ) |
240 |
ENDIF |
ENDIF |
241 |
IF ( aC + aCs .EQ. 0. ) THEN |
IF ( aC + aCs .EQ. 0. ) THEN |
244 |
pS(I,J,bi,bj) = |
pS(I,J,bi,bj) = |
245 |
& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 ) |
& -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 ) |
246 |
ENDIF |
ENDIF |
247 |
|
C- store solver main diagonal : |
248 |
|
aC2d(I,J,bi,bj) = aC |
249 |
C pC(I,J,bi,bj) = 1. |
C pC(I,J,bi,bj) = 1. |
250 |
C pW(I,J,bi,bj) = 0. |
C pW(I,J,bi,bj) = 0. |
251 |
C pS(I,J,bi,bj) = 0. |
C pS(I,J,bi,bj) = 0. |