Back to home page

MITgcm
 
 

    

File indexing completed on 2020-04-22 05:11:27 UTC

view on githubraw file Latest commit 07e78522 on 2020-04-21 13:33:29 UTC
5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 
                0005 CBOP
07e785229e dngo*0006       SUBROUTINE STREAMICE_CG_SOLVE(
5ca83cd8f7 Dani*0007      U                               cg_Uin,     ! x-velocities
                0008      U                               cg_Vin,     ! y-velocities
                0009      I                               cg_Bu,      ! force in x dir
                0010      I                               cg_Bv,      ! force in y dir
                0011      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0012      I                               A_uv,       ! section of matrix that multiplies v and projects on u
                0013      I                               A_vu,       ! section of matrix that multiplies u and projects on v
                0014      I                               A_vv,       ! section of matrix that multiplies v and projects on v
07e785229e dngo*0015      I                               tolerance,
5ca83cd8f7 Dani*0016      O                               iters,
d2cdb9260d Dani*0017      I                               maxIter,
5ca83cd8f7 Dani*0018      I                               myThid )
                0019 C     /============================================================\
07e785229e dngo*0020 C     | SUBROUTINE                                                 |
5ca83cd8f7 Dani*0021 C     | o                                                          |
                0022 C     |============================================================|
                0023 C     |                                                            |
                0024 C     \============================================================/
                0025       IMPLICIT NONE
                0026 
                0027 #include "SIZE.h"
                0028 #include "EEPARAMS.h"
                0029 #include "PARAMS.h"
                0030 #include "STREAMICE.h"
                0031 #include "STREAMICE_CG.h"
                0032 
07e785229e dngo*0033 c#ifdef ALLOW_PETSC
                0034 c#include "finclude/petsc.h"
                0035 c UNCOMMENT IF V3.0
                0036 c#include "finclude/petscvec.h"
                0037 c#include "finclude/petscmat.h"
                0038 c#include "finclude/petscksp.h"
                0039 c#include "finclude/petscpc.h"
                0040 c#endif
5ca83cd8f7 Dani*0041 C     === Global variables ===
                0042 
                0043 C     !INPUT/OUTPUT ARGUMENTS
                0044 C     cg_Uin, cg_Vin - input and output velocities
                0045 C     cg_Bu, cg_Bv - driving stress
                0046       INTEGER myThid
                0047       INTEGER iters
d2cdb9260d Dani*0048       INTEGER maxIter
5ca83cd8f7 Dani*0049       _RL tolerance
                0050       _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0051       _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052       _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
07e785229e dngo*0054       _RL
5ca83cd8f7 Dani*0055      & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0056      & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0057      & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0058      & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
                0059 
                0060 C     LOCAL VARIABLES
                0061       INTEGER i, j, bi, bj, cg_halo, conv_flag
07e785229e dngo*0062       INTEGER iter, is, js, ie, je, colx, coly
5ca83cd8f7 Dani*0063       _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
                0064       _RL dot_p1_tile (nSx,nSy)
                0065       _RL dot_p2_tile (nSx,nSy)
                0066       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0067 
07e785229e dngo*0068 c#ifdef ALLOW_PETSC
                0069 c      INTEGER indices(2*(sNx*nSx*sNy*nSy))
                0070 c      INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
                0071 c      _RL rhs_values(2*(sNx*nSx*sNy*nSy))
                0072 c      _RL solution_values(2*(sNx*nSx*sNy*nSy))
                0073 c      _RL mat_values (2*Nx*Ny,2*(sNx*nSx*sNy*nSy))
                0074 c      _RL mat_values (18,1), mat_val_return(1)
                0075 c      INTEGER indices_col(18)
                0076 c      INTEGER local_dofs, global_dofs, dof_index, dof_index_col
                0077 c      INTEGER local_offset
                0078 c      Mat matrix
                0079 c      KSP ksp
                0080 c      PC  pc
                0081 c      Vec rhs
                0082 c      Vec solution
                0083 c      PetscErrorCode ierr
                0084 c#ifdef ALLOW_USE_MPI
                0085 c      integer mpiRC, mpiMyWid
                0086 c#endif
                0087 c#endif
5ca83cd8f7 Dani*0088 
                0089 #ifdef ALLOW_STREAMICE
                0090 
                0091       CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
                0092 #ifndef STREAMICE_SERIAL_TRISOLVE
                0093 
                0094 #ifdef ALLOW_PETSC
                0095 
07e785229e dngo*0096 #ifndef ALLOW_OPENAD
                0097       if (streamice_use_petsc.and.petscFlag.eq.1) then
                0098 #else
18a089944d Dani*0099       if (streamice_use_petsc) then
07e785229e dngo*0100 #endif
18a089944d Dani*0101 
67b4a2b594 Dani*0102       CALL STREAMICE_CG_SOLVE_PETSC(
                0103      U         cg_Uin,     ! x-velocities
                0104      U         cg_Vin,     ! y-velocities
                0105      I         cg_Bu,      ! force in x dir
                0106      I         cg_Bv,      ! force in y dir
                0107      I         A_uu,       ! section of matrix that multiplies u and projects on u
                0108      I         A_uv,       ! section of matrix that multiplies v and projects on u
                0109      I         A_vu,       ! section of matrix that multiplies u and projects on v
                0110      I         A_vv,       ! section of matrix that multiplies v and projects on v
                0111      I         tolerance,
78b0c9990f Dani*0112      I         iters,
                0113      O         maxiter,
67b4a2b594 Dani*0114      I         myThid )
5ca83cd8f7 Dani*0115 
18a089944d Dani*0116       else
                0117 
                0118 #endif  /* ALLOW_PETSC */
5ca83cd8f7 Dani*0119 
d2cdb9260d Dani*0120       iters = maxIter
5ca83cd8f7 Dani*0121       conv_flag = 0
                0122 
                0123       DO bj = myByLo(myThid), myByHi(myThid)
                0124        DO bi = myBxLo(myThid), myBxHi(myThid)
                0125         DO j=1,sNy
                0126          DO i=1,sNx
                0127           Zu_SI (i,j,bi,bj) = 0. _d 0
                0128           Zv_SI (i,j,bi,bj) = 0. _d 0
                0129           Ru_SI (i,j,bi,bj) = 0. _d 0
                0130           Rv_SI (i,j,bi,bj) = 0. _d 0
                0131           Au_SI (i,j,bi,bj) = 0. _d 0
                0132           Av_SI (i,j,bi,bj) = 0. _d 0
                0133           Du_SI (i,j,bi,bj) = 0. _d 0
                0134           Dv_SI (i,j,bi,bj) = 0. _d 0
                0135          ENDDO
                0136         ENDDO
                0137        ENDDO
                0138       ENDDO
                0139 
07e785229e dngo*0140 C     FIND INITIAL RESIDUAL, and initialize r
5ca83cd8f7 Dani*0141 
07e785229e dngo*0142 c #ifdef STREAMICE_CONSTRUCT_MATRIX
5ca83cd8f7 Dani*0143 
                0144         DO bj = myByLo(myThid), myByHi(myThid)
                0145          DO bi = myBxLo(myThid), myBxHi(myThid)
                0146           DO j=1,sNy
                0147            DO i=1,sNx
                0148             DO colx=-1,1
                0149              DO coly=-1,1
07e785229e dngo*0150               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
5ca83cd8f7 Dani*0151      &         A_uu(i,j,bi,bj,colx,coly)*
                0152      &         cg_Uin(i+colx,j+coly,bi,bj)+
07e785229e dngo*0153      &         A_uv(i,j,bi,bj,colx,coly)*
5ca83cd8f7 Dani*0154      &         cg_Vin(i+colx,j+coly,bi,bj)
                0155 
07e785229e dngo*0156               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
5ca83cd8f7 Dani*0157      &         A_vu(i,j,bi,bj,colx,coly)*
                0158      &         cg_Uin(i+colx,j+coly,bi,bj)+
07e785229e dngo*0159      &         A_vv(i,j,bi,bj,colx,coly)*
5ca83cd8f7 Dani*0160      &         cg_Vin(i+colx,j+coly,bi,bj)
                0161              ENDDO
                0162             ENDDO
                0163            ENDDO
                0164           ENDDO
                0165          ENDDO
                0166         ENDDO
                0167 
                0168       _EXCH_XY_RL( Au_SI, myThid )
                0169       _EXCH_XY_RL( Av_SI, myThid )
                0170 
                0171       DO bj = myByLo(myThid), myByHi(myThid)
                0172        DO bi = myBxLo(myThid), myBxHi(myThid)
                0173         DO j=1-OLy,sNy+OLy
                0174          DO i=1-OLx,sNx+OLx
                0175           Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
                0176      &     Au_SI(i,j,bi,bj)
                0177           Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
                0178      &     Av_SI(i,j,bi,bj)
                0179          ENDDO
                0180         ENDDO
                0181         dot_p1_tile(bi,bj) = 0. _d 0
07e785229e dngo*0182         dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0183        ENDDO
                0184       ENDDO
                0185 
                0186       DO bj = myByLo(myThid), myByHi(myThid)
                0187        DO bi = myBxLo(myThid), myBxHi(myThid)
                0188         DO j=1,sNy
                0189          DO i=1,sNx
                0190           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0191      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
                0192           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0193      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
                0194          ENDDO
                0195         ENDDO
                0196        ENDDO
                0197       ENDDO
                0198 
07e785229e dngo*0199       CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
5ca83cd8f7 Dani*0200       resid_0 = sqrt(dot_p1)
                0201 
223e7922c4 Dani*0202       DO bj = myByLo(myThid), myByHi(myThid)
                0203        DO bi = myBxLo(myThid), myBxHi(myThid)
                0204 
                0205        WRITE(msgBuf,'(A,I1,I1,E14.7)') 'CONJ GRAD INIT RESID LOCAL, ',
                0206      &         bi,bj, dot_p1_tile(bi,bj)
                0207        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0208      &                     SQUEEZE_RIGHT , 1)
                0209 
                0210        enddo
                0211       enddo
                0212 
5ca83cd8f7 Dani*0213       WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
                0214      &         resid_0
                0215        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0216      &                     SQUEEZE_RIGHT , 1)
                0217 
                0218 C    CCCCCCCCCCCCCCCCCCCC
                0219 
                0220       DO bj = myByLo(myThid), myByHi(myThid)
                0221        DO bi = myBxLo(myThid), myBxHi(myThid)
                0222         DO j=1-OLy,sNy+OLy
                0223          DO i=1-OLx,sNx+OLx
                0224           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0225      &      Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
                0226           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0227      &      Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
                0228          ENDDO
                0229         ENDDO
                0230        ENDDO
                0231       ENDDO
                0232 
                0233       cg_halo = min(OLx-1,OLy-1)
                0234       conv_flag = 0
                0235 
                0236       DO bj = myByLo(myThid), myByHi(myThid)
                0237        DO bi = myBxLo(myThid), myBxHi(myThid)
                0238         DO j=1-OLy,sNy+OLy
                0239          DO i=1-OLx,sNx+OLx
07e785229e dngo*0240           Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
5ca83cd8f7 Dani*0241           Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
                0242          ENDDO
                0243         ENDDO
                0244        ENDDO
                0245       ENDDO
                0246 
                0247       resid = resid_0
                0248       iters = 0
07e785229e dngo*0249 
5ca83cd8f7 Dani*0250 c  !!!!!!!!!!!!!!!!!!
                0251 c  !!              !!
                0252 c  !! MAIN CG LOOP !!
                0253 c  !!              !!
                0254 c  !!!!!!!!!!!!!!!!!!
                0255 
                0256 c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
                0257 
                0258       WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
                0259        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0260      &                     SQUEEZE_RIGHT , 1)
                0261 
07e785229e dngo*0262 c       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
5ca83cd8f7 Dani*0263 
d2cdb9260d Dani*0264       do iter = 1, maxIter
5ca83cd8f7 Dani*0265        if (resid .gt. tolerance*resid_0) then
                0266 
                0267 c      to avoid using "exit"
                0268        iters = iters + 1
                0269 
07e785229e dngo*0270        is = 1 - cg_halo
                0271        ie = sNx + cg_halo
                0272        js = 1 - cg_halo
5ca83cd8f7 Dani*0273        je = sNy + cg_halo
07e785229e dngo*0274 
5ca83cd8f7 Dani*0275        DO bj = myByLo(myThid), myByHi(myThid)
                0276         DO bi = myBxLo(myThid), myBxHi(myThid)
                0277          DO j=1-OLy,sNy+OLy
07e785229e dngo*0278           DO i=1-OLx,sNx+OLx
5ca83cd8f7 Dani*0279            Au_SI(i,j,bi,bj) = 0. _d 0
                0280            Av_SI(i,j,bi,bj) = 0. _d 0
                0281           ENDDO
                0282          ENDDO
                0283         ENDDO
                0284        ENDDO
                0285 
07e785229e dngo*0286 c        IF (STREAMICE_construct_matrix) THEN
5ca83cd8f7 Dani*0287 
07e785229e dngo*0288 c #ifdef STREAMICE_CONSTRUCT_MATRIX
5ca83cd8f7 Dani*0289 
                0290         DO bj = myByLo(myThid), myByHi(myThid)
                0291          DO bi = myBxLo(myThid), myBxHi(myThid)
                0292           DO j=js,je
                0293            DO i=is,ie
                0294             DO colx=-1,1
                0295              DO coly=-1,1
07e785229e dngo*0296               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
5ca83cd8f7 Dani*0297      &         A_uu(i,j,bi,bj,colx,coly)*
                0298      &         Du_SI(i+colx,j+coly,bi,bj)+
07e785229e dngo*0299      &         A_uv(i,j,bi,bj,colx,coly)*
5ca83cd8f7 Dani*0300      &         Dv_SI(i+colx,j+coly,bi,bj)
07e785229e dngo*0301               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
5ca83cd8f7 Dani*0302      &         A_vu(i,j,bi,bj,colx,coly)*
                0303      &         Du_SI(i+colx,j+coly,bi,bj)+
07e785229e dngo*0304      &         A_vv(i,j,bi,bj,colx,coly)*
5ca83cd8f7 Dani*0305      &         Dv_SI(i+colx,j+coly,bi,bj)
                0306              ENDDO
                0307             ENDDO
                0308            ENDDO
                0309           ENDDO
                0310          ENDDO
                0311         ENDDO
                0312 
07e785229e dngo*0313 c        else
                0314 c #else
                0315 c
                0316 c         CALL STREAMICE_CG_ACTION( myThid,
                0317 c      O     Au_SI,
                0318 c      O     Av_SI,
                0319 c      I     Du_SI,
                0320 c      I     Dv_SI,
                0321 c      I     is,ie,js,je)
                0322 c
                0323 c !        ENDIF
                0324 c
                0325 c #endif
5ca83cd8f7 Dani*0326 
                0327        DO bj = myByLo(myThid), myByHi(myThid)
                0328         DO bi = myBxLo(myThid), myBxHi(myThid)
                0329          dot_p1_tile(bi,bj) = 0. _d 0
07e785229e dngo*0330          dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0331         ENDDO
                0332        ENDDO
                0333 
                0334        DO bj = myByLo(myThid), myByHi(myThid)
                0335         DO bi = myBxLo(myThid), myBxHi(myThid)
                0336          DO j=1,sNy
                0337           DO i=1,sNx
                0338            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0339            dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
                0340      &        Ru_SI(i,j,bi,bj)
                0341             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
                0342      &        Au_SI(i,j,bi,bj)
                0343            ENDIF
                0344            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0345             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
                0346      &        Rv_SI(i,j,bi,bj)
                0347             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
                0348      &        Av_SI(i,j,bi,bj)
                0349            ENDIF
                0350           ENDDO
                0351          ENDDO
                0352         ENDDO
                0353        ENDDO
                0354 
                0355        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
07e785229e dngo*0356        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
5ca83cd8f7 Dani*0357        alpha_k = dot_p1/dot_p2
                0358 
                0359        DO bj = myByLo(myThid), myByHi(myThid)
                0360         DO bi = myBxLo(myThid), myBxHi(myThid)
                0361          DO j=1-OLy,sNy+OLy
                0362           DO i=1-OLx,sNx+OLx
                0363 
                0364            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0365             cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
                0366      &       alpha_k*Du_SI(i,j,bi,bj)
                0367             Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
                0368             Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
                0369             Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
                0370      &       alpha_k*Au_SI(i,j,bi,bj)
07e785229e dngo*0371             Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
5ca83cd8f7 Dani*0372      &       DIAGu_SI(i,j,bi,bj)
                0373            ENDIF
                0374 
                0375            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0376             cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
                0377      &       alpha_k*Dv_SI(i,j,bi,bj)
                0378             Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
                0379             Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
                0380             Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
                0381      &       alpha_k*Av_SI(i,j,bi,bj)
07e785229e dngo*0382             Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
5ca83cd8f7 Dani*0383      &       DIAGv_SI(i,j,bi,bj)
                0384 
                0385            ENDIF
                0386           ENDDO
                0387          ENDDO
                0388         ENDDO
                0389        ENDDO
                0390 
                0391        DO bj = myByLo(myThid), myByHi(myThid)
                0392         DO bi = myBxLo(myThid), myBxHi(myThid)
                0393          dot_p1_tile(bi,bj) = 0. _d 0
07e785229e dngo*0394          dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0395         ENDDO
                0396        ENDDO
                0397 
                0398        DO bj = myByLo(myThid), myByHi(myThid)
                0399         DO bi = myBxLo(myThid), myBxHi(myThid)
                0400          DO j=1,sNy
                0401           DO i=1,sNx
                0402 
                0403            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0404             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
                0405      &        Ru_SI(i,j,bi,bj)
                0406             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
                0407      &        Ru_old_SI(i,j,bi,bj)
                0408            ENDIF
                0409 
                0410            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0411             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
                0412      &        Rv_SI(i,j,bi,bj)
                0413             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
                0414      &        Rv_old_SI(i,j,bi,bj)
                0415            ENDIF
                0416 
                0417           ENDDO
                0418          ENDDO
                0419         ENDDO
                0420        ENDDO
                0421 
                0422        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
                0423        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
07e785229e dngo*0424 
5ca83cd8f7 Dani*0425        beta_k = dot_p1/dot_p2
                0426 
                0427        DO bj = myByLo(myThid), myByHi(myThid)
                0428         DO bi = myBxLo(myThid), myBxHi(myThid)
                0429          DO j=1-OLy,sNy+OLy
                0430           DO i=1-OLx,sNx+OLx
07e785229e dngo*0431            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
5ca83cd8f7 Dani*0432      &      Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
                0433      &      Zu_SI(i,j,bi,bj)
07e785229e dngo*0434            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
5ca83cd8f7 Dani*0435      &      Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
                0436      &      Zv_SI(i,j,bi,bj)
                0437           ENDDO
                0438          ENDDO
                0439         ENDDO
                0440        ENDDO
                0441 
                0442        DO bj = myByLo(myThid), myByHi(myThid)
                0443         DO bi = myBxLo(myThid), myBxHi(myThid)
                0444          dot_p1_tile(bi,bj) = 0. _d 0
                0445         ENDDO
                0446        ENDDO
                0447 
                0448        DO bj = myByLo(myThid), myByHi(myThid)
                0449         DO bi = myBxLo(myThid), myBxHi(myThid)
                0450          DO j=1,sNy
                0451           DO i=1,sNx
                0452            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0453      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
                0454            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0455      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
                0456           ENDDO
                0457          ENDDO
                0458         ENDDO
                0459        ENDDO
                0460 
07e785229e dngo*0461        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
5ca83cd8f7 Dani*0462        resid = sqrt(dot_p1)
                0463 
07e785229e dngo*0464 c        IF (iter .eq. 1) then
                0465 c         print *, alpha_k, beta_k, resid
                0466 c        ENDIF
5ca83cd8f7 Dani*0467 
                0468        cg_halo = cg_halo - 1
                0469 
                0470        if (cg_halo .eq. 0) then
                0471         cg_halo = min(OLx-1,OLy-1)
                0472         _EXCH_XY_RL( Du_SI, myThid )
                0473         _EXCH_XY_RL( Dv_SI, myThid )
                0474         _EXCH_XY_RL( Ru_SI, myThid )
                0475         _EXCH_XY_RL( Rv_SI, myThid )
                0476         _EXCH_XY_RL( cg_Uin, myThid )
                0477         _EXCH_XY_RL( cg_Vin, myThid )
                0478        endif
                0479 
                0480        endif
                0481       enddo ! end of CG loop
07e785229e dngo*0482 
5ca83cd8f7 Dani*0483 c     to avoid using "exit"
                0484 c     if iters has reached max_iters there is no convergence
07e785229e dngo*0485 
d2cdb9260d Dani*0486       IF (iters .lt. maxIter) THEN
5ca83cd8f7 Dani*0487        conv_flag = 1
                0488       ENDIF
95afe7199b Dani*0489       PRINT *, "GOT HERE CG ITERATIONS", iters
5ca83cd8f7 Dani*0490 
07e785229e dngo*0491 c       DO bj = myByLo(myThid), myByHi(myThid)
                0492 c        DO bi = myBxLo(myThid), myBxHi(myThid)
                0493 c         DO j=1-OLy,sNy+OLy
                0494 c          DO i=1-OLy,sNx+OLy
                0495 c           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
                0496 c      &     cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
                0497 c           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
                0498 c      &     cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
                0499 c          ENDDO
                0500 c         ENDDO
                0501 c        ENDDO
                0502 c       ENDDO
                0503 c
                0504 c       _EXCH_XY_RL( cg_Uin, myThid )
                0505 c       _EXCH_XY_RL( cg_Vin, myThid )
18a089944d Dani*0506 
                0507 #ifdef ALLOW_PETSC
                0508       endif   !if (streamice_use_petsc)
                0509 #endif
5ca83cd8f7 Dani*0510 
                0511 #else /* STREAMICE_SERIAL_TRISOLVE */
                0512 
af8b28c42b Dani*0513       iters = 0
                0514 
07e785229e dngo*0515       CALL STREAMICE_TRIDIAG_SOLVE(
5ca83cd8f7 Dani*0516      U                               cg_Uin,     ! x-velocities
                0517      U                               cg_Vin,
                0518      U                               cg_Bu,      ! force in x dir
                0519      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0520      I                               STREAMICE_umask,
                0521      I                               myThid )
                0522 
                0523 #endif
                0524 
                0525       CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
                0526 
                0527 #endif
                0528       RETURN
                0529       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!