Back to home page

MITgcm
 
 

    

File indexing completed on 2023-09-21 05:10:48 UTC

view on githubraw file Latest commit 96b00645 on 2023-09-20 15:15:14 UTC
5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 
                0005 CBOP
96b006450c dngo*0006       SUBROUTINE STREAMICE_CG_SOLVE_MATFREE (
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
96b006450c dngo*0011      I                               tolerance,
5ca83cd8f7 Dani*0012      O                               iters,
                0013      I                               myThid )
                0014 C     /============================================================\
96b006450c dngo*0015 C     | SUBROUTINE                                                 |
5ca83cd8f7 Dani*0016 C     | o                                                          |
                0017 C     |============================================================|
                0018 C     |                                                            |
                0019 C     \============================================================/
                0020       IMPLICIT NONE
                0021 
                0022 C     === Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "STREAMICE.h"
                0027 #include "STREAMICE_CG.h"
                0028 
                0029 C     !INPUT/OUTPUT ARGUMENTS
                0030 C     cg_Uin, cg_Vin - input and output velocities
                0031 C     cg_Bu, cg_Bv - driving stress
                0032       INTEGER myThid
                0033       INTEGER iters
                0034       _RL tolerance
                0035       _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0036       _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0037       _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038       _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0039 
                0040 C     LOCAL VARIABLES
                0041       INTEGER i, j, bi, bj, cg_halo, conv_flag
96b006450c dngo*0042       INTEGER iter, is, js, ie, je
                0043 c     INTEGER colx, coly
5ca83cd8f7 Dani*0044       _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
                0045       _RL dot_p1_tile (nSx,nSy)
                0046       _RL dot_p2_tile (nSx,nSy)
                0047       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0048 
                0049       iters = streamice_max_cg_iter
                0050 
                0051 #ifdef ALLOW_STREAMICE
                0052 
                0053       conv_flag = 0
                0054 
                0055       DO bj = myByLo(myThid), myByHi(myThid)
                0056        DO bi = myBxLo(myThid), myBxHi(myThid)
                0057         DO j=1,sNy
                0058          DO i=1,sNx
                0059           Zu_SI (i,j,bi,bj) = 0. _d 0
                0060           Zv_SI (i,j,bi,bj) = 0. _d 0
                0061           Ru_SI (i,j,bi,bj) = 0. _d 0
                0062           Rv_SI (i,j,bi,bj) = 0. _d 0
                0063           Au_SI (i,j,bi,bj) = 0. _d 0
                0064           Av_SI (i,j,bi,bj) = 0. _d 0
                0065           Du_SI (i,j,bi,bj) = 0. _d 0
                0066           Dv_SI (i,j,bi,bj) = 0. _d 0
                0067          ENDDO
                0068         ENDDO
                0069        ENDDO
                0070       ENDDO
96b006450c dngo*0071 
5ca83cd8f7 Dani*0072 C     FIND INITIAL RESIDUAL, and initialize r
                0073 
96b006450c dngo*0074 c #ifdef STREAMICE_CONSTRUCT_MATRIX
                0075 
                0076 c         DO bj = myByLo(myThid), myByHi(myThid)
                0077 c          DO bi = myBxLo(myThid), myBxHi(myThid)
                0078 c           DO j=js,je
                0079 c            DO i=is,ie
                0080 c             DO colx=-1,1
                0081 c              DO coly=-1,1
                0082 c               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
                0083 c      &         streamice_cg_A1(i,j,bi,bj,colx,coly)*
                0084 c      &         cg_Uin(i+colx,j+coly,bi,bj)+
                0085 c      &         streamice_cg_A2(i,j,bi,bj,colx,coly)*
                0086 c      &         cg_Vin(i+colx,j+coly,bi,bj)
                0087 c               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
                0088 c      &         streamice_cg_A3(i,j,bi,bj,colx,coly)*
                0089 c      &         cg_Uin(i+colx,j+coly,bi,bj)+
                0090 c      &         streamice_cg_A4(i,j,bi,bj,colx,coly)*
                0091 c      &         cg_Vin(i+colx,j+coly,bi,bj)
                0092 c              ENDDO
                0093 c             ENDDO
                0094 c            ENDDO
                0095 c           ENDDO
                0096 c          ENDDO
                0097 c         ENDDO
                0098 
                0099 c #else
                0100 
                0101       CALL STREAMICE_CG_ACTION( myThid,
                0102      O    Au_SI,
                0103      O    Av_SI,
                0104      I    cg_Uin,
                0105      I    cg_Vin,
5ca83cd8f7 Dani*0106      I    0, sNx+1, 0, sNy+1 )
                0107 
96b006450c dngo*0108 c #endif
5ca83cd8f7 Dani*0109 
                0110       _EXCH_XY_RL( Au_SI, myThid )
                0111       _EXCH_XY_RL( Av_SI, myThid )
                0112 
                0113       DO bj = myByLo(myThid), myByHi(myThid)
                0114        DO bi = myBxLo(myThid), myBxHi(myThid)
                0115         DO j=1-OLy,sNy+OLy
                0116          DO i=1-OLx,sNx+OLx
                0117           Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
                0118      &     Au_SI(i,j,bi,bj)
                0119           Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
                0120      &     Av_SI(i,j,bi,bj)
                0121          ENDDO
                0122         ENDDO
                0123         dot_p1_tile(bi,bj) = 0. _d 0
96b006450c dngo*0124         dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0125        ENDDO
                0126       ENDDO
                0127 
                0128       DO bj = myByLo(myThid), myByHi(myThid)
                0129        DO bi = myBxLo(myThid), myBxHi(myThid)
                0130         DO j=1,sNy
                0131          DO i=1,sNx
                0132           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0133      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
                0134           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0135      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
                0136          ENDDO
                0137         ENDDO
                0138        ENDDO
                0139       ENDDO
                0140 
96b006450c dngo*0141       CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
5ca83cd8f7 Dani*0142       resid_0 = sqrt(dot_p1)
                0143 
                0144 C    CCCCCCCCCCCCCCCCCCCC
                0145 
                0146       DO bj = myByLo(myThid), myByHi(myThid)
                0147        DO bi = myBxLo(myThid), myBxHi(myThid)
                0148         DO j=1-OLy,sNy+OLy
                0149          DO i=1-OLx,sNx+OLx
                0150           IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0151      &      Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
                0152           IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0153      &      Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
                0154          ENDDO
                0155         ENDDO
                0156        ENDDO
                0157       ENDDO
                0158 
                0159       cg_halo = min(OLx-1,OLy-1)
                0160       conv_flag = 0
                0161 
                0162       DO bj = myByLo(myThid), myByHi(myThid)
                0163        DO bi = myBxLo(myThid), myBxHi(myThid)
                0164         DO j=1-OLy,sNy+OLy
                0165          DO i=1-OLx,sNx+OLx
96b006450c dngo*0166           Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
5ca83cd8f7 Dani*0167           Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
                0168          ENDDO
                0169         ENDDO
                0170        ENDDO
                0171       ENDDO
                0172 
                0173       resid = resid_0
                0174       iters = 0
96b006450c dngo*0175 
5ca83cd8f7 Dani*0176 c  !!!!!!!!!!!!!!!!!!
                0177 c  !!              !!
                0178 c  !! MAIN CG LOOP !!
                0179 c  !!              !!
                0180 c  !!!!!!!!!!!!!!!!!!
                0181 
                0182 c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
                0183 
                0184        WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
                0185        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0186      &                     SQUEEZE_RIGHT , 1)
                0187 
96b006450c dngo*0188 c       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
5ca83cd8f7 Dani*0189 
                0190       do iter = 1, streamice_max_cg_iter
                0191        if (resid .gt. tolerance*resid_0) then
                0192 
                0193 c      to avoid using "exit"
                0194        iters = iters + 1
                0195 
96b006450c dngo*0196        is = 1 - cg_halo
                0197        ie = sNx + cg_halo
                0198        js = 1 - cg_halo
5ca83cd8f7 Dani*0199        je = sNy + cg_halo
96b006450c dngo*0200 
5ca83cd8f7 Dani*0201        DO bj = myByLo(myThid), myByHi(myThid)
                0202         DO bi = myBxLo(myThid), myBxHi(myThid)
                0203          DO j=1-OLy,sNy+OLy
96b006450c dngo*0204           DO i=1-OLx,sNx+OLx
5ca83cd8f7 Dani*0205            Au_SI(i,j,bi,bj) = 0. _d 0
                0206            Av_SI(i,j,bi,bj) = 0. _d 0
                0207           ENDDO
                0208          ENDDO
                0209         ENDDO
                0210        ENDDO
                0211 
96b006450c dngo*0212 c        IF (STREAMICE_construct_matrix) THEN
                0213 
                0214 c #ifdef STREAMICE_CONSTRUCT_MATRIX
                0215 c
                0216 c         DO bj = myByLo(myThid), myByHi(myThid)
                0217 c          DO bi = myBxLo(myThid), myBxHi(myThid)
                0218 c           DO j=js,je
                0219 c            DO i=is,ie
                0220 c             DO colx=-1,1
                0221 c              DO coly=-1,1
                0222 c               Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
                0223 c      &         streamice_cg_A1(i,j,bi,bj,colx,coly)*
                0224 c      &         Du_SI(i+colx,j+coly,bi,bj)+
                0225 c      &         streamice_cg_A2(i,j,bi,bj,colx,coly)*
                0226 c      &         Dv_SI(i+colx,j+coly,bi,bj)
                0227 c               Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
                0228 c      &         streamice_cg_A3(i,j,bi,bj,colx,coly)*
                0229 c      &         Du_SI(i+colx,j+coly,bi,bj)+
                0230 c      &         streamice_cg_A4(i,j,bi,bj,colx,coly)*
                0231 c      &         Dv_SI(i+colx,j+coly,bi,bj)
                0232 c              ENDDO
                0233 c             ENDDO
                0234 c            ENDDO
                0235 c           ENDDO
                0236 c          ENDDO
                0237 c         ENDDO
                0238 c
                0239 c !        else
                0240 c #else
                0241 
                0242         CALL STREAMICE_CG_ACTION( myThid,
                0243      O     Au_SI,
                0244      O     Av_SI,
                0245      I     Du_SI,
                0246      I     Dv_SI,
5ca83cd8f7 Dani*0247      I     is,ie,js,je)
                0248 
96b006450c dngo*0249 c        ENDIF
5ca83cd8f7 Dani*0250 
96b006450c dngo*0251 c #endif
5ca83cd8f7 Dani*0252 
                0253        DO bj = myByLo(myThid), myByHi(myThid)
                0254         DO bi = myBxLo(myThid), myBxHi(myThid)
                0255          dot_p1_tile(bi,bj) = 0. _d 0
96b006450c dngo*0256          dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0257         ENDDO
                0258        ENDDO
                0259 
                0260        DO bj = myByLo(myThid), myByHi(myThid)
                0261         DO bi = myBxLo(myThid), myBxHi(myThid)
                0262          DO j=1,sNy
                0263           DO i=1,sNx
                0264            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0265            dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
                0266      &        Ru_SI(i,j,bi,bj)
                0267             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
                0268      &        Au_SI(i,j,bi,bj)
                0269            ENDIF
                0270            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0271             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
                0272      &        Rv_SI(i,j,bi,bj)
                0273             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
                0274      &        Av_SI(i,j,bi,bj)
                0275            ENDIF
                0276           ENDDO
                0277          ENDDO
                0278         ENDDO
                0279        ENDDO
                0280 
                0281        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
96b006450c dngo*0282        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
5ca83cd8f7 Dani*0283        alpha_k = dot_p1/dot_p2
                0284 
                0285        DO bj = myByLo(myThid), myByHi(myThid)
                0286         DO bi = myBxLo(myThid), myBxHi(myThid)
                0287          DO j=1-OLy,sNy+OLy
                0288           DO i=1-OLx,sNx+OLx
                0289 
                0290            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0291             cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
                0292      &       alpha_k*Du_SI(i,j,bi,bj)
                0293             Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
                0294             Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
                0295             Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
                0296      &       alpha_k*Au_SI(i,j,bi,bj)
96b006450c dngo*0297             Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
5ca83cd8f7 Dani*0298      &       DIAGu_SI(i,j,bi,bj)
                0299            ENDIF
                0300 
                0301            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0302             cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
                0303      &       alpha_k*Dv_SI(i,j,bi,bj)
                0304             Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
                0305             Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
                0306             Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
                0307      &       alpha_k*Av_SI(i,j,bi,bj)
96b006450c dngo*0308             Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
5ca83cd8f7 Dani*0309      &       DIAGv_SI(i,j,bi,bj)
                0310 
                0311            ENDIF
                0312           ENDDO
                0313          ENDDO
                0314         ENDDO
                0315        ENDDO
                0316 
                0317        DO bj = myByLo(myThid), myByHi(myThid)
                0318         DO bi = myBxLo(myThid), myBxHi(myThid)
                0319          dot_p1_tile(bi,bj) = 0. _d 0
96b006450c dngo*0320          dot_p2_tile(bi,bj) = 0. _d 0
5ca83cd8f7 Dani*0321         ENDDO
                0322        ENDDO
                0323 
                0324        DO bj = myByLo(myThid), myByHi(myThid)
                0325         DO bi = myBxLo(myThid), myBxHi(myThid)
                0326          DO j=1,sNy
                0327           DO i=1,sNx
                0328 
                0329            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
                0330             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
                0331      &        Ru_SI(i,j,bi,bj)
                0332             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
                0333      &        Ru_old_SI(i,j,bi,bj)
                0334            ENDIF
                0335 
                0336            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
                0337             dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
                0338      &        Rv_SI(i,j,bi,bj)
                0339             dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
                0340      &        Rv_old_SI(i,j,bi,bj)
                0341            ENDIF
                0342 
                0343           ENDDO
                0344          ENDDO
                0345         ENDDO
                0346        ENDDO
                0347 
                0348        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
                0349        CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
96b006450c dngo*0350 
5ca83cd8f7 Dani*0351        beta_k = dot_p1/dot_p2
                0352 
                0353        DO bj = myByLo(myThid), myByHi(myThid)
                0354         DO bi = myBxLo(myThid), myBxHi(myThid)
                0355          DO j=1-OLy,sNy+OLy
                0356           DO i=1-OLx,sNx+OLx
96b006450c dngo*0357            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
5ca83cd8f7 Dani*0358      &      Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
                0359      &      Zu_SI(i,j,bi,bj)
96b006450c dngo*0360            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
5ca83cd8f7 Dani*0361      &      Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
                0362      &      Zv_SI(i,j,bi,bj)
                0363           ENDDO
                0364          ENDDO
                0365         ENDDO
                0366        ENDDO
                0367 
                0368        DO bj = myByLo(myThid), myByHi(myThid)
                0369         DO bi = myBxLo(myThid), myBxHi(myThid)
                0370          dot_p1_tile(bi,bj) = 0. _d 0
                0371         ENDDO
                0372        ENDDO
                0373 
                0374        DO bj = myByLo(myThid), myByHi(myThid)
                0375         DO bi = myBxLo(myThid), myBxHi(myThid)
                0376          DO j=1,sNy
                0377           DO i=1,sNx
                0378            IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
                0379      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
                0380            IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
                0381      &      dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
                0382           ENDDO
                0383          ENDDO
                0384         ENDDO
                0385        ENDDO
                0386 
96b006450c dngo*0387        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
5ca83cd8f7 Dani*0388        resid = sqrt(dot_p1)
                0389 
96b006450c dngo*0390 c        IF (iter .eq. 1) then
                0391 c         print *, alpha_k, beta_k, resid
                0392 c        ENDIF
5ca83cd8f7 Dani*0393 
                0394        cg_halo = cg_halo - 1
                0395 
                0396        if (cg_halo .eq. 0) then
                0397         cg_halo = min(OLx-1,OLy-1)
                0398         _EXCH_XY_RL( Du_SI, myThid )
                0399         _EXCH_XY_RL( Dv_SI, myThid )
                0400         _EXCH_XY_RL( Ru_SI, myThid )
                0401         _EXCH_XY_RL( Rv_SI, myThid )
                0402         _EXCH_XY_RL( cg_Uin, myThid )
                0403         _EXCH_XY_RL( cg_Vin, myThid )
                0404        endif
                0405 
                0406        endif
                0407       enddo ! end of CG loop
96b006450c dngo*0408 
5ca83cd8f7 Dani*0409 c     to avoid using "exit"
                0410 c     if iters has reached max_iters there is no convergence
96b006450c dngo*0411 
5ca83cd8f7 Dani*0412       IF (iters .lt. streamice_max_cg_iter) THEN
                0413        conv_flag = 1
                0414       ENDIF
                0415 
96b006450c dngo*0416 c       DO bj = myByLo(myThid), myByHi(myThid)
                0417 c        DO bi = myBxLo(myThid), myBxHi(myThid)
                0418 c         DO j=1-OLy,sNy+OLy
                0419 c          DO i=1-OLy,sNx+OLy
                0420 c           IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
                0421 c      &     cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
                0422 c           IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
                0423 c      &     cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
                0424 c          ENDDO
                0425 c         ENDDO
                0426 c        ENDDO
                0427 c       ENDDO
                0428 c
                0429 c       _EXCH_XY_RL( cg_Uin, myThid )
                0430 c       _EXCH_XY_RL( cg_Vin, myThid )
5ca83cd8f7 Dani*0431 
                0432 #endif
                0433       RETURN
                0434       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!