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

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

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

revision 1.3 by dgoldberg, Mon May 14 16:52:29 2012 UTC revision 1.4 by dgoldberg, Thu Jul 19 18:46:56 2012 UTC
# Line 7  C---+----1----+----2----+----3----+----4 Line 7  C---+----1----+----2----+----3----+----4
7    
8  CBOP  CBOP
9        SUBROUTINE STREAMICE_CG_SOLVE(        SUBROUTINE STREAMICE_CG_SOLVE(
10       U                               cg_Uin,       U                               cg_Uin,     ! x-velocities
11       U                               cg_Vin,       U                               cg_Vin,     ! y-velocities
12       I                               cg_Bu,       I                               cg_Bu,      ! force in x dir
13       I                               cg_Bv,       I                               cg_Bv,      ! force in y dir
14         I                               A_uu,       ! section of matrix that multiplies u and projects on u
15         I                               A_uv,       ! section of matrix that multiplies v and projects on u
16         I                               A_vu,       ! section of matrix that multiplies u and projects on v
17         I                               A_vv,       ! section of matrix that multiplies v and projects on v
18       I                               tolerance,       I                               tolerance,
19       O                               iters,       O                               iters,
20       I                               myThid )       I                               myThid )
# Line 40  C     cg_Bu, cg_Bv - driving stress Line 44  C     cg_Bu, cg_Bv - driving stress
44        _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45        _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46        _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47          _RL
48         & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
49         & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
50         & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
51         & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
52    
53  C     LOCAL VARIABLES  C     LOCAL VARIABLES
54        INTEGER i, j, bi, bj, cg_halo, conv_flag        INTEGER i, j, bi, bj, cg_halo, conv_flag
55        INTEGER iter, is, js, ie, je, colx, coly, k        INTEGER iter, is, js, ie, je, colx, coly, k
56        _RL dot_p1, dot_p2, resid_0, alpha_k, beta_k, resid        _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
57        _RL dot_p1_tile (nSx,nSy)        _RL dot_p1_tile (nSx,nSy)
58        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
59    
# Line 52  C     LOCAL VARIABLES Line 61  C     LOCAL VARIABLES
61    
62  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
63    
64    
65        conv_flag = 0        conv_flag = 0
66    
67    
68        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
69         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
70          DO j=1,sNy          DO j=1,sNy
71           DO i=1,sNx           DO i=1,sNx
72            Zu_SI (i,j,bi,bj) = 0. _d 0            Zu_SI (i,j,bi,bj) = 0. _d 0
73            Zv_SI (i,j,bi,bj) = 0. _d 0            Zv_SI (i,j,bi,bj) = 0. _d 0
           DIAGu_SI (i,j,bi,bj) = 0. _d 0  
           DIAGv_SI (i,j,bi,bj) = 0. _d 0  
74            Ru_SI (i,j,bi,bj) = 0. _d 0            Ru_SI (i,j,bi,bj) = 0. _d 0
75            Rv_SI (i,j,bi,bj) = 0. _d 0            Rv_SI (i,j,bi,bj) = 0. _d 0
76            Au_SI (i,j,bi,bj) = 0. _d 0            Au_SI (i,j,bi,bj) = 0. _d 0
77            Av_SI (i,j,bi,bj) = 0. _d 0            Av_SI (i,j,bi,bj) = 0. _d 0
78            Du_SI (i,j,bi,bj) = 0. _d 0            Du_SI (i,j,bi,bj) = 0. _d 0
79            Dv_SI (i,j,bi,bj) = 0. _d 0            Dv_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  
80           ENDDO           ENDDO
81          ENDDO          ENDDO
82         ENDDO         ENDDO
83        ENDDO        ENDDO
84          
85    C     FIND INITIAL RESIDUAL, and initialize r
86    
87  C     PREAMBLE: get bdry contribution, find matrix diagonal,  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
 C     initialize iterates R (residual), Z, and D  
88    
89        CALL STREAMICE_CG_BOUND_VALS( myThid,              
      O    ubd_SI,  
      O    vbd_SI)  
90    
91        DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
92         DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
93          DO j=1-OLy,sNy+OLy            DO j=1,sNy
94           DO i=1-OLx,sNx+OLx             DO i=1,sNx
95            RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)              DO colx=-1,1
96       &     - ubd_SI(i,j,bi,bj)               DO coly=-1,1
97            RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
98       &     - vbd_SI(i,j,bi,bj)       &         A_uu(i,j,bi,bj,colx,coly)*
99           ENDDO       &         cg_Uin(i+colx,j+coly,bi,bj)+
100          ENDDO       &         A_uv(i,j,bi,bj,colx,coly)*    
101         ENDDO       &         cg_Vin(i+colx,j+coly,bi,bj)
       ENDDO  
         
       _EXCH_XY_RL( RHSu_SI, myThid )  
       _EXCH_XY_RL( RHSv_SI, myThid )  
102    
       CALL STREAMICE_CG_ADIAG( myThid,  
      O    DIAGu_SI,  
      O    DIAGv_SI)  
103    
104        _EXCH_XY_RL( DIAGu_SI, myThid )                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
105        _EXCH_XY_RL( DIAGv_SI, myThid )       &         A_vu(i,j,bi,bj,colx,coly)*
106               &         cg_Uin(i+colx,j+coly,bi,bj)+
107  C     FIX PROBLEM WITH PRECOND LATER       &         A_vv(i,j,bi,bj,colx,coly)*    
108         &         cg_Vin(i+colx,j+coly,bi,bj)
109  !       DO bj = myByLo(myThid), myByHi(myThid)               ENDDO
110  !        DO bi = myBxLo(myThid), myBxHi(myThid)              ENDDO
111  !         DO j=1-OLy,sNy+OLy             ENDDO
112  !          DO i=1-OLx,sNx+OLx            ENDDO
113  !           DIAGu_SI(i,j,bi,bj)=1.0           ENDDO
114  !           DIAGv_SI(i,j,bi,bj)=1.0          ENDDO
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
         
       CALL STREAMICE_CG_ACTION( myThid,  
      O    Au_SI,  
      O    Av_SI,  
      I    cg_Uin,  
      I    cg_Vin,  
      I    0, sNx+1, 0, sNy+1 )  
115    
116                
117        _EXCH_XY_RL( Au_SI, myThid )        _EXCH_XY_RL( Au_SI, myThid )
118        _EXCH_XY_RL( Av_SI, myThid )        _EXCH_XY_RL( Av_SI, myThid )
119    
120    
121        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
122         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
123          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
124           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
125            Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-            Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
126       &     Au_SI(i,j,bi,bj)       &     Au_SI(i,j,bi,bj)
127            Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-            Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
128       &     Av_SI(i,j,bi,bj)       &     Av_SI(i,j,bi,bj)
129           ENDDO           ENDDO
130          ENDDO          ENDDO
# Line 143  C     FIX PROBLEM WITH PRECOND LATER Line 133  C     FIX PROBLEM WITH PRECOND LATER
133         ENDDO         ENDDO
134        ENDDO        ENDDO
135    
136          
137    
138        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
139         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
140          DO j=1,sNy          DO j=1,sNy
# Line 159  C     FIX PROBLEM WITH PRECOND LATER Line 151  C     FIX PROBLEM WITH PRECOND LATER
151        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
152        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
153    
154    C    CCCCCCCCCCCCCCCCCCCC
155    
156        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
157         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
158          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 172  C     FIX PROBLEM WITH PRECOND LATER Line 166  C     FIX PROBLEM WITH PRECOND LATER
166         ENDDO         ENDDO
167        ENDDO        ENDDO
168    
   
169        cg_halo = min(OLx-1,OLy-1)        cg_halo = min(OLx-1,OLy-1)
170        conv_flag = 0        conv_flag = 0
171    
# Line 195  c  !!              !! Line 188  c  !!              !!
188  c  !! MAIN CG LOOP !!  c  !! MAIN CG LOOP !!
189  c  !!              !!  c  !!              !!
190  c  !!!!!!!!!!!!!!!!!!  c  !!!!!!!!!!!!!!!!!!
191    
192    
193        
194        
195  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
196    
197        print *, "BEGINNING MAIN CG LOOP"        print *, "BEGINNING MAIN CG LOOP"
198    
199  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
 #ifdef STREAMICE_CONSTRUCT_MATRIX  
   
       CALL STREAMICE_CG_MAKE_A(myThid)  
200    
 #endif  
201    
202        do iter = 1, streamice_max_cg_iter        do iter = 1, streamice_max_cg_iter
203         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
# Line 224  c      to avoid using "exit" Line 216  c      to avoid using "exit"
216            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
217             Au_SI(i,j,bi,bj) = 0. _d 0             Au_SI(i,j,bi,bj) = 0. _d 0
218             Av_SI(i,j,bi,bj) = 0. _d 0             Av_SI(i,j,bi,bj) = 0. _d 0
 !            Du_SI(i,j,bi,bj) = Real(i)  
 !            Dv_SI(i,j,bi,bj) = 0.0  
219            ENDDO            ENDDO
220           ENDDO           ENDDO
221          ENDDO          ENDDO
# Line 233  c      to avoid using "exit" Line 223  c      to avoid using "exit"
223    
224  !        IF (STREAMICE_construct_matrix) THEN  !        IF (STREAMICE_construct_matrix) THEN
225    
226  #ifdef STREAMICE_CONSTRUCT_MATRIX  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
227    
228          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
229           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 242  c      to avoid using "exit" Line 232  c      to avoid using "exit"
232              DO colx=-1,1              DO colx=-1,1
233               DO coly=-1,1               DO coly=-1,1
234                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
235       &         streamice_cg_A1(i,j,bi,bj,colx,coly)*       &         A_uu(i,j,bi,bj,colx,coly)*
236       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
237       &         streamice_cg_A2(i,j,bi,bj,colx,coly)*           &         A_uv(i,j,bi,bj,colx,coly)*    
238       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
239                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
240       &         streamice_cg_A3(i,j,bi,bj,colx,coly)*       &         A_vu(i,j,bi,bj,colx,coly)*
241       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
242       &         streamice_cg_A4(i,j,bi,bj,colx,coly)*           &         A_vv(i,j,bi,bj,colx,coly)*    
243       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
244               ENDDO               ENDDO
245              ENDDO              ENDDO
# Line 259  c      to avoid using "exit" Line 249  c      to avoid using "exit"
249          ENDDO          ENDDO
250    
251  !        else  !        else
252  #else  ! #else
253    !
254          CALL STREAMICE_CG_ACTION( myThid,  !         CALL STREAMICE_CG_ACTION( myThid,
255       O     Au_SI,  !      O     Au_SI,
256       O     Av_SI,  !      O     Av_SI,
257       I     Du_SI,  !      I     Du_SI,
258       I     Dv_SI,  !      I     Dv_SI,
259       I     is,ie,js,je)  !      I     is,ie,js,je)
260    !
261  !        ENDIF  ! !        ENDIF
262    !
263  #endif  ! #endif
264                
 ! !        if (iter.eq.1) then  
 ! !         CALL WRITE_FLD_XY_RL ("Au2","",Au_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Av2","",Av_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Du","",Du_SI,0,myThid)  
 ! !          
 ! !         CALL WRITE_FLD_XY_RL("Dv","",Au_SI,0,myThid)  
 ! !          
 ! !         CALL WRITE_FLD_XY_RL("DiagU1","",Diagu_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV1","",Diagv_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagU2","",  
 ! !      &    streamice_cg_A1(i,j,bi,bj,0,0),0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV2","",  
 ! !      &    streamice_cg_A4(i,j,bi,bj,0,0),0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL("DiagV2","",Diagv_SI,0,myThid)  
 ! !  
 ! !        endif  
   
   
   
   
   
 ! !        if (iter.eq.1) then  
 ! !         CALL WRITE_FLD_XY_RL ("Au1","",Au_SI,0,myThid)  
 ! !         CALL WRITE_FLD_XY_RL ("Av1","",Av_SI,0,myThid)  
 ! !          
 ! !        endif  
   
   
265    
266         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
267          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 454  c      to avoid using "exit" Line 416  c      to avoid using "exit"
416          _EXCH_XY_RL( cg_Vin, myThid )          _EXCH_XY_RL( cg_Vin, myThid )
417         endif         endif
418    
419    
420         endif         endif
421        enddo ! end of CG loop        enddo ! end of CG loop
422                
# Line 464  c     if iters has reached max_iters the Line 427  c     if iters has reached max_iters the
427         conv_flag = 1         conv_flag = 1
428        ENDIF        ENDIF
429    
430        DO bj = myByLo(myThid), myByHi(myThid)  !       DO bj = myByLo(myThid), myByHi(myThid)
431         DO bi = myBxLo(myThid), myBxHi(myThid)  !        DO bi = myBxLo(myThid), myBxHi(myThid)
432          DO j=1-OLy,sNy+OLy  !         DO j=1-OLy,sNy+OLy
433           DO i=1-OLy,sNx+OLy  !          DO i=1-OLy,sNx+OLy
434            IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
435       &     cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)  !      &     cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
436            IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
437       &     cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)  !      &     cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
438           ENDDO  !          ENDDO
439          ENDDO  !         ENDDO
440         ENDDO  !        ENDDO
441        ENDDO        !       ENDDO      
442    !
443        _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
444        _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
445    
446  #endif  #endif
447        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22