/[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.6 by dgoldberg, Thu Jul 26 16:22:58 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"        WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
204           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
205         &                     SQUEEZE_RIGHT , 1)
206    
207  !       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)  
208    
 #endif  
209    
210        do iter = 1, streamice_max_cg_iter        do iter = 1, streamice_max_cg_iter
211         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
# Line 224  c      to avoid using "exit" Line 224  c      to avoid using "exit"
224            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
225             Au_SI(i,j,bi,bj) = 0. _d 0             Au_SI(i,j,bi,bj) = 0. _d 0
226             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  
227            ENDDO            ENDDO
228           ENDDO           ENDDO
229          ENDDO          ENDDO
# Line 233  c      to avoid using "exit" Line 231  c      to avoid using "exit"
231    
232  !        IF (STREAMICE_construct_matrix) THEN  !        IF (STREAMICE_construct_matrix) THEN
233    
234  #ifdef STREAMICE_CONSTRUCT_MATRIX  ! #ifdef STREAMICE_CONSTRUCT_MATRIX
235    
236          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
237           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 242  c      to avoid using "exit" Line 240  c      to avoid using "exit"
240              DO colx=-1,1              DO colx=-1,1
241               DO coly=-1,1               DO coly=-1,1
242                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +                Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
243       &         streamice_cg_A1(i,j,bi,bj,colx,coly)*       &         A_uu(i,j,bi,bj,colx,coly)*
244       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
245       &         streamice_cg_A2(i,j,bi,bj,colx,coly)*           &         A_uv(i,j,bi,bj,colx,coly)*    
246       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
247                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +                Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
248       &         streamice_cg_A3(i,j,bi,bj,colx,coly)*       &         A_vu(i,j,bi,bj,colx,coly)*
249       &         Du_SI(i+colx,j+coly,bi,bj)+       &         Du_SI(i+colx,j+coly,bi,bj)+
250       &         streamice_cg_A4(i,j,bi,bj,colx,coly)*           &         A_vv(i,j,bi,bj,colx,coly)*    
251       &         Dv_SI(i+colx,j+coly,bi,bj)       &         Dv_SI(i+colx,j+coly,bi,bj)
252               ENDDO               ENDDO
253              ENDDO              ENDDO
# Line 259  c      to avoid using "exit" Line 257  c      to avoid using "exit"
257          ENDDO          ENDDO
258    
259  !        else  !        else
260  #else  ! #else
261    !
262          CALL STREAMICE_CG_ACTION( myThid,  !         CALL STREAMICE_CG_ACTION( myThid,
263       O     Au_SI,  !      O     Au_SI,
264       O     Av_SI,  !      O     Av_SI,
265       I     Du_SI,  !      I     Du_SI,
266       I     Dv_SI,  !      I     Dv_SI,
267       I     is,ie,js,je)  !      I     is,ie,js,je)
268    !
269  !        ENDIF  ! !        ENDIF
270    !
271  #endif  ! #endif
272                
 ! !        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  
   
   
273    
274         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
275          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 454  c      to avoid using "exit" Line 424  c      to avoid using "exit"
424          _EXCH_XY_RL( cg_Vin, myThid )          _EXCH_XY_RL( cg_Vin, myThid )
425         endif         endif
426    
427    
428         endif         endif
429        enddo ! end of CG loop        enddo ! end of CG loop
430                
# Line 464  c     if iters has reached max_iters the Line 435  c     if iters has reached max_iters the
435         conv_flag = 1         conv_flag = 1
436        ENDIF        ENDIF
437    
438        DO bj = myByLo(myThid), myByHi(myThid)  !       DO bj = myByLo(myThid), myByHi(myThid)
439         DO bi = myBxLo(myThid), myBxHi(myThid)  !        DO bi = myBxLo(myThid), myBxHi(myThid)
440          DO j=1-OLy,sNy+OLy  !         DO j=1-OLy,sNy+OLy
441           DO i=1-OLy,sNx+OLy  !          DO i=1-OLy,sNx+OLy
442            IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
443       &     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)
444            IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)  !           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
445       &     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)
446           ENDDO  !          ENDDO
447          ENDDO  !         ENDDO
448         ENDDO  !        ENDDO
449        ENDDO        !       ENDDO      
450    !
451        _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
452        _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
453    
454  #endif  #endif
455        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22