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

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

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

revision 1.5 by dgoldberg, Sat Apr 6 17:44:02 2013 UTC revision 1.6 by dgoldberg, Wed Aug 27 19:29:12 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 ADSTREAMICE_CG_SOLVE(        SUBROUTINE ADSTREAMICE_CG_SOLVE(
10       U                               U_state,    ! velocities - need to be recalc'ed       U                               U_state,    ! velocities - need to be recalc ed
11       I                               cg_Bu,      ! adjoint of vel (input)       I                               cg_Bu,      ! adjoint of vel (input)
12       U                               V_state,    ! velocities - need to be recalc'ed       U                               V_state,    ! velocities - need to be recalc ed
13       I                               cg_Bv,      ! adjoint of vel (input)       I                               cg_Bv,      ! adjoint of vel (input)
14       I                               Bu_state,   ! to recalc velocities       I                               Bu_state,   ! to recalc velocities
15       U                               cg_Uin,     ! adjoint of RHS (output)       U                               cg_Uin,     ! adjoint of RHS (output)
# Line 23  CBOP Line 23  CBOP
23       U                               adA_vu,     ! adjoint of matrix coeffs (output)       U                               adA_vu,     ! adjoint of matrix coeffs (output)
24       I                               A_vv,       ! section of matrix that multiplies v and projects on v       I                               A_vv,       ! section of matrix that multiplies v and projects on v
25       U                               adA_vv,     ! adjoint of matrix coeffs (output)       U                               adA_vv,     ! adjoint of matrix coeffs (output)
26       I                               tolerance,       I                               tolerance,
27       I                               iters,       I                               maxiters,
28       I                               myThid )       I                               myThid )
29  C     /============================================================\  C     *============================================================*
30  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |
31  C     | o                                                          |  C     | o                                                          |
32  C     |============================================================|  C     *============================================================*
33  C     |                                                            |  
34  C     \============================================================/  C     !USES:
35        IMPLICIT NONE        IMPLICIT NONE
36    
37  C     === Global variables ===  C     === Global variables ===
# Line 41  C     === Global variables === Line 41  C     === Global variables ===
41  #include "STREAMICE.h"  #include "STREAMICE.h"
42  #include "STREAMICE_CG.h"  #include "STREAMICE_CG.h"
43    
         
44  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
45  C     cg_Uin, cg_Vin - input and output velocities  C     cg_Uin, cg_Vin - input and output velocities
46  C     cg_Bu, cg_Bv - driving stress  C     cg_Bu, cg_Bv - driving stress
# Line 53  C     cg_Bu, cg_Bv - driving stress Line 52  C     cg_Bu, cg_Bv - driving stress
52        _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
53        _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54        _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55        _RL        _RL
56       & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),       & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
57       & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),       & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
58       & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),       & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
# Line 63  C     cg_Bu, cg_Bv - driving stress Line 62  C     cg_Bu, cg_Bv - driving stress
62       & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),       & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
63       & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)       & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
64        _RL tolerance        _RL tolerance
65        INTEGER iters        INTEGER maxiters
66        INTEGER myThid        INTEGER myThid
67    
68  C     LOCAL VARIABLES  C     !LOCAL VARIABLES
69        INTEGER i, j, bi, bj, cg_halo, conv_flag, tmpiter        INTEGER i, j, bi, bj, cg_halo, conv_flag, tmpiter
70        INTEGER iter, is, js, ie, je, colx, coly, k        INTEGER iter, is, js, ie, je, colx, coly, k
71        _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
# Line 78  C     LOCAL VARIABLES Line 77  C     LOCAL VARIABLES
77        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
78        _RL ad_tolerance        _RL ad_tolerance
79        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
80    CEOP
81    
82  !       iters = streamice_max_cg_iter  !       iters = streamice_max_cg_iter
83    
# Line 86  C     LOCAL VARIABLES Line 86  C     LOCAL VARIABLES
86        WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'        WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'
87         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
88       &                     SQUEEZE_RIGHT , 1)       &                     SQUEEZE_RIGHT , 1)
89          
90          print *, "GOT HERE mythid=", mythid, tolerance
91    
92        conv_flag = 0        conv_flag = 0
93        ad_tolerance = 1.e-14        ad_tolerance = 1.e-14
94    
95        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
96         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
97          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
98           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
99            Utemp (i,j,bi,bj) =            Utemp (i,j,bi,bj) =
100       &     cg_Uin (i,j,bi,bj)       &     cg_Uin (i,j,bi,bj)
101            Vtemp (i,j,bi,bj) =            Vtemp (i,j,bi,bj) =
# Line 108  C     LOCAL VARIABLES Line 109  C     LOCAL VARIABLES
109         ENDDO         ENDDO
110        ENDDO        ENDDO
111    
112                print *, "GOT HERE 2 mythid=", mythid, tolerance
113    
114        CALL STREAMICE_CG_SOLVE(        CALL STREAMICE_CG_SOLVE(
115       &  U_state,       &  U_state,
116       &  V_state,       &  V_state,
117       &  Bu_state,       &  Bu_state,
# Line 118  C     LOCAL VARIABLES Line 119  C     LOCAL VARIABLES
119       &  A_uu,       &  A_uu,
120       &  A_uv,       &  A_uv,
121       &  A_vu,       &  A_vu,
122       &  A_vv,       &  A_vv,
123       &  tolerance,       &  tolerance,
124       &  tmpiter,       &  tmpiter,
125         &  maxiters,
126       &  myThid )       &  myThid )
127    
128          print *, "GOT HERE 3 mythid=", mythid, tolerance
129    
130        tmpiter = 0        tmpiter = 0
131    
132        _EXCH_XY_RL( cg_Bu, myThid )        _EXCH_XY_RL( cg_Bu, myThid )
133        _EXCH_XY_RL( cg_Bv, myThid )        _EXCH_XY_RL( cg_Bv, myThid )
134    
135        CALL STREAMICE_CG_SOLVE(        CALL STREAMICE_CG_SOLVE(
136       &  cg_Uin,       &  cg_Uin,
137       &  cg_Vin,       &  cg_Vin,
138       &  cg_Bu,       &  cg_Bu,
# Line 136  C     LOCAL VARIABLES Line 140  C     LOCAL VARIABLES
140       &  A_uu,       &  A_uu,
141       &  A_uv,       &  A_uv,
142       &  A_vu,       &  A_vu,
143       &  A_vv,       &  A_vv,
144       &  ad_tolerance,       &  ad_tolerance,
145       &  tmpiter,       &  tmpiter,
146         &  maxiters,
147       &  myThid )       &  myThid )
148    
149        _EXCH_XY_RL( cg_Uin, myThid )        print *, "GOT HERE 4 mythid=", mythid, tolerance
       _EXCH_XY_RL( cg_Vin, myThid )  
       _EXCH_XY_RL( U_state, myThid )  
       _EXCH_XY_RL( V_state, myThid )    
150    
151          _EXCH_XY_RL( cg_Uin, myThid )
152          _EXCH_XY_RL( cg_Vin, myThid )
153          _EXCH_XY_RL( U_state, myThid )
154          _EXCH_XY_RL( V_state, myThid )
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,sNy
159           DO i=1-Olx,sNx+Olx           DO i=1,sNx
160            DO colx=-1,1            DO colx=-1,1
161             DO coly=-1,1             DO coly=-1,1
162    
163              if (STREAMICE_umask(i,j,bi,bj).eq.1) then              if (STREAMICE_umask(i,j,bi,bj).eq.1) then
164               if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then               if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
165                  adA_uu(i,j,bi,bj,colx,coly) =                  adA_uu(i,j,bi,bj,colx,coly) =
166       &           adA_uu(i,j,bi,bj,colx,coly) -       &           adA_uu(i,j,bi,bj,colx,coly) -
167       &           cg_Uin(i,j,bi,bj) *       &           cg_Uin(i,j,bi,bj) *
168       &           U_state(i+colx,j+coly,bi,bj)       &           U_state(i+colx,j+coly,bi,bj)
169    
170               endif               endif
171               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
172                  adA_uv(i,j,bi,bj,colx,coly) =                  adA_uv(i,j,bi,bj,colx,coly) =
173       &           adA_uv(i,j,bi,bj,colx,coly) -       &           adA_uv(i,j,bi,bj,colx,coly) -
174       &           cg_Uin(i,j,bi,bj) *       &           cg_Uin(i,j,bi,bj) *
175       &           V_state(i+colx,j+coly,bi,bj)       &           V_state(i+colx,j+coly,bi,bj)
176               endif               endif
# Line 173  C     LOCAL VARIABLES Line 179  C     LOCAL VARIABLES
179              if (STREAMICE_vmask(i,j,bi,bj).eq.1) then              if (STREAMICE_vmask(i,j,bi,bj).eq.1) then
180               if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then               if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
181                  adA_vu(i,j,bi,bj,colx,coly) =                  adA_vu(i,j,bi,bj,colx,coly) =
182       &           adA_vu(i,j,bi,bj,colx,coly) -       &           adA_vu(i,j,bi,bj,colx,coly) -
183       &           cg_Vin(i,j,bi,bj) *       &           cg_Vin(i,j,bi,bj) *
184       &           U_state(i+colx,j+coly,bi,bj)       &           U_state(i+colx,j+coly,bi,bj)
185               endif               endif
186               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
187                  adA_vv(i,j,bi,bj,colx,coly) =                  adA_vv(i,j,bi,bj,colx,coly) =
188       &           adA_vv(i,j,bi,bj,colx,coly) -       &           adA_vv(i,j,bi,bj,colx,coly) -
189       &           cg_Vin(i,j,bi,bj) *       &           cg_Vin(i,j,bi,bj) *
190       &           V_state(i+colx,j+coly,bi,bj)       &           V_state(i+colx,j+coly,bi,bj)
191               endif               endif
# Line 192  C     LOCAL VARIABLES Line 198  C     LOCAL VARIABLES
198         enddo         enddo
199        enddo        enddo
200    
   
   
201        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
202         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
203          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
204           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
205            if (i.lt.1.or.i.gt.sNx.or.            if (i.lt.1.or.i.gt.sNx.or.
206       &        j.lt.1.or.j.gt.sNy) then       &        j.lt.1.or.j.gt.sNy) then
207             cg_Uin (i,j,bi,bj) = 0.0             cg_Uin (i,j,bi,bj) = 0.0
208             cg_Vin (i,j,bi,bj) = 0.0             cg_Vin (i,j,bi,bj) = 0.0
209              
210             DO colx=-1,1             DO colx=-1,1
211              DO coly=-1,1                                                DO coly=-1,1
212               ada_uu(i,j,bi,bj,colx,coly)=0.0               ada_uu(i,j,bi,bj,colx,coly)=0.0
213               ada_uv(i,j,bi,bj,colx,coly)=0.0               ada_uv(i,j,bi,bj,colx,coly)=0.0
214               ada_vu(i,j,bi,bj,colx,coly)=0.0               ada_vu(i,j,bi,bj,colx,coly)=0.0
# Line 212  C     LOCAL VARIABLES Line 216  C     LOCAL VARIABLES
216              enddo              enddo
217             enddo             enddo
218    
                                     
219            endif            endif
220            cg_Uin (i,j,bi,bj) =            cg_Uin (i,j,bi,bj) =
221       &     cg_Uin (i,j,bi,bj) +       &     cg_Uin (i,j,bi,bj) +
222       &     Utemp (i,j,bi,bj)       &     Utemp (i,j,bi,bj)
223            cg_Vin (i,j,bi,bj) =            cg_Vin (i,j,bi,bj) =
224       &     cg_Vin (i,j,bi,bj) +       &     cg_Vin (i,j,bi,bj) +
225       &     Vtemp (i,j,bi,bj)       &     Vtemp (i,j,bi,bj)
226            cg_bu (i,j,bi,bj) = 0.            cg_bu (i,j,bi,bj) = 0.
227            cg_bv (i,j,bi,bj) = 0.            cg_bv (i,j,bi,bj) = 0.
# Line 231  C     LOCAL VARIABLES Line 234  C     LOCAL VARIABLES
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
   
237        WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter,        WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter,
238       & 'ITERS'       & 'ITERS'
239         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
240       &                     SQUEEZE_RIGHT , 1)       &                     SQUEEZE_RIGHT , 1)
         
241    
242  #endif  #endif
243        RETURN        RETURN
244        END        END
245    
   
         
   
   
   
   
   
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22