/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F
ViewVC logotype

Diff of /MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F

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

revision 1.1 by heimbach, Thu Mar 29 15:59:21 2012 UTC revision 1.6 by dgoldberg, Mon Jul 30 19:04:55 2012 UTC
# Line 21  C     === Global variables === Line 21  C     === Global variables ===
21  #include "PARAMS.h"  #include "PARAMS.h"
22  #include "STREAMICE.h"  #include "STREAMICE.h"
23  #include "STREAMICE_CG.h"  #include "STREAMICE_CG.h"
24    #ifdef ALLOW_AUTODIFF_TAMC
25    # include "tamc.h"
26    #endif
27    
28  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
29        INTEGER myThid        INTEGER myThid
# Line 39  C     LOCAL VARIABLES Line 42  C     LOCAL VARIABLES
42    
43        INTEGER conv_flag, i, j, k, l, iter, cg_iters, bi, bj        INTEGER conv_flag, i, j, k, l, iter, cg_iters, bi, bj
44        INTEGER iter_numconv        INTEGER iter_numconv
45        _RL err_max, err_tempu, err_tempv, err_init, area        INTEGER ikey_nl
46          _RL err_max, err_tempu, err_tempv, err_init, area, err_max_fp
47        _RL max_vel, tempu, tempv, err_lastchange, cgtol        _RL max_vel, tempu, tempv, err_lastchange, cgtol
48        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
49  !       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+Oly,nSx,nSy)  !       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+Oly,nSx,nSy)
# Line 67  C     LOCAL VARIABLES Line 71  C     LOCAL VARIABLES
71         ENDDO         ENDDO
72        ENDDO        ENDDO
73    
74    #ifdef ALLOW_AUTODIFF_TAMC
75    !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
76    !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
77    #endif
78    
79        CALL STREAMICE_VISC_BETA ( myThid )        CALL STREAMICE_VISC_BETA ( myThid )
80    
81        _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )        _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
# Line 85  C     LOCAL VARIABLES Line 94  C     LOCAL VARIABLES
94         ENDDO         ENDDO
95        ENDDO        ENDDO
96    
97    
98        CALL STREAMICE_CG_BOUND_VALS( myThid,            CALL STREAMICE_CG_BOUND_VALS( myThid,    
99       O    ubd_SI,       O    ubd_SI,
100       O    vbd_SI)       O    vbd_SI)
101    
102    #ifdef ALLOW_AUTODIFF_TAMC
103    !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
104    !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
105    #endif
106    
107        CALL STREAMICE_CG_ACTION( myThid,        CALL STREAMICE_CG_ACTION( myThid,
108       O    Au_SI,       O    Au_SI,
109       O    Av_SI,       O    Av_SI,
# Line 97  C     LOCAL VARIABLES Line 112  C     LOCAL VARIABLES
112       I    0, sNx+1, 0, sNy+1 )       I    0, sNx+1, 0, sNy+1 )
113    
114        err_init = 0. _d 0        err_init = 0. _d 0
       err_tempu = 0. _d 0  
       err_tempv = 0. _d 0  
115    
116        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
117         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
118    #ifdef ALLOW_AUTODIFF_TAMC
119    !$TAF STORE err_init = comlev1, key=ikey_dynamics
120    #endif
121          DO j=1,sNy          DO j=1,sNy
122           DO i=1,sNx           DO i=1,sNx
123              err_tempu = 0. _d 0
124              err_tempv = 0. _d 0
125            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
126             err_tempu =             err_tempu =
127       &      ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -       &      ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -
# Line 119  C     LOCAL VARIABLES Line 137  C     LOCAL VARIABLES
137          ENDDO          ENDDO
138         ENDDO         ENDDO
139        ENDDO        ENDDO
140    #ifdef ALLOW_AUTODIFF_TAMC
141    !$TAF STORE err_init = comlev1, key=ikey_dynamics
142    #endif
143    
144        CALL GLOBAL_MAX_R8 (err_init, myThid)        CALL GLOBAL_MAX_R8 (err_init, myThid)
145    
146          WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',
147         &                       err_init
148          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
149         &                    SQUEEZE_RIGHT , 1)
150    
151    
152        iter_numconv = 0        iter_numconv = 0
153        err_max = err_init        err_max = err_init
154          err_max_fp = streamice_nonlin_tol_fp * 10.
155        err_lastchange = err_init        err_lastchange = err_init
156    
157    C START NL ITER. LOOP
158    C-------------------------------------------------------------------
159    
160        DO iter=1,streamice_max_nl_iter        DO iter=1,streamice_max_nl_iter
161    
162  C     To avoid using "exit", loop goes through all iterations  C     To avoid using "exit", loop goes through all iterations
163  C       but after convergence loop does nothing  C       but after convergence loop does nothing
164    
165         IF (err_max .GT. streamice_nonlin_tol * err_init) THEN  #ifdef ALLOW_AUTODIFF_TAMC
166             ikey_nl = (ikey_dynamics-1)*streamice_max_nl + iter
167    #endif
168    #ifdef ALLOW_AUTODIFF_TAMC
169    !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
170    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
171    !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
172    !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
173    !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
174    !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
175    !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
176    !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
177    !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
178    !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
179    !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
180    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
181    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
182    !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
183    !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
184    !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
185    !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
186    #endif
187    
188           IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
189         &     (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
190    
191         iter_numconv = iter_numconv + 1         iter_numconv = iter_numconv + 1
192    
193         CALL STREAMICE_CG_SOLVE(  #ifdef ALLOW_AUTODIFF_TAMC
194       &  U_streamice,  !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
195       &  V_streamice,  !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
196       &  taudx_SI,  !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
197       &  taudy_SI,  !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
198       &  myThid , cgtol , cg_iters)  !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
199    !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
200    !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
201    !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
202    !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
203    !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
204    !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
205    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
206    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
207    !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
208    !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
209    !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
210    !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
211    #endif
212    
213    #ifdef ALLOW_AUTODIFF_TAMC
214    !       DO bj = myByLo(myThid), myByHi(myThid)
215    !        DO bi = myBxLo(myThid), myBxHi(myThid)
216    !         DO j=1-OLy,sNy+OLy
217    !          DO i=1-OLx,sNx+OLx
218    !           U_streamice (i,j,bi,bj) = 0. _d 0
219    !           V_streamice (i,j,bi,bj) = 0. _d 0
220    !          ENDDO
221    !         ENDDO
222    !        ENDDO
223    !       ENDDO
224    #endif
225    
226           CALL STREAMICE_CG_WRAPPER(
227         &       U_streamice,
228         &       V_streamice,
229         &       taudx_SI,
230         &       taudy_SI,
231         &       cgtol,
232         &       cg_iters,
233         &       myThid )
234    
235                
236    
# Line 151  C       but after convergence loop does Line 241  C       but after convergence loop does
241          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
242       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
243    
244    #ifdef ALLOW_AUTODIFF_TAMC
245    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
246    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
247    #endif
248    
249         CALL STREAMICE_VISC_BETA ( myThid )         CALL STREAMICE_VISC_BETA ( myThid )
250    
251    
252    #ifdef ALLOW_AUTODIFF_TAMC
253    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
254    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
255    #endif
256    
257         _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )         _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
258         _EXCH_XY_RL ( visc_streamice , myThid )         _EXCH_XY_RL ( visc_streamice , myThid )
259    
# Line 173  C       but after convergence loop does Line 274  C       but after convergence loop does
274       O    ubd_SI,       O    ubd_SI,
275       O    vbd_SI)       O    vbd_SI)
276    
277    #ifdef ALLOW_AUTODIFF_TAMC
278    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
279    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
280    #endif
281    
282         CALL STREAMICE_CG_ACTION( myThid,         CALL STREAMICE_CG_ACTION( myThid,
283       O    Au_SI,       O    Au_SI,
284       O    Av_SI,       O    Av_SI,
# Line 180  C       but after convergence loop does Line 286  C       but after convergence loop does
286       I    V_streamice,       I    V_streamice,
287       I    0, sNx+1, 0, sNy+1 )       I    0, sNx+1, 0, sNy+1 )
288    
289    #ifdef ALLOW_AUTODIFF_TAMC
290    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
291    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
292    #endif
293    
294           DO bj = myByLo(myThid), myByHi(myThid)
295            DO bi = myBxLo(myThid), myBxHi(myThid)
296             DO j=1-OLy,sNy+OLy
297              DO i=1-OLx,sNx+OLx
298               u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
299               v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
300              ENDDO
301             ENDDO
302            ENDDO
303           ENDDO
304    
305         err_max = 0. _d 0         err_max = 0. _d 0
        err_tempu = 0. _d 0  
        err_tempv = 0. _d 0  
306    
307    #ifdef ALLOW_AUTODIFF_TAMC
308    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
309    #endif
310         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
311          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
312    #ifdef ALLOW_AUTODIFF_TAMC
313    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
314    #endif
315           DO j=1,sNy           DO j=1,sNy
316            DO i=1,sNx            DO i=1,sNx
317               err_tempu = 0. _d 0
318               err_tempv = 0. _d 0
319             IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN             IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
320              err_tempu =              err_tempu =
321       &       ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -       &       ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -
# Line 206  C       but after convergence loop does Line 334  C       but after convergence loop does
334    
335         CALL GLOBAL_MAX_R8 (err_max, myThid)         CALL GLOBAL_MAX_R8 (err_max, myThid)
336    
337           DO bj = myByLo(myThid), myByHi(myThid)
338            DO bi = myBxLo(myThid), myBxHi(myThid)
339             DO j=1,sNy
340              DO i=1,sNx
341               err_tempu = 0. _d 0
342               err_tempv = 0. _d 0
343               IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
344                err_tempu =
345         &       ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))
346               ENDIF
347               IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
348                err_tempv = MAX( err_tempu,
349         &       ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))
350               ENDIF
351               IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv
352              ENDDO
353             ENDDO
354            ENDDO
355           ENDDO
356    
357           CALL GLOBAL_MAX_R8 (err_max_fp, myThid)
358    
359         WRITE(msgBuf,'(A,F11.7)') 'err/err_init',         WRITE(msgBuf,'(A,F11.7)') 'err/err_init',
360       &                       err_max/err_init       &                       err_max/err_init
361         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
# Line 215  C       but after convergence loop does Line 365  C       but after convergence loop does
365       &   STREAMICE_lower_cg_tol) THEN       &   STREAMICE_lower_cg_tol) THEN
366           cgtol = cgtol * 5.e-2           cgtol = cgtol * 5.e-2
367           err_lastchange = err_max           err_lastchange = err_max
368           WRITE(msgBuf,'(A,F11.7)') 'new cg tol: ',           WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
369       &                       cgtol       &                       cgtol
370           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
371       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
372         ENDIF         ENDIF
373    
374    
375         ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)         ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
376        ENDDO        ENDDO
377    
378    C END NL ITER. LOOP
379    C-------------------------------------------------------------------
380    
381        if (iter_numconv .lt. streamice_max_nl_iter) then        if (iter_numconv .lt. streamice_max_nl_iter) then
382         PRINT *, "VELOCITY SOLVE CONVERGED, ", iter_numconv,         WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
383       & " iterations"       &         iter_numconv, ' iterations'
384           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
385         &                     SQUEEZE_RIGHT , 1)
386        else        else
387         PRINT *, "VELOCITY SOLVE DID NOT CONVERGE IN ",         WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE NOT CONVERGED IN ',
388       & iter, " iterations"       &         iter_numconv, ' iterations'
389           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
390         &                     SQUEEZE_RIGHT , 1)
391        endif        endif
392    
393    !        DO bj = myByLo(myThid), myByHi(myThid)
394    !         DO bi = myBxLo(myThid), myBxHi(myThid)
395    !          DO j=1,sNy
396    !           DO i=1,sNx
397    !            U_streamice (i,j,bi,bj) = 0. _d 0
398    !            V_streamice (i,j,bi,bj) = 0. _d 0
399    !           ENDDO
400    !          ENDDO
401    !         ENDDO
402    !        ENDDO
403    !
404    !        CALL STREAMICE_CG_WRAPPER(
405    !      &       U_streamice,
406    !      &       V_streamice,
407    !      &       taudx_SI,
408    !      &       taudy_SI,
409    !      &       cgtol,
410    !      &       cg_iters,
411    !      &       myThid )
412    
413    
414        _EXCH_XY_RL (U_streamice, myThid)        _EXCH_XY_RL (U_streamice, myThid)
415        _EXCH_XY_RL (V_streamice, myThid)        _EXCH_XY_RL (V_streamice, myThid)
416    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22