/[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.9 by heimbach, Thu Sep 20 02:04:45 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        CALL STREAMICE_VISC_BETA ( myThid )  #ifdef ALLOW_AUTODIFF_TAMC
75    !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
76    !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
77    #ifdef STREAMICE_HYBRID_STRESS
78    !$TAF STORE streamice_taubx = comlev1, key=ikey_dynamics
79    !$TAF STORE streamice_tauby = comlev1, key=ikey_dynamics
80    !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
81    #endif
82    #endif
83    
84    #ifdef STREAMICE_HYBRID_STRESS
85           CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
86    #else
87           CALL STREAMICE_VISC_BETA ( myThid )
88    #endif
89    
90    #ifdef STREAMICE_HYBRID_STRESS
91    !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
92    #endif
93    
94        _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )        _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
95        _EXCH_XY_RL ( visc_streamice , myThid )        _EXCH_XY_RL ( visc_streamice , myThid )
# Line 85  C     LOCAL VARIABLES Line 107  C     LOCAL VARIABLES
107         ENDDO         ENDDO
108        ENDDO        ENDDO
109    
110    
111        CALL STREAMICE_CG_BOUND_VALS( myThid,            CALL STREAMICE_CG_BOUND_VALS( myThid,    
112       O    ubd_SI,       O    ubd_SI,
113       O    vbd_SI)       O    vbd_SI)
114    
115    #ifdef ALLOW_AUTODIFF_TAMC
116    !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
117    !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
118    #endif
119    
120        CALL STREAMICE_CG_ACTION( myThid,        CALL STREAMICE_CG_ACTION( myThid,
121       O    Au_SI,       O    Au_SI,
122       O    Av_SI,       O    Av_SI,
# Line 97  C     LOCAL VARIABLES Line 125  C     LOCAL VARIABLES
125       I    0, sNx+1, 0, sNy+1 )       I    0, sNx+1, 0, sNy+1 )
126    
127        err_init = 0. _d 0        err_init = 0. _d 0
       err_tempu = 0. _d 0  
       err_tempv = 0. _d 0  
128    
129        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
130         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
131    #ifdef ALLOW_AUTODIFF_TAMC
132    !$TAF STORE err_init = comlev1, key=ikey_dynamics
133    #endif
134          DO j=1,sNy          DO j=1,sNy
135           DO i=1,sNx           DO i=1,sNx
136              err_tempu = 0. _d 0
137              err_tempv = 0. _d 0
138            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
139             err_tempu =             err_tempu =
140       &      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) -
141       &           taudx_SI(i,j,bi,bj))       &           taudx_SI(i,j,bi,bj))
142    !            print *, "err_temp_u", err_tempu
143            ENDIF            ENDIF
144            IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN            IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
145             err_tempv = MAX( err_tempu,             err_tempv = MAX( err_tempu,
# Line 119  C     LOCAL VARIABLES Line 151  C     LOCAL VARIABLES
151          ENDDO          ENDDO
152         ENDDO         ENDDO
153        ENDDO        ENDDO
154    #ifdef ALLOW_AUTODIFF_TAMC
155    !$TAF STORE err_init = comlev1, key=ikey_dynamics
156    #endif
157    
158        CALL GLOBAL_MAX_R8 (err_init, myThid)        CALL GLOBAL_MAX_R8 (err_init, myThid)
159    
160          WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',
161         &                       err_init
162          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
163         &                    SQUEEZE_RIGHT , 1)
164    
165    
166        iter_numconv = 0        iter_numconv = 0
167        err_max = err_init        err_max = err_init
168          err_max_fp = streamice_nonlin_tol_fp * 10.
169        err_lastchange = err_init        err_lastchange = err_init
170    
171    C START NL ITER. LOOP
172    C-------------------------------------------------------------------
173    
174        DO iter=1,streamice_max_nl_iter        DO iter=1,streamice_max_nl_iter
175    
176  C     To avoid using "exit", loop goes through all iterations  C     To avoid using "exit", loop goes through all iterations
177  C       but after convergence loop does nothing  C       but after convergence loop does nothing
178    
179         IF (err_max .GT. streamice_nonlin_tol * err_init) THEN  #ifdef ALLOW_AUTODIFF_TAMC
180             ikey_nl = (ikey_dynamics-1)*streamice_max_nl + iter
181    #endif
182    #ifdef ALLOW_AUTODIFF_TAMC
183    !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
184    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
185    !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
186    !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
187    !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
188    !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
189    !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
190    !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
191    !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
192    !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
193    !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
194    !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
195    !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
196    !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
197    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
198    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
199    !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
200    !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
201    !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
202    !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
203    #endif
204    
205         iter_numconv = iter_numconv + 1         IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
206         &     (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
207    
208         CALL STREAMICE_CG_SOLVE(         iter_numconv = iter_numconv + 1
      &  U_streamice,  
      &  V_streamice,  
      &  taudx_SI,  
      &  taudy_SI,  
      &  myThid , cgtol , cg_iters)  
209    
210          #ifdef ALLOW_AUTODIFF_TAMC
211    !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
212    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
213    !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
214    !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
215    !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
216    !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
217    !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
218    !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
219    !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
220    !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
221    !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
222    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
223    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
224    !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
225    !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
226    !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
227    !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
228    #endif
229    
230    #ifdef ALLOW_AUTODIFF_TAMC
231    !       DO bj = myByLo(myThid), myByHi(myThid)
232    !        DO bi = myBxLo(myThid), myBxHi(myThid)
233    !         DO j=1-OLy,sNy+OLy
234    !          DO i=1-OLx,sNx+OLx
235    !           U_streamice (i,j,bi,bj) = 0. _d 0
236    !           V_streamice (i,j,bi,bj) = 0. _d 0
237    !          ENDDO
238    !         ENDDO
239    !        ENDDO
240    !       ENDDO
241    #endif
242    
243           CALL STREAMICE_CG_WRAPPER(
244         &       U_streamice,
245         &       V_streamice,
246         &       taudx_SI,
247         &       taudy_SI,
248         &       cgtol,
249         &       cg_iters,
250         &       myThid )
251    
252    #ifdef STREAMICE_HYBRID_STRESS
253    #ifdef ALLOW_AUTODIFF_TAMC
254    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
255    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
256    #endif
257    #endif
258    
259    #ifdef STREAMICE_HYBRID_STRESS
260            CALL STREAMICE_TAUB (myThid)
261    #endif
262    
263         WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',         WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
264       &                       iter, " ",       &                       iter, " ",
# Line 151  C       but after convergence loop does Line 267  C       but after convergence loop does
267          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
268       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
269    
270    #ifdef ALLOW_AUTODIFF_TAMC
271    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
272    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
273    #ifdef STREAMICE_HYBRID_STRESS
274    !$TAF STORE streamice_taubx = comlev1_stream_nl, key=ikey_nl
275    !$TAF STORE streamice_tauby = comlev1_stream_nl, key=ikey_nl
276    !$TAF STORE visc_streamice_full = comlev1_stream_nl, key=ikey_nl
277    #endif
278    #endif
279    
280    #ifdef STREAMICE_HYBRID_STRESS
281           CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
282    #else
283         CALL STREAMICE_VISC_BETA ( myThid )         CALL STREAMICE_VISC_BETA ( myThid )
284    #endif
285    
286    
287    #ifdef ALLOW_AUTODIFF_TAMC
288    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
289    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
290    #endif
291    
292         _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )         _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
293         _EXCH_XY_RL ( visc_streamice , myThid )         _EXCH_XY_RL ( visc_streamice , myThid )
294    
# Line 173  C       but after convergence loop does Line 309  C       but after convergence loop does
309       O    ubd_SI,       O    ubd_SI,
310       O    vbd_SI)       O    vbd_SI)
311    
312    #ifdef ALLOW_AUTODIFF_TAMC
313    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
314    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
315    #endif
316    
317         CALL STREAMICE_CG_ACTION( myThid,         CALL STREAMICE_CG_ACTION( myThid,
318       O    Au_SI,       O    Au_SI,
319       O    Av_SI,       O    Av_SI,
# Line 180  C       but after convergence loop does Line 321  C       but after convergence loop does
321       I    V_streamice,       I    V_streamice,
322       I    0, sNx+1, 0, sNy+1 )       I    0, sNx+1, 0, sNy+1 )
323    
324    #ifdef ALLOW_AUTODIFF_TAMC
325    !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
326    !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
327    #endif
328    
329         err_max = 0. _d 0         err_max = 0. _d 0
330         err_tempu = 0. _d 0         err_max_fp = 0. _d 0
        err_tempv = 0. _d 0  
331    
332    #ifdef ALLOW_AUTODIFF_TAMC
333    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
334    #endif
335         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
336          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
337    #ifdef ALLOW_AUTODIFF_TAMC
338    !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
339    #endif
340           DO j=1,sNy           DO j=1,sNy
341            DO i=1,sNx            DO i=1,sNx
342               err_tempu = 0. _d 0
343               err_tempv = 0. _d 0
344             IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN             IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
345              err_tempu =              err_tempu =
346       &       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 359  C       but after convergence loop does
359    
360         CALL GLOBAL_MAX_R8 (err_max, myThid)         CALL GLOBAL_MAX_R8 (err_max, myThid)
361    
362         WRITE(msgBuf,'(A,F11.7)') 'err/err_init',         DO bj = myByLo(myThid), myByHi(myThid)
363            DO bi = myBxLo(myThid), myBxHi(myThid)
364             DO j=1,sNy
365              DO i=1,sNx
366               err_tempu = 0. _d 0
367               err_tempv = 0. _d 0
368               IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
369                err_tempu =
370         &       ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))
371               ENDIF
372               IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
373                err_tempv = MAX( err_tempu,
374         &       ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))
375               ENDIF
376               IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv
377              ENDDO
378             ENDDO
379            ENDDO
380           ENDDO
381    
382           CALL GLOBAL_MAX_R8 (err_max_fp, myThid)
383           WRITE(msgBuf,'(A,E15.7)') '||x_i-x_{i-1}||_inf',
384         &                       err_max_fp
385           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
386         &                    SQUEEZE_RIGHT , 1)
387    
388           DO bj = myByLo(myThid), myByHi(myThid)
389            DO bi = myBxLo(myThid), myBxHi(myThid)
390             DO j=1-OLy,sNy+OLy
391              DO i=1-OLx,sNx+OLx
392               u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
393               v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
394              ENDDO
395             ENDDO
396            ENDDO
397           ENDDO
398    
399           WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
400       &                       err_max/err_init       &                       err_max/err_init
401         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
402       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
# Line 215  C       but after convergence loop does Line 405  C       but after convergence loop does
405       &   STREAMICE_lower_cg_tol) THEN       &   STREAMICE_lower_cg_tol) THEN
406           cgtol = cgtol * 5.e-2           cgtol = cgtol * 5.e-2
407           err_lastchange = err_max           err_lastchange = err_max
408           WRITE(msgBuf,'(A,F11.7)') 'new cg tol: ',           WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
409       &                       cgtol       &                       cgtol
410           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
411       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
412         ENDIF         ENDIF
413    
414    
415         ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)         ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
416        ENDDO        ENDDO
417    
418    C END NL ITER. LOOP
419    C-------------------------------------------------------------------
420    
421        if (iter_numconv .lt. streamice_max_nl_iter) then        if (iter_numconv .lt. streamice_max_nl_iter) then
422         PRINT *, "VELOCITY SOLVE CONVERGED, ", iter_numconv,         WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
423       & " iterations"       &         iter_numconv, ' iterations'
424           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
425         &                     SQUEEZE_RIGHT , 1)
426        else        else
427         PRINT *, "VELOCITY SOLVE DID NOT CONVERGE IN ",         WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE NOT CONVERGED IN ',
428       & iter, " iterations"       &         iter_numconv, ' iterations'
429           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
430         &                     SQUEEZE_RIGHT , 1)
431        endif        endif
432    
433        _EXCH_XY_RL (U_streamice, myThid)        _EXCH_XY_RL (U_streamice, myThid)
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22