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

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

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


Revision 1.5 - (hide annotations) (download)
Sat Apr 6 17:44:02 2013 UTC (12 years, 3 months ago) by dgoldberg
Branch: MAIN
Changes since 1.4: +3 -418 lines
clean up commented code

1 dgoldberg 1.5 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/adstreamice_cg_solve.F,v 1.4 2013/03/07 15:23:19 dgoldberg Exp $
2 dgoldberg 1.1 C $Name: $
3    
4     #include "STREAMICE_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9     SUBROUTINE ADSTREAMICE_CG_SOLVE(
10     U U_state, ! velocities - need to be recalc'ed
11     I cg_Bu, ! adjoint of vel (input)
12     U V_state, ! velocities - need to be recalc'ed
13     I cg_Bv, ! adjoint of vel (input)
14     I Bu_state, ! to recalc velocities
15     U cg_Uin, ! adjoint of RHS (output)
16     I Bv_state, ! to recalc velocities
17     U cg_Vin, ! adjoint of RHS (output)
18     I A_uu, ! section of matrix that multiplies u and projects on u
19     U adA_uu, ! adjoint of matrix coeffs (output)
20     I A_uv, ! section of matrix that multiplies v and projects on u
21     U adA_uv, ! adjoint of matrix coeffs (output)
22     I A_vu, ! section of matrix that multiplies u and projects on v
23     U adA_vu, ! adjoint of matrix coeffs (output)
24     I A_vv, ! section of matrix that multiplies v and projects on v
25     U adA_vv, ! adjoint of matrix coeffs (output)
26     I tolerance,
27     I iters,
28     I myThid )
29     C /============================================================\
30     C | SUBROUTINE |
31     C | o |
32     C |============================================================|
33     C | |
34     C \============================================================/
35     IMPLICIT NONE
36    
37     C === Global variables ===
38     #include "SIZE.h"
39     #include "EEPARAMS.h"
40     #include "PARAMS.h"
41     #include "STREAMICE.h"
42     #include "STREAMICE_CG.h"
43    
44    
45     C !INPUT/OUTPUT ARGUMENTS
46     C cg_Uin, cg_Vin - input and output velocities
47     C cg_Bu, cg_Bv - driving stress
48     _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49     _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50     _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51     _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52     _RL U_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
53     _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54     _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55     _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56     _RL
57     & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
58     & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
59     & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
60     & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
61     & adA_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
62     & adA_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
63     & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
64     & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
65     _RL tolerance
66     INTEGER iters
67     INTEGER myThid
68    
69     C LOCAL VARIABLES
70     INTEGER i, j, bi, bj, cg_halo, conv_flag, tmpiter
71     INTEGER iter, is, js, ie, je, colx, coly, k
72     _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73     _RL Vtemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
74     _RL UtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
75     _RL VtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76     _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
77     _RL dot_p1_tile (nSx,nSy)
78     _RL dot_p2_tile (nSx,nSy)
79     _RL ad_tolerance
80 dgoldberg 1.3 CHARACTER*(MAX_LEN_MBUF) msgBuf
81 dgoldberg 1.1
82     ! iters = streamice_max_cg_iter
83    
84     #ifdef ALLOW_STREAMICE
85    
86 dgoldberg 1.3 WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'
87     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
88     & SQUEEZE_RIGHT , 1)
89    
90 dgoldberg 1.1
91     conv_flag = 0
92     ad_tolerance = 1.e-14
93    
94     DO bj = myByLo(myThid), myByHi(myThid)
95     DO bi = myBxLo(myThid), myBxHi(myThid)
96     DO j=1-Oly,sNy+Oly
97     DO i=1-Olx,sNx+Olx
98     Utemp (i,j,bi,bj) =
99     & cg_Uin (i,j,bi,bj)
100     Vtemp (i,j,bi,bj) =
101     & cg_Vin (i,j,bi,bj)
102     UtempSt (i,j,bi,bj) =
103     & U_state (i,j,bi,bj)
104     VtempSt (i,j,bi,bj) =
105     & V_state (i,j,bi,bj)
106     ENDDO
107     ENDDO
108     ENDDO
109     ENDDO
110    
111 dgoldberg 1.5
112 dgoldberg 1.2
113 dgoldberg 1.1 CALL STREAMICE_CG_SOLVE(
114     & U_state,
115     & V_state,
116     & Bu_state,
117     & Bv_state,
118     & A_uu,
119     & A_uv,
120     & A_vu,
121     & A_vv,
122     & tolerance,
123     & tmpiter,
124     & myThid )
125    
126     tmpiter = 0
127    
128 dgoldberg 1.2 _EXCH_XY_RL( cg_Bu, myThid )
129     _EXCH_XY_RL( cg_Bv, myThid )
130    
131 dgoldberg 1.1 CALL STREAMICE_CG_SOLVE(
132     & cg_Uin,
133     & cg_Vin,
134     & cg_Bu,
135     & cg_Bv,
136     & A_uu,
137     & A_uv,
138     & A_vu,
139     & A_vv,
140     & ad_tolerance,
141     & tmpiter,
142     & myThid )
143    
144 dgoldberg 1.2 _EXCH_XY_RL( cg_Uin, myThid )
145     _EXCH_XY_RL( cg_Vin, myThid )
146     _EXCH_XY_RL( U_state, myThid )
147     _EXCH_XY_RL( V_state, myThid )
148    
149 dgoldberg 1.1
150     DO bj = myByLo(myThid), myByHi(myThid)
151     DO bi = myBxLo(myThid), myBxHi(myThid)
152     DO j=1-Oly,sNy+Oly
153     DO i=1-Olx,sNx+Olx
154     DO colx=-1,1
155     DO coly=-1,1
156    
157     if (STREAMICE_umask(i,j,bi,bj).eq.1) then
158     if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
159     adA_uu(i,j,bi,bj,colx,coly) =
160     & adA_uu(i,j,bi,bj,colx,coly) -
161     & cg_Uin(i,j,bi,bj) *
162     & U_state(i+colx,j+coly,bi,bj)
163 dgoldberg 1.5
164 dgoldberg 1.1 endif
165     if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
166     adA_uv(i,j,bi,bj,colx,coly) =
167     & adA_uv(i,j,bi,bj,colx,coly) -
168     & cg_Uin(i,j,bi,bj) *
169     & V_state(i+colx,j+coly,bi,bj)
170     endif
171     endif
172    
173     if (STREAMICE_vmask(i,j,bi,bj).eq.1) then
174     if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
175     adA_vu(i,j,bi,bj,colx,coly) =
176     & adA_vu(i,j,bi,bj,colx,coly) -
177     & cg_Vin(i,j,bi,bj) *
178     & U_state(i+colx,j+coly,bi,bj)
179     endif
180     if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
181     adA_vv(i,j,bi,bj,colx,coly) =
182     & adA_vv(i,j,bi,bj,colx,coly) -
183     & cg_Vin(i,j,bi,bj) *
184     & V_state(i+colx,j+coly,bi,bj)
185     endif
186     endif
187    
188     enddo
189     enddo
190     enddo
191     enddo
192     enddo
193     enddo
194    
195    
196    
197     DO bj = myByLo(myThid), myByHi(myThid)
198     DO bi = myBxLo(myThid), myBxHi(myThid)
199     DO j=1-Oly,sNy+Oly
200     DO i=1-Olx,sNx+Olx
201     if (i.lt.1.or.i.gt.sNx.or.
202     & j.lt.1.or.j.gt.sNy) then
203     cg_Uin (i,j,bi,bj) = 0.0
204     cg_Vin (i,j,bi,bj) = 0.0
205    
206     DO colx=-1,1
207     DO coly=-1,1
208     ada_uu(i,j,bi,bj,colx,coly)=0.0
209     ada_uv(i,j,bi,bj,colx,coly)=0.0
210     ada_vu(i,j,bi,bj,colx,coly)=0.0
211     ada_vv(i,j,bi,bj,colx,coly)=0.0
212     enddo
213     enddo
214    
215    
216     endif
217     cg_Uin (i,j,bi,bj) =
218     & cg_Uin (i,j,bi,bj) +
219     & Utemp (i,j,bi,bj)
220     cg_Vin (i,j,bi,bj) =
221     & cg_Vin (i,j,bi,bj) +
222     & Vtemp (i,j,bi,bj)
223     cg_bu (i,j,bi,bj) = 0.
224     cg_bv (i,j,bi,bj) = 0.
225     U_state (i,j,bi,bj) =
226     & UtempSt (i,j,bi,bj)
227     V_state (i,j,bi,bj) =
228     & VtempSt (i,j,bi,bj)
229     ENDDO
230     ENDDO
231     ENDDO
232     ENDDO
233    
234    
235 dgoldberg 1.3 WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter,
236     & 'ITERS'
237     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
238     & SQUEEZE_RIGHT , 1)
239    
240 dgoldberg 1.1
241     #endif
242     RETURN
243     END
244    
245    
246    
247    
248    
249    
250    
251    

  ViewVC Help
Powered by ViewVC 1.1.22