/[MITgcm]/MITgcm/pkg/seaice/seaice_jfnk.F
ViewVC logotype

Contents of /MITgcm/pkg/seaice/seaice_jfnk.F

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


Revision 1.27 - (show annotations) (download)
Mon Oct 20 03:20:57 2014 UTC (10 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65g
Changes since 1.26: +4 -1 lines
- ECCO_OPTIONS.h is needed when including ecco_cost.h, ecco.h
- AUTODIFF_OPTIONS.h is needed when including tamc.h, tamc_keys.h
- CTRL_OPTIONS.h is needed when including ctrl.h, etc

- pkg/seaice/seaice_cost*.F : clean up CPP brackets
- SEAICE_SIZE.h : replace ALLOW_AUTODIFF_TAMC with ALLOW_AUTODIFF to
  avoid needing AUTODIFF_OPTIONS.h anytime SEAICE_SIZE.h is included
  (it seems that THSICE_SIZE.h, PTRACERS_SIZE.h have the same issue...)

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.26 2014/03/20 09:24:49 mlosch Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_AUTODIFF
6 # include "AUTODIFF_OPTIONS.h"
7 #endif
8
9 C-- File seaice_jfnk.F: seaice jfnk dynamical solver S/R:
10 C-- Contents
11 C-- o SEAICE_JFNK
12 C-- o SEAICE_JFNK_UPDATE
13
14 CBOP
15 C !ROUTINE: SEAICE_JFNK
16 C !INTERFACE:
17 SUBROUTINE SEAICE_JFNK( myTime, myIter, myThid )
18
19 C !DESCRIPTION: \bv
20 C *==========================================================*
21 C | SUBROUTINE SEAICE_JFNK
22 C | o Ice dynamics using a Jacobian-free Newton-Krylov solver
23 C | following J.-F. Lemieux et al. Improving the numerical
24 C | convergence of viscous-plastic sea ice models with the
25 C | Jacobian-free Newton-Krylov method. J. Comp. Phys. 229,
26 C | 2840-2852 (2010).
27 C | o The logic follows JFs code.
28 C *==========================================================*
29 C | written by Martin Losch, Oct 2012
30 C *==========================================================*
31 C \ev
32
33 C !USES:
34 IMPLICIT NONE
35
36 C === Global variables ===
37 #include "SIZE.h"
38 #include "EEPARAMS.h"
39 #include "PARAMS.h"
40 #include "DYNVARS.h"
41 #include "GRID.h"
42 #include "SEAICE_SIZE.h"
43 #include "SEAICE_PARAMS.h"
44 #include "SEAICE.h"
45
46 #ifdef ALLOW_AUTODIFF_TAMC
47 # include "tamc.h"
48 #endif
49
50 C !INPUT/OUTPUT PARAMETERS:
51 C === Routine arguments ===
52 C myTime :: Simulation time
53 C myIter :: Simulation timestep number
54 C myThid :: my Thread Id. number
55 _RL myTime
56 INTEGER myIter
57 INTEGER myThid
58
59 #ifdef SEAICE_ALLOW_JFNK
60 C !FUNCTIONS:
61 LOGICAL DIFFERENT_MULTIPLE
62 EXTERNAL DIFFERENT_MULTIPLE
63
64 C !LOCAL VARIABLES:
65 C === Local variables ===
66 C i,j,bi,bj :: loop indices
67 INTEGER i,j,bi,bj
68 C loop indices
69 INTEGER newtonIter
70 INTEGER krylovIter, krylovFails
71 INTEGER totalKrylovItersLoc, totalNewtonItersLoc
72 C FGMRES flag that determines amount of output messages of fgmres
73 INTEGER iOutFGMRES
74 C FGMRES flag that indicates what fgmres wants us to do next
75 INTEGER iCode
76 _RL JFNKresidual
77 _RL JFNKresidualKm1
78 C parameters to compute convergence criterion
79 _RL JFNKgamma_lin
80 _RL FGMRESeps
81 _RL JFNKtol
82 C backward differences extrapolation factors
83 _RL bdfFac, bdfAlpha
84 C
85 _RL recip_deltaT
86 LOGICAL JFNKconverged, krylovConverged
87 LOGICAL writeNow
88 CHARACTER*(MAX_LEN_MBUF) msgBuf
89
90 C u/vIceRes :: residual of sea-ice momentum equations
91 _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
92 _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
93 C extra time level required for backward difference time stepping
94 _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
95 _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
96 C du/vIce :: ice velocity increment to be added to u/vIce
97 _RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
98 _RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
99 C precomputed (= constant per Newton iteration) versions of
100 C zeta, eta, and DWATN, press
101 _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
102 _RL etaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
103 _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
104 _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
105 CEOP
106
107 C Initialise
108 newtonIter = 0
109 krylovFails = 0
110 totalKrylovItersLoc = 0
111 JFNKconverged = .FALSE.
112 JFNKtol = 0. _d 0
113 JFNKresidual = 0. _d 0
114 JFNKresidualKm1 = 0. _d 0
115 FGMRESeps = 0. _d 0
116 recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn
117
118 iOutFGMRES=0
119 C with iOutFgmres=1, seaice_fgmres prints the residual at each iteration
120 IF ( debugLevel.GE.debLevC .AND.
121 & DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
122 & iOutFGMRES=1
123
124 C backward difference extrapolation factors
125 bdfFac = 0. _d 0
126 IF ( SEAICEuseBDF2 ) THEN
127 IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
128 bdfFac = 0. _d 0
129 ELSE
130 bdfFac = 0.5 _d 0
131 ENDIF
132 ENDIF
133 bdfAlpha = 1. _d 0 + bdfFac
134
135 DO bj=myByLo(myThid),myByHi(myThid)
136 DO bi=myBxLo(myThid),myBxHi(myThid)
137 DO J=1-OLy,sNy+OLy
138 DO I=1-OLx,sNx+OLx
139 uIceRes(I,J,bi,bj) = 0. _d 0
140 vIceRes(I,J,bi,bj) = 0. _d 0
141 duIce (I,J,bi,bj) = 0. _d 0
142 dvIce (I,J,bi,bj) = 0. _d 0
143 ENDDO
144 ENDDO
145 C cycle ice velocities
146 DO J=1-OLy,sNy+OLy
147 DO I=1-OLx,sNx+OLx
148 duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha
149 & + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac
150 dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha
151 & + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac
152 uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
153 vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
154 ENDDO
155 ENDDO
156 C As long as IMEX is not properly implemented leave this commented out
157 CML IF ( .NOT.SEAICEuseIMEX ) THEN
158 C Compute things that do no change during the Newton iteration:
159 C sea-surface tilt and wind stress:
160 C FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT
161 DO J=1-OLy,sNy+OLy
162 DO I=1-OLx,sNx+OLx
163 FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj)
164 & + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT
165 FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
166 & + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
167 ENDDO
168 ENDDO
169 CML ENDIF
170 ENDDO
171 ENDDO
172 C Start nonlinear Newton iteration: outer loop iteration
173 DO WHILE ( newtonIter.LT.SEAICEnewtonIterMax .AND.
174 & .NOT.JFNKconverged )
175 newtonIter = newtonIter + 1
176 C Compute initial residual F(u), (includes computation of global
177 C variables DWATN, zeta, and eta)
178 IF ( newtonIter .EQ. 1 ) CALL SEAICE_JFNK_UPDATE(
179 I duIce, dvIce,
180 U uIce, vIce, JFNKresidual,
181 O uIceRes, vIceRes,
182 I newtonIter, myTime, myIter, myThid )
183 C local copies of precomputed coefficients that are to stay
184 C constant for the preconditioner
185 DO bj=myByLo(myThid),myByHi(myThid)
186 DO bi=myBxLo(myThid),myBxHi(myThid)
187 DO j=1-OLy,sNy+OLy
188 DO i=1-OLx,sNx+OLx
189 zetaPre(I,J,bi,bj) = zeta(I,J,bi,bj)
190 etaPre(I,J,bi,bj) = eta(I,J,bi,bj)
191 etaZPre(I,J,bi,bj) = etaZ(I,J,bi,bj)
192 dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj)
193 ENDDO
194 ENDDO
195 ENDDO
196 ENDDO
197 C compute convergence criterion for linear preconditioned FGMRES
198 JFNKgamma_lin = JFNKgamma_lin_max
199 IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter
200 & .AND.JFNKresidual.LT.JFNKres_t ) THEN
201 C Eisenstat and Walker (1996), eq.(2.6)
202 JFNKgamma_lin = SEAICE_JFNKphi
203 & *( JFNKresidual/JFNKresidualKm1 )**SEAICE_JFNKalpha
204 JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin)
205 JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin)
206 ENDIF
207 C save the residual for the next iteration
208 JFNKresidualKm1 = JFNKresidual
209
210 C The Krylov iteration using FGMRES, the preconditioner is LSOR
211 C for now. The code is adapted from SEAICE_LSR, but heavily stripped
212 C down.
213 C krylovIter is mapped into "its" in seaice_fgmres and is incremented
214 C in that routine
215 krylovIter = 0
216 iCode = 0
217
218 JFNKconverged = JFNKresidual.LT.JFNKtol
219
220 C do Krylov loop only if convergence is not reached
221
222 IF ( .NOT.JFNKconverged ) THEN
223
224 C start Krylov iteration (FGMRES)
225
226 krylovConverged = .FALSE.
227 FGMRESeps = JFNKgamma_lin * JFNKresidual
228 DO WHILE ( .NOT.krylovConverged )
229 C solution vector sol = du/vIce
230 C residual vector (rhs) Fu = u/vIceRes
231 C output work vectors wk1, -> input work vector wk2
232
233 CALL SEAICE_FGMRES_DRIVER(
234 I uIceRes, vIceRes,
235 U duIce, dvIce, iCode,
236 I FGMRESeps, iOutFGMRES,
237 I newtonIter, krylovIter, myTime, myIter, myThid )
238 C FGMRES returns iCode either asking for an new preconditioned vector
239 C or product of matrix (Jacobian) times vector. For iCode = 0, terminate
240 C iteration
241 IF (iCode.EQ.1) THEN
242 C Call preconditioner
243 IF ( SOLV_MAX_ITERS .GT. 0 )
244 & CALL SEAICE_PRECONDITIONER(
245 U duIce, dvIce,
246 I zetaPre, etaPre, etaZpre, dwatPre,
247 I newtonIter, krylovIter, myTime, myIter, myThid )
248 ELSEIF (iCode.GE.2) THEN
249 C Compute Jacobian times vector
250 CALL SEAICE_JACVEC(
251 I uIce, vIce, uIceRes, vIceRes,
252 U duIce, dvIce,
253 I newtonIter, krylovIter, myTime, myIter, myThid )
254 ENDIF
255 krylovConverged = iCode.EQ.0
256 C End of Krylov iterate
257 ENDDO
258 totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
259 C some output diagnostics
260 IF ( debugLevel.GE.debLevA ) THEN
261 _BEGIN_MASTER( myThid )
262 totalNewtonItersLoc =
263 & SEAICEnewtonIterMax*(myIter-nIter0)+newtonIter
264 WRITE(msgBuf,'(2A,2(1XI6),2E12.5)')
265 & ' S/R SEAICE_JFNK: Newton iterate / total, ',
266 & 'JFNKgamma_lin, initial norm = ',
267 & newtonIter, totalNewtonItersLoc,
268 & JFNKgamma_lin,JFNKresidual
269 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
270 & SQUEEZE_RIGHT, myThid )
271 WRITE(msgBuf,'(3(A,I6))')
272 & ' S/R SEAICE_JFNK: Newton iterate / total = ',newtonIter,
273 & ' / ', totalNewtonItersLoc,
274 & ', Nb. of FGMRES iterations = ', krylovIter
275 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
276 & SQUEEZE_RIGHT, myThid )
277 _END_MASTER( myThid )
278 ENDIF
279 IF ( krylovIter.EQ.SEAICEkrylovIterMax ) THEN
280 krylovFails = krylovFails + 1
281 ENDIF
282 C Set the stopping criterion for the Newton iteration and the
283 C criterion for the transition from accurate to approximate FGMRES
284 IF ( newtonIter .EQ. 1 ) THEN
285 JFNKtol=JFNKgamma_nonlin*JFNKresidual
286 IF ( JFNKres_tFac .NE. UNSET_RL )
287 & JFNKres_t = JFNKresidual * JFNKres_tFac
288 ENDIF
289 C Update linear solution vector and return to Newton iteration
290 C Do a linesearch if necessary, and compute a new residual.
291 C Note that it should be possible to do the following operations
292 C at the beginning of the Newton iteration, thereby saving us from
293 C the extra call of seaice_jfnk_update, but unfortunately that
294 C changes the results, so we leave the stuff here for now.
295 CALL SEAICE_JFNK_UPDATE(
296 I duIce, dvIce,
297 U uIce, vIce, JFNKresidual,
298 O uIceRes, vIceRes,
299 I newtonIter, myTime, myIter, myThid )
300 C reset du/vIce here instead of setting sol = 0 in seaice_fgmres_driver
301 DO bj=myByLo(myThid),myByHi(myThid)
302 DO bi=myBxLo(myThid),myBxHi(myThid)
303 DO J=1-OLy,sNy+OLy
304 DO I=1-OLx,sNx+OLx
305 duIce(I,J,bi,bj)= 0. _d 0
306 dvIce(I,J,bi,bj)= 0. _d 0
307 ENDDO
308 ENDDO
309 ENDDO
310 ENDDO
311 ENDIF
312 C end of Newton iterate
313 ENDDO
314
315 C-- Output diagnostics
316
317 IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
318 C Count iterations
319 totalJFNKtimeSteps = totalJFNKtimeSteps + 1
320 totalNewtonIters = totalNewtonIters + newtonIter
321 totalKrylovIters = totalKrylovIters + totalKrylovItersLoc
322 C Record failure
323 totalKrylovFails = totalKrylovFails + krylovFails
324 IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
325 totalNewtonFails = totalNewtonFails + 1
326 ENDIF
327 ENDIF
328 C Decide whether it is time to dump and reset the counter
329 writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
330 & myTime+deltaTClock, deltaTClock)
331 #ifdef ALLOW_CAL
332 IF ( useCAL ) THEN
333 CALL CAL_TIME2DUMP(
334 I zeroRL, SEAICE_monFreq, deltaTClock,
335 U writeNow,
336 I myTime+deltaTclock, myIter+1, myThid )
337 ENDIF
338 #endif
339 IF ( writeNow ) THEN
340 _BEGIN_MASTER( myThid )
341 WRITE(msgBuf,'(A)')
342 &' // ======================================================='
343 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
344 & SQUEEZE_RIGHT, myThid )
345 WRITE(msgBuf,'(A)') ' // Begin JFNK statistics'
346 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
347 & SQUEEZE_RIGHT, myThid )
348 WRITE(msgBuf,'(A)')
349 &' // ======================================================='
350 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
351 & SQUEEZE_RIGHT, myThid )
352 WRITE(msgBuf,'(A,I10)')
353 & ' %JFNK_MON: time step = ', myIter+1
354 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
355 & SQUEEZE_RIGHT, myThid )
356 WRITE(msgBuf,'(A,I10)')
357 & ' %JFNK_MON: Nb. of time steps = ', totalJFNKtimeSteps
358 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
359 & SQUEEZE_RIGHT, myThid )
360 WRITE(msgBuf,'(A,I10)')
361 & ' %JFNK_MON: Nb. of Newton steps = ', totalNewtonIters
362 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
363 & SQUEEZE_RIGHT, myThid )
364 WRITE(msgBuf,'(A,I10)')
365 & ' %JFNK_MON: Nb. of Krylov steps = ', totalKrylovIters
366 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
367 & SQUEEZE_RIGHT, myThid )
368 WRITE(msgBuf,'(A,I10)')
369 & ' %JFNK_MON: Nb. of Newton failures = ', totalNewtonFails
370 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
371 & SQUEEZE_RIGHT, myThid )
372 WRITE(msgBuf,'(A,I10)')
373 & ' %JFNK_MON: Nb. of Krylov failures = ', totalKrylovFails
374 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
375 & SQUEEZE_RIGHT, myThid )
376 WRITE(msgBuf,'(A)')
377 &' // ======================================================='
378 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
379 & SQUEEZE_RIGHT, myThid )
380 WRITE(msgBuf,'(A)') ' // End JFNK statistics'
381 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
382 & SQUEEZE_RIGHT, myThid )
383 WRITE(msgBuf,'(A)')
384 &' // ======================================================='
385 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
386 & SQUEEZE_RIGHT, myThid )
387 _END_MASTER( myThid )
388 C reset and start again
389 totalJFNKtimeSteps = 0
390 totalNewtonIters = 0
391 totalKrylovIters = 0
392 totalKrylovFails = 0
393 totalNewtonFails = 0
394 ENDIF
395
396 C Print more debugging information
397 IF ( debugLevel.GE.debLevA ) THEN
398 IF ( newtonIter .EQ. SEAICEnewtonIterMax ) THEN
399 _BEGIN_MASTER( myThid )
400 WRITE(msgBuf,'(A,I10)')
401 & ' S/R SEAICE_JFNK: JFNK did not converge in timestep ',
402 & myIter+1
403 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
404 & SQUEEZE_RIGHT, myThid )
405 _END_MASTER( myThid )
406 ENDIF
407 IF ( krylovFails .GT. 0 ) THEN
408 _BEGIN_MASTER( myThid )
409 WRITE(msgBuf,'(A,I4,A,I10)')
410 & ' S/R SEAICE_JFNK: FGMRES did not converge ',
411 & krylovFails, ' times in timestep ', myIter+1
412 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
413 & SQUEEZE_RIGHT, myThid )
414 _END_MASTER( myThid )
415 ENDIF
416 _BEGIN_MASTER( myThid )
417 WRITE(msgBuf,'(A,I6,A,I10)')
418 & ' S/R SEAICE_JFNK: Total number FGMRES iterations = ',
419 & totalKrylovItersLoc, ' in timestep ', myIter+1
420 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
421 & SQUEEZE_RIGHT, myThid )
422 _END_MASTER( myThid )
423 ENDIF
424
425 RETURN
426 END
427
428 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
429 CBOP
430 C !ROUTINE: SEAICE_JFNK_UPDATE
431 C !INTERFACE:
432
433 SUBROUTINE SEAICE_JFNK_UPDATE(
434 I duIce, dvIce,
435 U uIce, vIce, JFNKresidual,
436 O uIceRes, vIceRes,
437 I newtonIter, myTime, myIter, myThid )
438
439 C !DESCRIPTION: \bv
440 C *==========================================================*
441 C | SUBROUTINE SEAICE_JFNK_UPDATE
442 C | o Update velocities with incremental solutions of FGMRES
443 C | o compute residual of updated solutions and do
444 C | o linesearch:
445 C | reduce update until residual is smaller than previous
446 C | one (input)
447 C *==========================================================*
448 C | written by Martin Losch, Jan 2013
449 C *==========================================================*
450 C \ev
451
452 C !USES:
453 IMPLICIT NONE
454
455 C === Global variables ===
456 #include "SIZE.h"
457 #include "EEPARAMS.h"
458 #include "PARAMS.h"
459 #include "SEAICE_SIZE.h"
460 #include "SEAICE_PARAMS.h"
461
462 C !INPUT/OUTPUT PARAMETERS:
463 C === Routine arguments ===
464 C myTime :: Simulation time
465 C myIter :: Simulation timestep number
466 C myThid :: my Thread Id. number
467 C newtonIter :: current iterate of Newton iteration
468 _RL myTime
469 INTEGER myIter
470 INTEGER myThid
471 INTEGER newtonIter
472 C JFNKresidual :: Residual at the beginning of the FGMRES iteration,
473 C changes with newtonIter (updated)
474 _RL JFNKresidual
475 C du/vIce :: ice velocity increment to be added to u/vIce (input)
476 _RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
477 _RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
478 C u/vIce :: ice velocity increment to be added to u/vIce (updated)
479 _RL uIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
480 _RL vIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
481 C u/vIceRes :: residual of sea-ice momentum equations (output)
482 _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
483 _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
484
485 C !LOCAL VARIABLES:
486 C === Local variables ===
487 C i,j,bi,bj :: loop indices
488 INTEGER i,j,bi,bj
489 INTEGER l
490 _RL resLoc, facLS
491 LOGICAL doLineSearch
492 C nVec :: size of the input vector(s)
493 C resTmp :: vector version of the residuals
494 INTEGER nVec
495 PARAMETER ( nVec = 2*sNx*sNy )
496 _RL resTmp (nVec,1,nSx,nSy)
497
498 CHARACTER*(MAX_LEN_MBUF) msgBuf
499 CEOP
500
501 C Initialise some local variables
502 l = 0
503 resLoc = JFNKresidual
504 facLS = 1. _d 0
505 doLineSearch = .TRUE.
506 DO WHILE ( doLineSearch )
507 C Create update
508 DO bj=myByLo(myThid),myByHi(myThid)
509 DO bi=myBxLo(myThid),myBxHi(myThid)
510 DO J=1-OLy,sNy+OLy
511 DO I=1-OLx,sNx+OLx
512 uIce(I,J,bi,bj) = uIce(I,J,bi,bj)+facLS*duIce(I,J,bi,bj)
513 vIce(I,J,bi,bj) = vIce(I,J,bi,bj)+facLS*dvIce(I,J,bi,bj)
514 ENDDO
515 ENDDO
516 ENDDO
517 ENDDO
518 C Compute current residual F(u), (includes re-computation of global
519 C variables DWATN, zeta, and eta, i.e. they are different after this)
520 CALL SEAICE_CALC_RESIDUAL(
521 I uIce, vIce,
522 O uIceRes, vIceRes,
523 I newtonIter, 0, myTime, myIter, myThid )
524 C Important: Compute the norm of the residual using the same scalar
525 C product that SEAICE_FGMRES does
526 CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,resTmp,.TRUE.,myThid)
527 CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,resLoc,myThid)
528 resLoc = SQRT(resLoc)
529 C Determine, if we need more iterations
530 doLineSearch = resLoc .GE. JFNKresidual
531 C Limit the maximum number of iterations arbitrarily to four
532 doLineSearch = doLineSearch .AND. l .LT. 4
533 C For the first iteration du/vIce = 0 and there will be no
534 C improvement of the residual possible, so we do only the first
535 C iteration
536 IF ( newtonIter .EQ. 1 ) doLineSearch = .FALSE.
537 C Only start a linesearch after some Newton iterations
538 IF ( newtonIter .LE. SEAICE_JFNK_lsIter ) doLineSearch = .FALSE.
539 C Increment counter
540 l = l + 1
541 C some output diagnostics
542 IF ( debugLevel.GE.debLevA .AND. doLineSearch ) THEN
543 _BEGIN_MASTER( myThid )
544 WRITE(msgBuf,'(2A,2(1XI6),3E12.5)')
545 & ' S/R SEAICE_JFNK_UPDATE: Newton iter, LSiter, ',
546 & 'facLS, JFNKresidual, resLoc = ',
547 & newtonIter, l, facLS, JFNKresidual, resLoc
548 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
549 & SQUEEZE_RIGHT, myThid )
550 _END_MASTER( myThid )
551 ENDIF
552 C Get ready for the next iteration: after adding du/vIce in the first
553 C iteration, we substract 0.5*du/vIce from u/vIce in the next
554 C iterations, 0.25*du/vIce in the second, etc.
555 facLS = - 0.5 _d 0 * ABS(facLS)
556 ENDDO
557 C This is the new residual
558 JFNKresidual = resLoc
559
560 #endif /* SEAICE_ALLOW_JFNK */
561
562 RETURN
563 END

  ViewVC Help
Powered by ViewVC 1.1.22