/[MITgcm]/MITgcm/model/src/ini_cg2d.F
ViewVC logotype

Annotation of /MITgcm/model/src/ini_cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.39 - (hide annotations) (download)
Fri Jun 21 18:36:05 2002 UTC (21 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint46l_post, checkpoint46g_pre, checkpoint46f_post, checkpoint46b_post, checkpoint46l_pre, checkpoint46d_pre, checkpoint45d_post, checkpoint46j_pre, checkpoint46a_post, checkpoint46j_post, checkpoint46k_post, checkpoint46e_pre, checkpoint46b_pre, checkpoint46c_pre, checkpoint46, checkpoint46h_pre, checkpoint46a_pre, checkpoint46g_post, checkpoint46i_post, checkpoint46c_post, checkpoint46e_post, checkpoint46h_post, checkpoint46d_post
Changes since 1.38: +5 -5 lines 
Added new parameter: deltaTfreesurf

Previously, the free-surface equation was intergrated forward
synchronously with the momentum equations. It is more consistent
to use the tracer time-step. This increases the number of
iterations required but strengthens the damping.

We *SHOULD* make the default time-step equal to the tracer time-step.
However, we don't for backward compatibility. At some point in the
future we need to change the default behaviour.

It turns out that the reason for the "reduced stability" encountered
in large-scale runs seems to be related to excess variability in
the free surface which in turn happens when the waves aren't damped.
Using a longer time-step fixes this.
1 adcroft 1.39 C $Header: /u/gcmpack/MITgcm/model/src/ini_cg2d.F,v 1.38 2002/02/10 20:04:11 jmc Exp $
2 jmc 1.30 C $Name: $
3 cnh 1.1
4 adcroft 1.13 #include "CPP_OPTIONS.h"
5 cnh 1.1
6 cnh 1.37 CBOP
7     C !ROUTINE: INI_CG2D
8     C !INTERFACE:
9 cnh 1.1 SUBROUTINE INI_CG2D( myThid )
10 cnh 1.37
11     C !DESCRIPTION: \bv
12     C *==========================================================*
13     C | SUBROUTINE INI_CG2D
14     C | o Initialise 2d conjugate gradient solver operators.
15     C *==========================================================*
16     C | These arrays are purely a function of the basin geom.
17     C | We set then here once and them use then repeatedly.
18     C *==========================================================*
19     C \ev
20    
21     C !USES:
22 adcroft 1.23 IMPLICIT NONE
23 cnh 1.1 C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28 adcroft 1.34 c#include "DYNVARS.h"
29 jmc 1.31 #include "SURFACE.h"
30 adcroft 1.34 #include "CG2D.h"
31 adcroft 1.26 #ifdef ALLOW_OBCS
32 adcroft 1.22 #include "OBCS.h"
33 adcroft 1.26 #endif
34 cnh 1.1
35 cnh 1.37 C !INPUT/OUTPUT PARAMETERS:
36 cnh 1.1 C === Routine arguments ===
37     C myThid - Thread no. that called this routine.
38     INTEGER myThid
39    
40 cnh 1.37 C !LOCAL VARIABLES:
41 cnh 1.1 C === Local variables ===
42     C xG, yG - Global coordinate location.
43     C zG
44     C iG, jG - Global coordinate index
45     C bi,bj - Loop counters
46     C faceArea - Temporary used to hold cell face areas.
47     C I,J,K
48 adcroft 1.34 C myNorm - Work variable used in calculating normalisation factor
49     C sumArea - Work variable used to compute the total Domain Area
50 cnh 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
51     INTEGER bi, bj
52     INTEGER I, J, K
53 cnh 1.7 _RL faceArea
54 cnh 1.15 _RS myNorm
55 adcroft 1.34 _RL sumArea
56 cnh 1.4 _RL aC, aCw, aCs
57 cnh 1.37 CEOP
58 jmc 1.38
59     C-- Initialize arrays in common blocs (CG2D.h) ; not really necessary
60     C but safer when EXCH do not fill all the overlap regions.
61     DO bj=myByLo(myThid),myByHi(myThid)
62     DO bi=myBxLo(myThid),myBxHi(myThid)
63     DO J=1-OLy,sNy+OLy
64     DO I=1-OLx,sNx+OLx
65     aW2d(I,J,bi,bj) = 0. _d 0
66     aS2d(I,J,bi,bj) = 0. _d 0
67     pW(I,J,bi,bj) = 0. _d 0
68     pS(I,J,bi,bj) = 0. _d 0
69     pC(I,J,bi,bj) = 0. _d 0
70     cg2d_q(I,J,bi,bj) = 0. _d 0
71     ENDDO
72     ENDDO
73     DO J=1-1,sNy+1
74     DO I=1-1,sNx+1
75     cg2d_r(I,J,bi,bj) = 0. _d 0
76     cg2d_s(I,J,bi,bj) = 0. _d 0
77     ENDDO
78     ENDDO
79     ENDDO
80     ENDDO
81 jmc 1.31
82 cnh 1.1 C-- Initialise laplace operator
83     C aW2d: integral in Z Ax/dX
84     C aS2d: integral in Z Ay/dY
85     myNorm = 0. _d 0
86     DO bj=myByLo(myThid),myByHi(myThid)
87     DO bi=myBxLo(myThid),myBxHi(myThid)
88     DO J=1,sNy
89     DO I=1,sNx
90     aW2d(I,J,bi,bj) = 0. _d 0
91     aS2d(I,J,bi,bj) = 0. _d 0
92     ENDDO
93     ENDDO
94 cnh 1.17 DO K=1,Nr
95 cnh 1.1 DO J=1,sNy
96     DO I=1,sNx
97 cnh 1.20 faceArea = _dyG(I,J,bi,bj)*drF(K)
98     & *_hFacW(I,J,K,bi,bj)
99 cnh 1.1 aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
100 jmc 1.32 & implicSurfPress*implicDiv2DFlow
101     & *faceArea*recip_dxC(I,J,bi,bj)
102 cnh 1.20 faceArea = _dxG(I,J,bi,bj)*drF(K)
103     & *_hFacS(I,J,K,bi,bj)
104 cnh 1.1 aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
105 jmc 1.32 & implicSurfPress*implicDiv2DFlow
106     & *faceArea*recip_dyC(I,J,bi,bj)
107 cnh 1.1 ENDDO
108     ENDDO
109     ENDDO
110 adcroft 1.26 #ifdef ALLOW_OBCS
111 adcroft 1.28 IF (useOBCS) THEN
112 adcroft 1.22 DO I=1,sNx
113     IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
114     IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
115     IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
116     IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
117     ENDDO
118     DO J=1,sNy
119     IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
120     IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
121     IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
122     IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
123     ENDDO
124     ENDIF
125 adcroft 1.26 #endif
126 cnh 1.1 DO J=1,sNy
127     DO I=1,sNx
128     myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
129     myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)
130     ENDDO
131     ENDDO
132     ENDDO
133     ENDDO
134 adcroft 1.25 _GLOBAL_MAX_R4( myNorm, myThid )
135     IF ( myNorm .NE. 0. _d 0 ) THEN
136     myNorm = 1. _d 0/myNorm
137 cnh 1.1 ELSE
138     myNorm = 1. _d 0
139     ENDIF
140 cnh 1.5 cg2dNorm = myNorm
141 cnh 1.1 _BEGIN_MASTER( myThid )
142     CcnhDebugStarts
143 cnh 1.20 WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
144     & cg2dNorm
145 cnh 1.3 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
146     WRITE(msgBuf,*) ' '
147 cnh 1.1 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
148     CcnhDebugEnds
149     _END_MASTER( myThid )
150     DO bj=myByLo(myThid),myByHi(myThid)
151     DO bi=myBxLo(myThid),myBxHi(myThid)
152     DO J=1,sNy
153     DO I=1,sNx
154     aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm
155     aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm
156     ENDDO
157     ENDDO
158     ENDDO
159     ENDDO
160    
161     C-- Update overlap regions
162     CcnhDebugStarts
163 cnh 1.14 C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
164     C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
165 cnh 1.1 CcnhDebugEnds
166 adcroft 1.34 c _EXCH_XY_R4(aW2d, myThid)
167     c _EXCH_XY_R4(aS2d, myThid)
168     CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
169 cnh 1.1 CcnhDebugStarts
170 adcroft 1.24 C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
171     C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
172 cnh 1.1 CcnhDebugEnds
173    
174 adcroft 1.34 C-- Define the solver tolerance in the appropriate Unit :
175     cg2dNormaliseRHS = cg2dTargetResWunit.LE.0
176     IF (cg2dNormaliseRHS) THEN
177     C- when using a normalisation of RHS, tolerance has no unit => no conversion
178     cg2dTolerance = cg2dTargetResidual
179     ELSE
180     C- compute the total Area of the domain :
181     sumArea = 0.
182     DO bj=myByLo(myThid),myByHi(myThid)
183     DO bi=myBxLo(myThid),myBxHi(myThid)
184     DO j=1,sNy
185     DO i=1,sNx
186     IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN
187     sumArea = sumArea + rA(i,j,bi,bj)
188     ENDIF
189     ENDDO
190     ENDDO
191     ENDDO
192     ENDDO
193 adcroft 1.35 c WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea
194 adcroft 1.34 _GLOBAL_SUM_R8( sumArea, myThid )
195     C- convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
196     cg2dTolerance = cg2dNorm * cg2dTargetResWunit
197 adcroft 1.39 & * sumArea / deltaTfreesurf
198 adcroft 1.35 WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ',
199 adcroft 1.34 & 'sumArea,cg2dTolerance =', sumArea,cg2dTolerance
200     ENDIF
201    
202 cnh 1.1 C-- Initialise preconditioner
203 cnh 1.4 C Note. 20th May 1998
204     C I made a weird discovery! In the model paper we argue
205     C for the form of the preconditioner used here ( see
206     C A Finite-volume, Incompressible Navier-Stokes Model
207     C ...., Marshall et. al ). The algebra gives a simple
208     C 0.5 factor for the averaging of ac and aCw to get a
209     C symmettric pre-conditioner. By using a factor of 0.51
210     C i.e. scaling the off-diagonal terms in the
211     C preconditioner down slightly I managed to get the
212     C number of iterations for convergence in a test case to
213     C drop form 192 -> 134! Need to investigate this further!
214     C For now I have introduced a parameter cg2dpcOffDFac which
215     C defaults to 0.51 but can be set at runtime.
216 cnh 1.1 DO bj=myByLo(myThid),myByHi(myThid)
217     DO bi=myBxLo(myThid),myBxHi(myThid)
218     DO J=1,sNy
219     DO I=1,sNx
220     pC(I,J,bi,bj) = 1. _d 0
221 cnh 1.4 aC = -(
222     & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
223 adcroft 1.24 & +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
224 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
225 adcroft 1.39 & rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
226 cnh 1.4 & )
227     aCs = -(
228     & aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
229     & +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
230 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
231 adcroft 1.39 & rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
232 cnh 1.4 & )
233     aCw = -(
234     & aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
235     & +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
236 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
237 adcroft 1.39 & rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
238 cnh 1.4 & )
239     IF ( aC .EQ. 0. ) THEN
240 adcroft 1.27 pC(I,J,bi,bj) = 1. _d 0
241 cnh 1.4 ELSE
242     pC(I,J,bi,bj) = 1. _d 0 / aC
243     ENDIF
244     IF ( aC + aCw .EQ. 0. ) THEN
245     pW(I,J,bi,bj) = 0.
246     ELSE
247     pW(I,J,bi,bj) =
248     & -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
249     ENDIF
250     IF ( aC + aCs .EQ. 0. ) THEN
251     pS(I,J,bi,bj) = 0.
252     ELSE
253     pS(I,J,bi,bj) =
254     & -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
255     ENDIF
256 cnh 1.6 C pC(I,J,bi,bj) = 1.
257     C pW(I,J,bi,bj) = 0.
258     C pS(I,J,bi,bj) = 0.
259 cnh 1.1 ENDDO
260     ENDDO
261     ENDDO
262     ENDDO
263     C-- Update overlap regions
264     _EXCH_XY_R4(pC, myThid)
265 adcroft 1.34 c _EXCH_XY_R4(pW, myThid)
266     c _EXCH_XY_R4(pS, myThid)
267     CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
268 cnh 1.18 CcnhDebugStarts
269 adcroft 1.24 C CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid )
270     C CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid )
271     C CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid )
272 cnh 1.18 CcnhDebugEnds
273 cnh 1.1
274     RETURN
275     END
  ViewVC Help
Powered by ViewVC 1.1.22