Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:09 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
be47244872 Jean*0002 #include "CPP_OPTIONS.h"
6d54cf9ca1 Ed H*0003 
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: UPDATE_CG2D
                0006 C     !INTERFACE:
be47244872 Jean*0007       SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid )
9366854e02 Chri*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
4606c28752 Jean*0010 C     | SUBROUTINE UPDATE_CG2D
                0011 C     | o Update 2d conjugate gradient solver operators
                0012 C     |   account for Free-Surf effect on total column thickness
9366854e02 Chri*0013 C     *==========================================================*
4606c28752 Jean*0014 C     | This routine is based on INI_CG2D, and simplified. It is
a8260bbda9 Jean*0015 C     | used when the non-linear free surface mode is activated
                0016 C     | or when bottom depth is part of the control vector.
9366854e02 Chri*0017 C     *==========================================================*
                0018 C     \ev
be47244872 Jean*0019 
9366854e02 Chri*0020 C     !USES:
                0021       IMPLICIT NONE
be47244872 Jean*0022 C     === Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "SURFACE.h"
                0028 #include "CG2D.h"
855d57fc61 Jean*0029 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
                0030 # include "DYNVARS.h"
                0031 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
be47244872 Jean*0032 
9366854e02 Chri*0033 C     !INPUT/OUTPUT PARAMETERS:
be47244872 Jean*0034 C     === Routine arguments ===
caa9902797 Jean*0035 C     myTime :: Current time in simulation
                0036 C     myIter :: Current iteration number in simulation
                0037 C     myThid :: Thread number for this instance of the routine.
be47244872 Jean*0038       _RL myTime
                0039       INTEGER myIter
                0040       INTEGER myThid
9366854e02 Chri*0041 
                0042 C     !LOCAL VARIABLES:
be47244872 Jean*0043 C-- Note : compared to "INI_CG2D", no needs to compute again
a8260bbda9 Jean*0044 C   the solver normalisation factor or the solver tolerance
be47244872 Jean*0045 C     === Local variables ===
4606c28752 Jean*0046 C     bi,bj  :: tile indices
caa9902797 Jean*0047 C     i,j,k  :: Loop counters
4606c28752 Jean*0048 C     faceArea :: Temporary used to hold cell face areas.
be47244872 Jean*0049       INTEGER bi, bj
caa9902797 Jean*0050       INTEGER i, j, k, ks
be47244872 Jean*0051       _RL     faceArea
29083c3b2a Jean*0052       _RL     pW_tmp, pS_tmp
fe587f155a Jean*0053       LOGICAL updatePreCond
9366854e02 Chri*0054 CEOP
be47244872 Jean*0055 
fe587f155a Jean*0056 C--   Decide when to update cg2d Preconditioner :
                0057       IF ( cg2dPreCondFreq.EQ.0 ) THEN
                0058         updatePreCond = .FALSE.
                0059       ELSE
                0060         updatePreCond = ( myIter.EQ.nIter0 )
                0061         IF ( MOD(myIter,cg2dPreCondFreq).EQ.0 ) updatePreCond=.TRUE.
                0062       ENDIF
                0063 
be47244872 Jean*0064 C--   Initialise laplace operator
                0065 C     aW2d: integral in Z Ax/dX
                0066 C     aS2d: integral in Z Ay/dY
                0067       DO bj=myByLo(myThid),myByHi(myThid)
                0068        DO bi=myBxLo(myThid),myBxHi(myThid)
b169b06b17 Jean*0069         DO j=1-OLy,sNy+OLy
                0070          DO i=1-OLx,sNx+OLx
caa9902797 Jean*0071           aW2d(i,j,bi,bj) = 0. _d 0
                0072           aS2d(i,j,bi,bj) = 0. _d 0
b169b06b17 Jean*0073 #ifdef ALLOW_AUTODIFF
                0074           aC2d(i,j,bi,bj) = 0. _d 0
                0075 #endif /* ALLOW_AUTODIFF */
be47244872 Jean*0076          ENDDO
                0077         ENDDO
855d57fc61 Jean*0078 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
                0079         IF ( selectImplicitDrag.EQ.2 ) THEN
                0080          DO k=1,Nr
                0081           DO j=1,sNy+1
                0082            DO i=1,sNx+1
                0083             faceArea = _dyG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
                0084      &                *_hFacW(i,j,k,bi,bj)
                0085             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
                0086      &               + faceArea*dU_psFacX(i,j,k,bi,bj)
                0087      &                         *recip_dxC(i,j,bi,bj)
                0088             faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k)
                0089      &                *_hFacS(i,j,k,bi,bj)
                0090             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
                0091      &               + faceArea*dV_psFacY(i,j,k,bi,bj)
                0092      &                         *recip_dyC(i,j,bi,bj)
                0093            ENDDO
                0094           ENDDO
                0095          ENDDO
                0096         ELSE
                0097 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
                0098          DO k=1,Nr
                0099           DO j=1,sNy+1
                0100            DO i=1,sNx+1
4606c28752 Jean*0101 C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
855d57fc61 Jean*0102             faceArea = _dyG(i,j,bi,bj)*drF(k)
                0103      &                *_hFacW(i,j,k,bi,bj)
                0104             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
                0105      &               + faceArea*recip_dxC(i,j,bi,bj)
                0106             faceArea = _dxG(i,j,bi,bj)*drF(k)
                0107      &                *_hFacS(i,j,k,bi,bj)
                0108             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
                0109      &               + faceArea*recip_dyC(i,j,bi,bj)
                0110            ENDDO
be47244872 Jean*0111           ENDDO
                0112          ENDDO
855d57fc61 Jean*0113 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
                0114         ENDIF
                0115 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
caa9902797 Jean*0116         DO j=1,sNy+1
                0117          DO i=1,sNx+1
                0118           aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
                0119      &                     *implicSurfPress*implicDiv2DFlow
be47244872 Jean*0120 #ifdef ALLOW_OBCS
caa9902797 Jean*0121      &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
be47244872 Jean*0122 #endif
caa9902797 Jean*0123           aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
be47244872 Jean*0124      &                     *implicSurfPress*implicDiv2DFlow
caa9902797 Jean*0125 #ifdef ALLOW_OBCS
                0126      &                  *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0127 #endif
be47244872 Jean*0128          ENDDO
                0129         ENDDO
4606c28752 Jean*0130 C--   compute matrix main diagonal :
                0131         IF ( deepAtmosphere ) THEN
caa9902797 Jean*0132          DO j=1,sNy
                0133           DO i=1,sNx
a8260bbda9 Jean*0134            ks = kSurfC(i,j,bi,bj)
caa9902797 Jean*0135            aC2d(i,j,bi,bj) = -(
                0136      &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
                0137      &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
                0138      &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
855d57fc61 Jean*0139      &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
4606c28752 Jean*0140      &                        )
                0141           ENDDO
29083c3b2a Jean*0142          ENDDO
4606c28752 Jean*0143         ELSE
caa9902797 Jean*0144          DO j=1,sNy
                0145           DO i=1,sNx
                0146            aC2d(i,j,bi,bj) = -(
                0147      &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
                0148      &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
                0149      &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)
855d57fc61 Jean*0150      &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
4606c28752 Jean*0151      &                        )
                0152           ENDDO
fe587f155a Jean*0153          ENDDO
                0154         ENDIF
                0155 C-    end bi,bj loops
be47244872 Jean*0156        ENDDO
                0157       ENDDO
4606c28752 Jean*0158 
fe587f155a Jean*0159       IF ( updatePreCond ) THEN
be47244872 Jean*0160 C--   Update overlap regions
4606c28752 Jean*0161       CALL EXCH_XY_RS(aC2d, myThid)
be47244872 Jean*0162 
                0163 C--   Initialise preconditioner
                0164       DO bj=myByLo(myThid),myByHi(myThid)
                0165        DO bi=myBxLo(myThid),myBxHi(myThid)
caa9902797 Jean*0166         DO j=1,sNy+1
                0167          DO i=1,sNx+1
                0168           IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
                0169             pC(i,j,bi,bj) = 1. _d 0
be47244872 Jean*0170           ELSE
caa9902797 Jean*0171            pC(i,j,bi,bj) =  1. _d 0 / aC2d(i,j,bi,bj)
be47244872 Jean*0172           ENDIF
caa9902797 Jean*0173           pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
29083c3b2a Jean*0174           IF ( pW_tmp .EQ. 0. ) THEN
caa9902797 Jean*0175            pW(i,j,bi,bj) = 0.
be47244872 Jean*0176           ELSE
caa9902797 Jean*0177            pW(i,j,bi,bj) =
                0178      &     -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 )
be47244872 Jean*0179           ENDIF
caa9902797 Jean*0180           pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
29083c3b2a Jean*0181           IF ( pS_tmp .EQ. 0. ) THEN
caa9902797 Jean*0182            pS(i,j,bi,bj) = 0.
be47244872 Jean*0183           ELSE
caa9902797 Jean*0184            pS(i,j,bi,bj) =
                0185      &     -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 )
be47244872 Jean*0186           ENDIF
                0187          ENDDO
                0188         ENDDO
                0189        ENDDO
                0190       ENDDO
fe587f155a Jean*0191 C-    if update Preconditioner : end
                0192       ENDIF
be47244872 Jean*0193 
                0194       RETURN
                0195       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!