/[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.5 by dgoldberg, Thu Jul 26 16:13:18 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          CHARACTER*(MAX_LEN_MBUF) msgBuf
60    
61        iters = streamice_max_cg_iter        iters = streamice_max_cg_iter
62    
63  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
64    
65    
66        conv_flag = 0        conv_flag = 0
67    
68    
69        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
70         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
71          DO j=1,sNy          DO j=1,sNy
72           DO i=1,sNx           DO i=1,sNx
73            Zu_SI (i,j,bi,bj) = 0. _d 0            Zu_SI (i,j,bi,bj) = 0. _d 0
74            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  
75            Ru_SI (i,j,bi,bj) = 0. _d 0            Ru_SI (i,j,bi,bj) = 0. _d 0
76            Rv_SI (i,j,bi,bj) = 0. _d 0            Rv_SI (i,j,bi,bj) = 0. _d 0
77            Au_SI (i,j,bi,bj) = 0. _d 0            Au_SI (i,j,bi,bj) = 0. _d 0
78            Av_SI (i,j,bi,bj) = 0. _d 0            Av_SI (i,j,bi,bj) = 0. _d 0
79            Du_SI (i,j,bi,bj) = 0. _d 0            Du_SI (i,j,bi,bj) = 0. _d 0
80            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  
81           ENDDO           ENDDO
82          ENDDO          ENDDO
83         ENDDO         ENDDO
84        ENDDO        ENDDO
85          
86    C     FIND INITIAL RESIDUAL, and initialize r
87    
88  C     PREAMBLE: get bdry contribution, find matrix diagonal,  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
 C     initialize iterates R (residual), Z, and D  
89    
90        CALL STREAMICE_CG_BOUND_VALS( myThid,              
      O    ubd_SI,  
      O    vbd_SI)  
91    
92        DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
93         DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
94          DO j=1-OLy,sNy+OLy            DO j=1,sNy
95           DO i=1-OLx,sNx+OLx             DO i=1,sNx
96            RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)              DO colx=-1,1
97       &     - ubd_SI(i,j,bi,bj)               DO coly=-1,1
98            RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
99       &     - vbd_SI(i,j,bi,bj)       &         A_uu(i,j,bi,bj,colx,coly)*
100           ENDDO       &         cg_Uin(i+colx,j+coly,bi,bj)+
101          ENDDO       &         A_uv(i,j,bi,bj,colx,coly)*    
102         ENDDO       &         cg_Vin(i+colx,j+coly,bi,bj)
       ENDDO  
         
       _EXCH_XY_RL( RHSu_SI, myThid )  
       _EXCH_XY_RL( RHSv_SI, myThid )  
103    
       CALL STREAMICE_CG_ADIAG( myThid,  
      O    DIAGu_SI,  
      O    DIAGv_SI)  
104    
105        _EXCH_XY_RL( DIAGu_SI, myThid )                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
106        _EXCH_XY_RL( DIAGv_SI, myThid )       &         A_vu(i,j,bi,bj,colx,coly)*
107               &         cg_Uin(i+colx,j+coly,bi,bj)+
108  C     FIX PROBLEM WITH PRECOND LATER       &         A_vv(i,j,bi,bj,colx,coly)*    
109         &         cg_Vin(i+colx,j+coly,bi,bj)
110  !       DO bj = myByLo(myThid), myByHi(myThid)               ENDDO
111  !        DO bi = myBxLo(myThid), myBxHi(myThid)              ENDDO
112  !         DO j=1-OLy,sNy+OLy             ENDDO
113  !          DO i=1-OLx,sNx+OLx            ENDDO
114  !           DIAGu_SI(i,j,bi,bj)=1.0           ENDDO
115  !           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 )  
116    
117                
118        _EXCH_XY_RL( Au_SI, myThid )        _EXCH_XY_RL( Au_SI, myThid )
119        _EXCH_XY_RL( Av_SI, myThid )        _EXCH_XY_RL( Av_SI, myThid )
120    
121    
122        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
123         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
124          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
125           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
126            Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-            Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
127       &     Au_SI(i,j,bi,bj)       &     Au_SI(i,j,bi,bj)
128            Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-            Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
129       &     Av_SI(i,j,bi,bj)       &     Av_SI(i,j,bi,bj)
130           ENDDO           ENDDO
131          ENDDO          ENDDO
# Line 143  C     FIX PROBLEM WITH PRECOND LATER Line 134  C     FIX PROBLEM WITH PRECOND LATER
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136    
137          
138    
139        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
140         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
141          DO j=1,sNy          DO j=1,sNy
# Line 159  C     FIX PROBLEM WITH PRECOND LATER Line 152  C     FIX PROBLEM WITH PRECOND LATER
152        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
153        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
154    
155          WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
156         &         resid_0
157           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
158         &                     SQUEEZE_RIGHT , 1)
159    
160    C    CCCCCCCCCCCCCCCCCCCC
161    
162        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
163         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
164          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 172  C     FIX PROBLEM WITH PRECOND LATER Line 172  C     FIX PROBLEM WITH PRECOND LATER
172         ENDDO         ENDDO
173        ENDDO        ENDDO
174    
   
175        cg_halo = min(OLx-1,OLy-1)        cg_halo = min(OLx-1,OLy-1)
176        conv_flag = 0        conv_flag = 0
177    
# Line 195  c  !!              !! Line 194  c  !!              !!
194  c  !! MAIN CG LOOP !!  c  !! MAIN CG LOOP !!
195  c  !!              !!  c  !!              !!
196  c  !!!!!!!!!!!!!!!!!!  c  !!!!!!!!!!!!!!!!!!
197    
198    
199        
200        
201  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!!)
202    
203        print *, "BEGINNING MAIN CG LOOP"        print *, "BEGINNING MAIN CG LOOP"
204    
205  !       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)  
206    
 #endif  
207    
208        do iter = 1, streamice_max_cg_iter        do iter = 1, streamice_max_cg_iter
209         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
# Line 224  c      to avoid using "exit" Line 222  c      to avoid using "exit"
222            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
223             Au_SI(i,j,bi,bj) = 0. _d 0             Au_SI(i,j,bi,bj) = 0. _d 0
224             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  
225            ENDDO            ENDDO
226           ENDDO           ENDDO
227          ENDDO          ENDDO
# Line 233  c      to avoid using "exit" Line 229  c      to avoid using "exit"
229    
230  !        IF (STREAMICE_construct_matrix) THEN  !        IF (STREAMICE_construct_matrix) THEN
231    
232  #ifdef STREAMICE_CONSTRUCT_MATRIX  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
233    
234          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
235           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 242  c      to avoid using "exit" Line 238  c      to avoid using "exit"
238              DO colx=-1,1              DO colx=-1,1
239               DO coly=-1,1               DO coly=-1,1
240                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
241       &         streamice_cg_A1(i,j,bi,bj,colx,coly)*       &         A_uu(i,j,bi,bj,colx,coly)*
242       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
243       &         streamice_cg_A2(i,j,bi,bj,colx,coly)*           &         A_uv(i,j,bi,bj,colx,coly)*    
244       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
245                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
246       &         streamice_cg_A3(i,j,bi,bj,colx,coly)*       &         A_vu(i,j,bi,bj,colx,coly)*
247       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
248       &         streamice_cg_A4(i,j,bi,bj,colx,coly)*           &         A_vv(i,j,bi,bj,colx,coly)*    
249       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
250               ENDDO               ENDDO
251              ENDDO              ENDDO
# Line 259  c      to avoid using "exit" Line 255  c      to avoid using "exit"
255          ENDDO          ENDDO
256    
257  !        else  !        else
258  #else  ! #else
259    !
260          CALL STREAMICE_CG_ACTION( myThid,  !         CALL STREAMICE_CG_ACTION( myThid,
261       O     Au_SI,  !      O     Au_SI,
262       O     Av_SI,  !      O     Av_SI,
263       I     Du_SI,  !      I     Du_SI,
264       I     Dv_SI,  !      I     Dv_SI,
265       I     is,ie,js,je)  !      I     is,ie,js,je)
266    !
267  !        ENDIF  ! !        ENDIF
268    !
269  #endif  ! #endif
270                
 ! !        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  
   
   
271    
272         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
273          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 454  c      to avoid using "exit" Line 422  c      to avoid using "exit"
422          _EXCH_XY_RL( cg_Vin, myThid )          _EXCH_XY_RL( cg_Vin, myThid )
423         endif         endif
424    
425    
426         endif         endif
427        enddo ! end of CG loop        enddo ! end of CG loop
428                
# Line 464  c     if iters has reached max_iters the Line 433  c     if iters has reached max_iters the
433         conv_flag = 1         conv_flag = 1
434        ENDIF        ENDIF
435    
436        DO bj = myByLo(myThid), myByHi(myThid)  !       DO bj = myByLo(myThid), myByHi(myThid)
437         DO bi = myBxLo(myThid), myBxHi(myThid)  !        DO bi = myBxLo(myThid), myBxHi(myThid)
438          DO j=1-OLy,sNy+OLy  !         DO j=1-OLy,sNy+OLy
439           DO i=1-OLy,sNx+OLy  !          DO i=1-OLy,sNx+OLy
440            IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
441       &     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)
442            IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
443       &     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)
444           ENDDO  !          ENDDO
445          ENDDO  !         ENDDO
446         ENDDO  !        ENDDO
447        ENDDO        !       ENDDO      
448    !
449        _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
450        _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
451    
452  #endif  #endif
453        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22