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