/[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.4 by dgoldberg, Thu Mar 7 15:23:19 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          print *, "GOT HERE 4 mythid=", mythid, tolerance
150    
151        _EXCH_XY_RL( cg_Uin, myThid )        _EXCH_XY_RL( cg_Uin, myThid )
152        _EXCH_XY_RL( cg_Vin, myThid )        _EXCH_XY_RL( cg_Vin, myThid )
153        _EXCH_XY_RL( U_state, myThid )        _EXCH_XY_RL( U_state, myThid )
154        _EXCH_XY_RL( V_state, myThid )          _EXCH_XY_RL( V_state, myThid )
   
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1,sNy  
 !          DO i=1,sNx  
 !           Zu_SI (i,j,bi,bj) = 0. _d 0  
 !           Zv_SI (i,j,bi,bj) = 0. _d 0  
 !           Ru_SI (i,j,bi,bj) = 0. _d 0  
 !           Rv_SI (i,j,bi,bj) = 0. _d 0  
 !           Au_SI (i,j,bi,bj) = 0. _d 0  
 !           Av_SI (i,j,bi,bj) = 0. _d 0  
 !           Du_SI (i,j,bi,bj) = 0. _d 0  
 !           Dv_SI (i,j,bi,bj) = 0. _d 0  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
 !        
 ! C     FIND INITIAL RESIDUAL, and initialize r  
 !  
 ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX  
 !  
 !          
 !  
 !         DO bj = myByLo(myThid), myByHi(myThid)  
 !          DO bi = myBxLo(myThid), myBxHi(myThid)  
 !           DO j=js,je  
 !            DO i=is,ie  
 !             DO colx=-1,1  
 !              DO coly=-1,1  
 !               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +  
 !      &         A_uu(i,j,bi,bj,colx,coly)*  
 !      &         cg_Uin(i+colx,j+coly,bi,bj)+  
 !      &         A_uv(i,j,bi,bj,colx,coly)*      
 !      &         cg_Vin(i+colx,j+coly,bi,bj)  
 !               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +  
 !      &         A_vu(i,j,bi,bj,colx,coly)*  
 !      &         cg_Uin(i+colx,j+coly,bi,bj)+  
 !      &         A_vv(i,j,bi,bj,colx,coly)*      
 !      &         cg_Vin(i+colx,j+coly,bi,bj)  
 !              ENDDO  
 !             ENDDO  
 !            ENDDO  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !  
 !  
 !  
 ! ! #else  
 ! !  
 ! !       CALL STREAMICE_CG_ACTION( myThid,  
 ! !      O    Au_SI,  
 ! !      O    Av_SI,  
 ! !      I    cg_Uin,  
 ! !      I    cg_Vin,  
 ! !      I    0, sNx+1, 0, sNy+1 )  
 ! !  
 ! ! #endif  
 !  
 !       _EXCH_XY_RL( Au_SI, myThid )  
 !       _EXCH_XY_RL( Av_SI, myThid )  
 !  
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1-OLy,sNy+OLy  
 !          DO i=1-OLx,sNx+OLx  
 !           Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-  
 !      &     Au_SI(i,j,bi,bj)  
 !           Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-  
 !      &     Av_SI(i,j,bi,bj)  
 !          ENDDO  
 !         ENDDO  
 !         dot_p1_tile(bi,bj) = 0. _d 0  
 !         dot_p2_tile(bi,bj) = 0. _d 0    
 !        ENDDO  
 !       ENDDO  
 !  
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1,sNy  
 !          DO i=1,sNx  
 !           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)  
 !      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2  
 !           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)  
 !      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
 !  
 !       CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )  
 !       resid_0 = sqrt(dot_p1)  
 !  
 ! C    CCCCCCCCCCCCCCCCCCCC  
 !  
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1-OLy,sNy+OLy  
 !          DO i=1-OLx,sNx+OLx  
 !           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)  
 !      &      Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)  
 !           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)  
 !      &      Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
 !  
 !       cg_halo = min(OLx-1,OLy-1)  
 !       conv_flag = 0  
 !  
 !       DO bj = myByLo(myThid), myByHi(myThid)  
 !        DO bi = myBxLo(myThid), myBxHi(myThid)  
 !         DO j=1-OLy,sNy+OLy  
 !          DO i=1-OLx,sNx+OLx  
 !           Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)  
 !           Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !       ENDDO  
 !  
 !       resid = resid_0  
 !       iters = 0  
 !        
 ! c  !!!!!!!!!!!!!!!!!!  
 ! c  !!              !!  
 ! c  !! MAIN CG LOOP !!  
 ! c  !!              !!  
 ! c  !!!!!!!!!!!!!!!!!!  
 !  
 !  
 !        
 !    
 ! c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)  
 !  
 !        
 !  
 ! !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  
 !  
 !  
 !       do iter = 1, streamice_max_cg_iter  
 !        if (resid .gt. tolerance*resid_0) then  
 !  
 ! c      to avoid using "exit"  
 !        iters = iters + 1  
 !  
 !        is = 1 - cg_halo  
 !        ie = sNx + cg_halo  
 !        js = 1 - cg_halo  
 !        je = sNy + cg_halo  
 !        
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1-OLy,sNy+OLy  
 !           DO i=1-OLx,sNx+OLx  
 !            Au_SI(i,j,bi,bj) = 0. _d 0  
 !            Av_SI(i,j,bi,bj) = 0. _d 0  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 ! !        IF (STREAMICE_construct_matrix) THEN  
 !  
 ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX  
 !  
 !         DO bj = myByLo(myThid), myByHi(myThid)  
 !          DO bi = myBxLo(myThid), myBxHi(myThid)  
 !           DO j=js,je  
 !            DO i=is,ie  
 !             DO colx=-1,1  
 !              DO coly=-1,1  
 !               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +  
 !      &         A_uu(i,j,bi,bj,colx,coly)*  
 !      &         Du_SI(i+colx,j+coly,bi,bj)+  
 !      &         A_uv(i,j,bi,bj,colx,coly)*      
 !      &         Dv_SI(i+colx,j+coly,bi,bj)  
 !               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +  
 !      &         A_vu(i,j,bi,bj,colx,coly)*  
 !      &         Du_SI(i+colx,j+coly,bi,bj)+  
 !      &         A_vv(i,j,bi,bj,colx,coly)*      
 !      &         Dv_SI(i+colx,j+coly,bi,bj)  
 !              ENDDO  
 !             ENDDO  
 !            ENDDO  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !  
 ! !        else  
 ! ! #else  
 ! !  
 ! !         CALL STREAMICE_CG_ACTION( myThid,  
 ! !      O     Au_SI,  
 ! !      O     Av_SI,  
 ! !      I     Du_SI,  
 ! !      I     Dv_SI,  
 ! !      I     is,ie,js,je)  
 ! !  
 ! ! !        ENDIF  
 ! !  
 ! ! #endif  
 !          
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          dot_p1_tile(bi,bj) = 0. _d 0  
 !          dot_p2_tile(bi,bj) = 0. _d 0    
 !         ENDDO  
 !        ENDDO  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1,sNy  
 !           DO i=1,sNx  
 !            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN  
 !            dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*  
 !      &        Ru_SI(i,j,bi,bj)  
 !             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*  
 !      &        Au_SI(i,j,bi,bj)  
 !            ENDIF  
 !            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN  
 !             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*  
 !      &        Rv_SI(i,j,bi,bj)  
 !             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*  
 !      &        Av_SI(i,j,bi,bj)  
 !            ENDIF  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 !        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )  
 !        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )              
 !        alpha_k = dot_p1/dot_p2  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1-OLy,sNy+OLy  
 !           DO i=1-OLx,sNx+OLx  
 !  
 !            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN  
 !             cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+  
 !      &       alpha_k*Du_SI(i,j,bi,bj)  
 !             Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)  
 !             Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)  
 !             Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-  
 !      &       alpha_k*Au_SI(i,j,bi,bj)  
 !             Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /  
 !      &       DIAGu_SI(i,j,bi,bj)  
 !            ENDIF  
 !  
 !            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN  
 !             cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+  
 !      &       alpha_k*Dv_SI(i,j,bi,bj)  
 !             Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)  
 !             Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)  
 !             Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-  
 !      &       alpha_k*Av_SI(i,j,bi,bj)  
 !             Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /  
 !      &       DIAGv_SI(i,j,bi,bj)  
 !  
 !            ENDIF  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          dot_p1_tile(bi,bj) = 0. _d 0  
 !          dot_p2_tile(bi,bj) = 0. _d 0    
 !         ENDDO  
 !        ENDDO  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1,sNy  
 !           DO i=1,sNx  
 !  
 !            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN  
 !             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*  
 !      &        Ru_SI(i,j,bi,bj)  
 !             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*  
 !      &        Ru_old_SI(i,j,bi,bj)  
 !            ENDIF  
 !  
 !            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN  
 !             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*  
 !      &        Rv_SI(i,j,bi,bj)  
 !             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*  
 !      &        Rv_old_SI(i,j,bi,bj)  
 !            ENDIF  
 !  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 !        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )  
 !        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )  
 !                      
 !        beta_k = dot_p1/dot_p2  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1-OLy,sNy+OLy  
 !           DO i=1-OLx,sNx+OLx  
 !            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)  
 !      &      Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+  
 !      &      Zu_SI(i,j,bi,bj)  
 !            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)  
 !      &      Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+  
 !      &      Zv_SI(i,j,bi,bj)  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          dot_p1_tile(bi,bj) = 0. _d 0  
 !         ENDDO  
 !        ENDDO  
 !  
 !        DO bj = myByLo(myThid), myByHi(myThid)  
 !         DO bi = myBxLo(myThid), myBxHi(myThid)  
 !          DO j=1,sNy  
 !           DO i=1,sNx  
 !            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)  
 !      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2  
 !            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)  
 !      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2  
 !           ENDDO  
 !          ENDDO  
 !         ENDDO  
 !        ENDDO  
 !  
 !        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )  
 !        resid = sqrt(dot_p1)  
 !  
 ! !        IF (iter .eq. 1) then  
 ! !         print *, alpha_k, beta_k, resid  
 ! !        ENDIF  
 !  
 !        cg_halo = cg_halo - 1  
 !  
 !        if (cg_halo .eq. 0) then  
 !         cg_halo = min(OLx-1,OLy-1)  
 !         _EXCH_XY_RL( Du_SI, myThid )  
 !         _EXCH_XY_RL( Dv_SI, myThid )  
 !         _EXCH_XY_RL( Ru_SI, myThid )  
 !         _EXCH_XY_RL( Rv_SI, myThid )  
 !         _EXCH_XY_RL( cg_Uin, myThid )  
 !         _EXCH_XY_RL( cg_Vin, myThid )  
 !        endif  
 !  
 !        endif  
 !       enddo ! end of CG loop  
   
         
         
 c     to avoid using "exit"  
 c     if iters has reached max_iters there is no convergence  
         
 !       IF (iters .lt. streamice_max_cg_iter) THEN  
 !        conv_flag = 1  
 !       ENDIF  
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  !                 print *,"ADA",cg_Uin(i,j,bi,bj)  
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 542  c     if iters has reached max_iters the Line 179  c     if iters has reached max_iters the
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 556  c     if iters has reached max_iters the Line 193  c     if iters has reached max_iters the
193    
194             enddo             enddo
195            enddo            enddo
 !           if (i.eq.1.and.j.eq.1) then  
 !            print *, "adA_uu", adA_uu(i,j,bi,bj,-1,-1),  
 !      &                        adA_uu(i,j,bi,bj,-1,0),  
 !      &                        adA_uu(i,j,bi,bj,-1,1),  
 !      &                        cg_Uin(i,j,bi,bj)  
 !           endif  
196           enddo           enddo
197          enddo          enddo
198         enddo         enddo
199        enddo        enddo
200    
 !       print *,"MATRIX 1"  
 !       do i=1,sNx  
 !        do j=1,sNy  
 !         print *, ada_uu(i,j,1,1,-1,-1), ada_uu(i,j,1,1,0,-1),  
 !      &  ada_uu(i,j,1,1,1,-1), ada_uu(i,j,1,1,-1,0),  
 !      &  ada_uu(i,j,1,1,0,0), ada_uu(i,j,1,1,1,0),  
 !      &  ada_uu(i,j,1,1,-1,1), ada_uu(i,j,1,1,0,1),ada_uu(i,j,1,1,1,1)  
 !        enddo  
 !       enddo  
 !  
 !       print *,"MATRIX 2"  
 !       do i=1,sNx  
 !        do j=1,sNy  
 !         print *, ada_uv(i,j,1,1,-1,-1), ada_uv(i,j,1,1,0,-1),  
 !      &  ada_uv(i,j,1,1,1,-1), ada_uv(i,j,1,1,-1,0),  
 !      &  ada_uv(i,j,1,1,0,0), ada_uv(i,j,1,1,1,0),  
 !      &  ada_uv(i,j,1,1,-1,1), ada_uv(i,j,1,1,0,1),ada_uv(i,j,1,1,1,1)  
 !        enddo  
 !       enddo  
 !  
 !       print *,"MATRIX 3"  
 !       do i=1,sNx  
 !        do j=1,sNy  
 !         print *, ada_vu(i,j,1,1,-1,-1), ada_vu(i,j,1,1,0,-1),  
 !      &  ada_vu(i,j,1,1,1,-1), ada_vu(i,j,1,1,-1,0),  
 !      &  ada_vu(i,j,1,1,0,0), ada_vu(i,j,1,1,1,0),  
 !      &  ada_vu(i,j,1,1,-1,1), ada_vu(i,j,1,1,0,1),ada_vu(i,j,1,1,1,1)  
 !        enddo  
 !       enddo  
 !  
 !       print *,"MATRIX 4"  
 !       do i=1,sNx  
 !        do j=1,sNy  
 !         print *, ada_vv(i,j,1,1,-1,-1), ada_vv(i,j,1,1,0,-1),  
 !      &  ada_vv(i,j,1,1,1,-1), ada_vv(i,j,1,1,-1,0),  
 !      &  ada_vv(i,j,1,1,0,0), ada_vv(i,j,1,1,1,0),  
 !      &  ada_vv(i,j,1,1,-1,1), ada_vv(i,j,1,1,0,1),ada_vv(i,j,1,1,1,1)  
 !        enddo  
 !       enddo  
           
   
   
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 627  c     if iters has reached max_iters the Line 216  c     if iters has reached max_iters the
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 646  c     if iters has reached max_iters the Line 234  c     if iters has reached max_iters the
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.4  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22