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

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

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

Revision 1.42 - (show annotations) (download)
Thu Oct 9 04:19:18 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint57g_pre, checkpoint57b_post, checkpoint52d_pre, checkpoint57g_post, checkpoint56b_post, checkpoint52j_pre, checkpoint51o_pre, checkpoint54d_post, checkpoint54e_post, checkpoint51l_post, checkpoint57d_post, checkpoint57i_post, checkpoint52l_post, checkpoint52k_post, checkpoint55, checkpoint54, checkpoint57, checkpoint56, checkpoint53, checkpoint52, checkpoint52f_post, checkpoint54f_post, checkpoint51t_post, checkpoint51n_post, checkpoint55i_post, checkpoint57l_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55c_post, checkpoint52e_pre, checkpoint57f_post, checkpoint52e_post, checkpoint51n_pre, checkpoint53d_post, checkpoint57a_post, checkpoint57h_pre, checkpoint52b_pre, checkpoint54b_post, checkpoint57h_post, checkpoint51l_pre, checkpoint52m_post, checkpoint55g_post, checkpoint51q_post, checkpoint52b_post, checkpoint52c_post, checkpoint57c_post, checkpoint52f_pre, checkpoint55d_post, checkpoint54a_pre, checkpoint53c_post, checkpoint55d_pre, checkpoint57c_pre, checkpoint55j_post, checkpoint54a_post, checkpoint55h_post, checkpoint51r_post, checkpoint51i_post, checkpoint57e_post, checkpoint55b_post, checkpoint53a_post, checkpoint55f_post, checkpoint52d_post, checkpoint53g_post, eckpoint57e_pre, checkpoint52a_pre, checkpoint52i_post, checkpoint52h_pre, checkpoint56a_post, checkpoint53f_post, checkpoint57h_done, checkpoint52j_post, checkpoint57j_post, checkpoint57f_pre, branch-netcdf, checkpoint52n_post, checkpoint53b_pre, checkpoint56c_post, checkpoint57a_pre, checkpoint55a_post, checkpoint51o_post, checkpoint57k_post, checkpoint53b_post, checkpoint52a_post, ecco_c52_e35, checkpoint51m_post, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.41: +2 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

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