/[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.42 - (hide 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 edhill 1.42 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 jmc 1.30 C $Name: $
3 cnh 1.1
4 edhill 1.42 #include "PACKAGES_CONFIG.h"
5 adcroft 1.13 #include "CPP_OPTIONS.h"
6 cnh 1.1
7 cnh 1.37 CBOP
8     C !ROUTINE: INI_CG2D
9     C !INTERFACE:
10 cnh 1.1 SUBROUTINE INI_CG2D( myThid )
11 cnh 1.37
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 adcroft 1.23 IMPLICIT NONE
24 cnh 1.1 C === Global variables ===
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "PARAMS.h"
28     #include "GRID.h"
29 adcroft 1.34 c#include "DYNVARS.h"
30 jmc 1.31 #include "SURFACE.h"
31 adcroft 1.34 #include "CG2D.h"
32 adcroft 1.26 #ifdef ALLOW_OBCS
33 adcroft 1.22 #include "OBCS.h"
34 adcroft 1.26 #endif
35 cnh 1.1
36 cnh 1.37 C !INPUT/OUTPUT PARAMETERS:
37 cnh 1.1 C === Routine arguments ===
38     C myThid - Thread no. that called this routine.
39     INTEGER myThid
40    
41 cnh 1.37 C !LOCAL VARIABLES:
42 cnh 1.1 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 adcroft 1.34 C myNorm - Work variable used in calculating normalisation factor
50     C sumArea - Work variable used to compute the total Domain Area
51 cnh 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
52     INTEGER bi, bj
53     INTEGER I, J, K
54 cnh 1.7 _RL faceArea
55 cnh 1.15 _RS myNorm
56 adcroft 1.34 _RL sumArea
57 cnh 1.4 _RL aC, aCw, aCs
58 cnh 1.37 CEOP
59 jmc 1.38
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 jmc 1.31
83 cnh 1.1 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 cnh 1.17 DO K=1,Nr
96 cnh 1.1 DO J=1,sNy
97     DO I=1,sNx
98 cnh 1.20 faceArea = _dyG(I,J,bi,bj)*drF(K)
99     & *_hFacW(I,J,K,bi,bj)
100 cnh 1.1 aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
101 jmc 1.32 & implicSurfPress*implicDiv2DFlow
102     & *faceArea*recip_dxC(I,J,bi,bj)
103 cnh 1.20 faceArea = _dxG(I,J,bi,bj)*drF(K)
104     & *_hFacS(I,J,K,bi,bj)
105 cnh 1.1 aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
106 jmc 1.32 & implicSurfPress*implicDiv2DFlow
107     & *faceArea*recip_dyC(I,J,bi,bj)
108 cnh 1.1 ENDDO
109     ENDDO
110     ENDDO
111 adcroft 1.26 #ifdef ALLOW_OBCS
112 adcroft 1.28 IF (useOBCS) THEN
113 adcroft 1.22 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 adcroft 1.26 #endif
127 cnh 1.1 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 adcroft 1.25 _GLOBAL_MAX_R4( myNorm, myThid )
136     IF ( myNorm .NE. 0. _d 0 ) THEN
137     myNorm = 1. _d 0/myNorm
138 cnh 1.1 ELSE
139     myNorm = 1. _d 0
140     ENDIF
141 cnh 1.5 cg2dNorm = myNorm
142 cnh 1.1 _BEGIN_MASTER( myThid )
143     CcnhDebugStarts
144 cnh 1.20 WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
145     & cg2dNorm
146 cnh 1.3 CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
147     WRITE(msgBuf,*) ' '
148 cnh 1.1 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 cnh 1.14 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 cnh 1.1 CcnhDebugEnds
167 adcroft 1.34 c _EXCH_XY_R4(aW2d, myThid)
168     c _EXCH_XY_R4(aS2d, myThid)
169     CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
170 cnh 1.1 CcnhDebugStarts
171 adcroft 1.24 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 cnh 1.1 CcnhDebugEnds
174    
175 adcroft 1.34 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 adcroft 1.35 c WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea
195 adcroft 1.34 _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 jmc 1.41 & * sumArea / deltaTmom
199 jmc 1.40 _BEGIN_MASTER( myThid )
200     WRITE(standardMessageUnit,'(2A,1P2E22.14)') ' ini_cg2d: ',
201 adcroft 1.34 & 'sumArea,cg2dTolerance =', sumArea,cg2dTolerance
202 jmc 1.40 _END_MASTER( myThid )
203 adcroft 1.34 ENDIF
204    
205 cnh 1.1 C-- Initialise preconditioner
206 cnh 1.4 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 cnh 1.1 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 cnh 1.4 aC = -(
225     & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
226 adcroft 1.24 & +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
227 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
228 adcroft 1.39 & rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
229 cnh 1.4 & )
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 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
234 adcroft 1.39 & rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
235 cnh 1.4 & )
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 jmc 1.32 & +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
240 adcroft 1.39 & rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
241 cnh 1.4 & )
242     IF ( aC .EQ. 0. ) THEN
243 adcroft 1.27 pC(I,J,bi,bj) = 1. _d 0
244 cnh 1.4 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 cnh 1.6 C pC(I,J,bi,bj) = 1.
260     C pW(I,J,bi,bj) = 0.
261     C pS(I,J,bi,bj) = 0.
262 cnh 1.1 ENDDO
263     ENDDO
264     ENDDO
265     ENDDO
266     C-- Update overlap regions
267     _EXCH_XY_R4(pC, myThid)
268 adcroft 1.34 c _EXCH_XY_R4(pW, myThid)
269     c _EXCH_XY_R4(pS, myThid)
270     CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
271 cnh 1.18 CcnhDebugStarts
272 adcroft 1.24 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 cnh 1.18 CcnhDebugEnds
276 cnh 1.1
277     RETURN
278     END
  ViewVC Help
Powered by ViewVC 1.1.22