/[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.5 by dgoldberg, Sat Apr 6 17:44:02 2013 UTC
# Line 108  C     LOCAL VARIABLES Line 108  C     LOCAL VARIABLES
108         ENDDO         ENDDO
109        ENDDO        ENDDO
110    
111          
112    
113        CALL STREAMICE_CG_SOLVE(        CALL STREAMICE_CG_SOLVE(
114       &  U_state,       &  U_state,
# Line 146  C     LOCAL VARIABLES Line 146  C     LOCAL VARIABLES
146        _EXCH_XY_RL( U_state, myThid )        _EXCH_XY_RL( U_state, myThid )
147        _EXCH_XY_RL( V_state, myThid )          _EXCH_XY_RL( V_state, myThid )  
148    
 !       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  
149    
150        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
151         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 529  c     if iters has reached max_iters the Line 160  c     if iters has reached max_iters the
160       &           adA_uu(i,j,bi,bj,colx,coly) -       &           adA_uu(i,j,bi,bj,colx,coly) -
161       &           cg_Uin(i,j,bi,bj) *       &           cg_Uin(i,j,bi,bj) *
162       &           U_state(i+colx,j+coly,bi,bj)       &           U_state(i+colx,j+coly,bi,bj)
163  !                 print *,"ADA",cg_Uin(i,j,bi,bj)  
164               endif               endif
165               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then               if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
166                  adA_uv(i,j,bi,bj,colx,coly) =                  adA_uv(i,j,bi,bj,colx,coly) =
# Line 556  c     if iters has reached max_iters the Line 187  c     if iters has reached max_iters the
187    
188             enddo             enddo
189            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  
190           enddo           enddo
191          enddo          enddo
192         enddo         enddo
193        enddo        enddo
194    
 !       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  
           
195    
196    
197        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22