/[MITgcm]/MITgcm/pkg/streamice/streamice_vel_solve.F
ViewVC logotype

Diff of /MITgcm/pkg/streamice/streamice_vel_solve.F

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

revision 1.7 by dgoldberg, Thu Apr 24 12:01:50 2014 UTC revision 1.8 by dgoldberg, Fri Sep 5 14:25:11 2014 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8  CBOP  CBOP
9        SUBROUTINE STREAMICE_VEL_SOLVE( myThid )        SUBROUTINE STREAMICE_VEL_SOLVE( myThid, maxNLIter, maxCGiter,
10         &                                myiter )
11  C     /============================================================\  C     /============================================================\
12  C     | SUBROUTINE                                                 |  C     | SUBROUTINE                                                 |
13  C     | o                                                          |  C     | o                                                          |
# Line 31  C     === Global variables === Line 32  C     === Global variables ===
32    
33  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
34        INTEGER myThid        INTEGER myThid
35          INTEGER maxNLIter
36          INTEGER maxCGIter
37          INTEGER myIter
38    
39  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
40    
# Line 44  C     LOCAL VARIABLES Line 48  C     LOCAL VARIABLES
48  !                         isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub  !                         isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
49  !   real                     :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi, rhow  !   real                     :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi, rhow
50    
51        INTEGER conv_flag, i, j, k, l, iter, cg_iters, bi, bj        INTEGER i, j, k, l, iter, cg_iters, bi, bj
52        INTEGER iter_numconv        INTEGER iter_numconv
53        INTEGER ikey_nl, myThidTemp        INTEGER ikey_nl, myThidTemp
54        _RL err_max, err_tempu, err_tempv, err_init, area, err_max_fp        _RL err_max, err_tempu, err_tempv, err_init, area, err_max_fp
55        _RL max_vel, tempu, tempv, err_lastchange, cgtol        _RL max_vel, tempu, tempv, err_lastchange, cgtol
56        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
57          LOGICAL CONVERGED
58  #ifdef ALLOW_PETSC  #ifdef ALLOW_PETSC
59  !      myThidTemp = myThid  !      myThidTemp = myThid
60  !      call streamice_initialize_petsc (myThidTemp)  !      call streamice_initialize_petsc (myThidTemp)
# Line 76  C     LOCAL VARIABLES Line 81  C     LOCAL VARIABLES
81        CALL STREAMICE_FORCED_BUTTRESS (myThid)        CALL STREAMICE_FORCED_BUTTRESS (myThid)
82  #endif  #endif
83    
84                CALL TIMER_START ('STREAMICE_VEL_SOLVE',myThid)
85    
86    
87        cgtol = streamice_cg_tol        cgtol = streamice_cg_tol
88          CONVERGED = .false.
89          err_max = 0.
90          err_max_fp = 0.
91    
 !       CALL WRITE_FULLARRAY_RL ("taudy_SI",taudy_SI,1,0,0,1,0,myThid)  
92    
93        _EXCH_XY_RL( taudx_SI , myThid )        _EXCH_XY_RL( taudx_SI , myThid )
94        _EXCH_XY_RL( taudy_SI , myThid )        _EXCH_XY_RL( taudy_SI , myThid )
95    
 !       CALL WRITE_FULLARRAY_RL ("taudy_SI_2",taudy_SI,1,0,0,1,0,myThid)  
96    
97        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
98         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 121  C     LOCAL VARIABLES Line 128  C     LOCAL VARIABLES
128        _EXCH_XY_RL( tau_beta_eff_streamice , myThid )        _EXCH_XY_RL( tau_beta_eff_streamice , myThid )
129        _EXCH_XY_RL( visc_streamice , myThid )        _EXCH_XY_RL( visc_streamice , myThid )
130    
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO j=1,sNy  
          DO i=1,sNx  
           Au_SI (i,j,bi,bj) = 0. _d 0  
           Av_SI (i,j,bi,bj) = 0. _d 0  
           ubd_SI (i,j,bi,bj) = 0. _d 0  
           vbd_SI (i,j,bi,bj) = 0. _d 0  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
   
131    
132        CALL STREAMICE_CG_BOUND_VALS( myThid,        if (STREAMICE_chkresidconvergence .or.
133       O    ubd_SI,       &    (streamice_maxnliter_cpl.eq.0 .and. myIter.eq.0)) then
      O    vbd_SI)  
134    
135            CALL STREAMICE_GET_VEL_RESID_ERR ( err_init, myThid )
 !      CALL WRITE_FLD_XY_RL("u_bound_cont","",ubd_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("v_bound_cont","",vbd_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("taudx_u","",taudx_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("taudx_v","",taudy_SI,0,myThid)  
136    
137  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 !$TAF STORE U_streamice = comlev1, key=ikey_dynamics  
 !$TAF STORE V_streamice = comlev1, key=ikey_dynamics  
 #endif  
   
       CALL STREAMICE_CG_ACTION( myThid,  
      O    Au_SI,  
      O    Av_SI,  
      I    U_streamice,  
      I    V_streamice,  
      I    0, sNx+1, 0, sNy+1 )  
   
   
   
       err_init = 0. _d 0  
   
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
138  !$TAF STORE err_init = comlev1, key=ikey_dynamics  !$TAF STORE err_init = comlev1, key=ikey_dynamics
139  #endif  #endif
         DO j=1,sNy  
          DO i=1,sNx  
           err_tempu = 0. _d 0  
           err_tempv = 0. _d 0  
           IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN  
            err_tempu =  
      &      ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -  
      &           taudx_SI(i,j,bi,bj))  
 !            print *, "err_temp_u", err_tempu  
           ENDIF  
           IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN  
            err_tempv = MAX( err_tempu,  
      &      ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -  
      &           taudy_SI(i,j,bi,bj)))  
           ENDIF  
           IF (err_tempv .ge. err_init) THEN  
             err_init = err_tempv  
           ENDIF  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
 #ifdef ALLOW_AUTODIFF_TAMC  
 !$TAF STORE err_init = comlev1, key=ikey_dynamics  
 #endif  
   
       CALL GLOBAL_MAX_R8 (err_init, myThid)  
140    
141        WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',         WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',
142       &                       err_init       &                       err_init
143        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
144       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
145    
146          endif !STREAMICE_chkresidconvergence
147    
148    
149        iter_numconv = 0        iter_numconv = 0
150        err_max = err_init        err_max = err_init
# Line 208  C     LOCAL VARIABLES Line 154  C     LOCAL VARIABLES
154  C START NL ITER. LOOP  C START NL ITER. LOOP
155  C-------------------------------------------------------------------  C-------------------------------------------------------------------
156    
157        DO iter=1,streamice_max_nl_iter        DO iter=1,maxNLIter
158    
159  C     To avoid using "exit", loop goes through all iterations  C     To avoid using "exit", loop goes through all iterations
160  C       but after convergence loop does nothing  C       but after convergence loop does nothing
# Line 239  C       but after convergence loop does Line 185  C       but after convergence loop does
185  !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl  !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
186  #endif  #endif
187    
188         IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.  !        IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
189       &     (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN  !      &     (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
190    
191    #ifdef ALLOW_AUTODIFF_TAMC
192    !$TAF STORE CONVERGED = comlev1_stream_nl, key=ikey_nl
193    #endif
194          
195           IF (.not.CONVERGED) THEN
196    
197    
198         iter_numconv = iter_numconv + 1         iter_numconv = iter_numconv + 1
199    
# Line 264  C       but after convergence loop does Line 217  C       but after convergence loop does
217  !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl  !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
218  #endif  #endif
219    
 #ifdef ALLOW_AUTODIFF_TAMC  
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1-OLy,sNy+OLy  
 !          DO i=1-OLx,sNx+OLx  
 !           U_streamice (i,j,bi,bj) = 0. _d 0  
 !           V_streamice (i,j,bi,bj) = 0. _d 0  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
 #endif  
220    
221         CALL STREAMICE_CG_WRAPPER(         CALL STREAMICE_CG_WRAPPER(
222       &       U_streamice,       &       U_streamice,
# Line 284  C       but after convergence loop does Line 225  C       but after convergence loop does
225       &       taudy_SI,       &       taudy_SI,
226       &       cgtol,       &       cgtol,
227       &       cg_iters,       &       cg_iters,
228         &       maxCGIter,
229       &       myThid )       &       myThid )
230    
231  #ifdef STREAMICE_HYBRID_STRESS  #ifdef STREAMICE_HYBRID_STRESS
# Line 297  C       but after convergence loop does Line 239  C       but after convergence loop does
239          CALL STREAMICE_TAUB (myThid)          CALL STREAMICE_TAUB (myThid)
240  #endif  #endif
241    
242    !-----------------------------------------------------------------------------
243         WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',         WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
244       &                       iter, " ",       &                       iter, " ",
245       &                       cg_iters,       &                       cg_iters,
246       &                       ' iterations '       &                       ' iterations '
247          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
248       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
249    !-----------------------------------------------------------------------------
250    
251  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
252  !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl  !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
# Line 321  C       but after convergence loop does Line 265  C       but after convergence loop does
265  #endif  #endif
266    
267    
268  #ifdef ALLOW_AUTODIFF_TAMC        _EXCH_XY_RL( tau_beta_eff_streamice , myThid )
269  !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl        _EXCH_XY_RL( visc_streamice , myThid )
 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl  
 #endif  
270    
271         _EXCH_XY_RL( tau_beta_eff_streamice , myThid )  !----------------------CONVERGENCE TESTS-------------------------------
        _EXCH_XY_RL( visc_streamice , myThid )  
272    
273         DO bj = myByLo(myThid), myByHi(myThid)        if (STREAMICE_chkresidconvergence .or.
274          DO bi = myBxLo(myThid), myBxHi(myThid)       &    (streamice_maxnliter_cpl.eq.0 .and. myIter.eq.0)) then
          DO j=1,sNy  
           DO i=1,sNx  
            Au_SI (i,j,bi,bj) = 0. _d 0  
            Av_SI (i,j,bi,bj) = 0. _d 0  
            ubd_SI (i,j,bi,bj) = 0. _d 0  
            vbd_SI (i,j,bi,bj) = 0. _d 0  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
275    
276         CALL STREAMICE_CG_BOUND_VALS( myThid,          CALL STREAMICE_GET_VEL_RESID_ERR ( err_max, myThid )
      O    ubd_SI,  
      O    vbd_SI)  
277    
278  #ifdef ALLOW_AUTODIFF_TAMC         WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
279  !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl       &                       err_max/err_init
280  !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
281  #endif       &                    SQUEEZE_RIGHT , 1)
282    
283         CALL STREAMICE_CG_ACTION( myThid,         IF (err_max .LE. streamice_nonlin_tol * err_init) THEN
284       O    Au_SI,          CONVERGED = .true.
285       O    Av_SI,         ENDIF
      I    U_streamice,  
      I    V_streamice,  
      I    0, sNx+1, 0, sNy+1 )  
   
 !      if (iter .eq. streamice_max_nl_iter) then  
 !      CALL WRITE_FLD_XY_RL("u_bound_cont_A","",ubd_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("v_bound_cont_A","",vbd_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("u_bound_cont_B","",Au_SI,0,myThid)  
 !      CALL WRITE_FLD_XY_RL("v_bound_cont_B","",Av_SI,0,myThid)  
 !      endif  
286    
 #ifdef ALLOW_AUTODIFF_TAMC  
 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl  
 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl  
 #endif  
287    
288         err_max = 0. _d 0        endif
        err_max_fp = 0. _d 0  
289    
 #ifdef ALLOW_AUTODIFF_TAMC  
 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl  
 #endif  
        DO bj = myByLo(myThid), myByHi(myThid)  
         DO bi = myBxLo(myThid), myBxHi(myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl  
 #endif  
          DO j=1,sNy  
           DO i=1,sNx  
            err_tempu = 0. _d 0  
            err_tempv = 0. _d 0  
            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN  
             err_tempu =  
      &       ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -  
      &            taudx_SI(i,j,bi,bj))  
            ENDIF  
            IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN  
             err_tempv = MAX( err_tempu,  
      &       ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -  
      &            taudy_SI(i,j,bi,bj)))  
            ENDIF  
 !           if (err_tempu.ge.1.e2.or.err_tempv.ge.1.e2) THEN  
 !            print *, "FOUND MAX ", i,j,err_tempu,err_tempv,  
 !     &      ubd_SI(i,j,bi,bj),vbd_SI(i,j,bi,bj)  
 !           endif  
            IF (err_tempv .ge. err_max) THEN  
             err_max = err_tempv  
            ENDIF  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
290    
        CALL GLOBAL_MAX_R8 (err_max, myThid)  
291    
292         DO bj = myByLo(myThid), myByHi(myThid)        if (STREAMICE_chkfixedptconvergence .or.
293          DO bi = myBxLo(myThid), myBxHi(myThid)       &    (streamice_maxnliter_cpl.eq.0 .and. myIter.eq.0)) then
          DO j=1,sNy  
           DO i=1,sNx  
            err_tempu = 0. _d 0  
            err_tempv = 0. _d 0  
            IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN  
             err_tempu =  
      &       ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))  
            ENDIF  
            IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN  
             err_tempv = MAX( err_tempu,  
      &       ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))  
            ENDIF  
            IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
294    
295         CALL GLOBAL_MAX_R8 (err_max_fp, myThid)          CALL STREAMICE_GET_VEL_FP_ERR ( err_max_fp, myThid )
        WRITE(msgBuf,'(A,1PE22.14)') 'STREAMICE_FP_ERROR =',  
      &                       err_max_fp  
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,  
      &                    SQUEEZE_RIGHT , 1)  
296    
297         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
298          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 446  C       but after convergence loop does Line 305  C       but after convergence loop does
305          ENDDO          ENDDO
306         ENDDO         ENDDO
307    
308         WRITE(msgBuf,'(A,E15.7)') 'err/err_init',         WRITE(msgBuf,'(A,1PE22.14)') 'STREAMICE_FP_ERROR =',
309       &                       err_max/err_init       &                       err_max_fp
310         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
311       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
312    
313           IF (err_max_fp .LE. streamice_nonlin_tol_fp) THEN
314            CONVERGED = .true.
315           ENDIF
316    
317          endif
318    
319    
320    
321    
322    !----------------------END CONVERGENCE TESTS-------------------------------
323    
324    
325    
326    
327         IF (err_max<err_lastchange*1.e-2 .and.         IF (err_max<err_lastchange*1.e-2 .and.
328       &   STREAMICE_lower_cg_tol) THEN       &   STREAMICE_lower_cg_tol) THEN
329           cgtol = cgtol * 5.e-2           cgtol = cgtol * 5.e-2
# Line 487  C--------------------------------------- Line 360  C---------------------------------------
360    
361        _EXCH_XY_RL(U_streamice, myThid)        _EXCH_XY_RL(U_streamice, myThid)
362        _EXCH_XY_RL(V_streamice, myThid)        _EXCH_XY_RL(V_streamice, myThid)
363  !      CALL WRITE_FLD_XY_RL("taubx","",streamice_taubx,0,myThid)  
364  !      CALL WRITE_FLD_XY_RL("tauby","",streamice_tauby,0,myThid)  
365          CALL TIMER_STOP ('STREAMICE_VEL_SOLVE',myThid)
366    
367  #endif  #endif
368        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22