C $Header: /u/gcmpack/MITgcm/model/src/update_cg2d.F,v 1.11 2016/11/28 23:05:05 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: UPDATE_CG2D C !INTERFACE: SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE UPDATE_CG2D C | o Update 2d conjugate gradient solver operators C | account for Free-Surf effect on total column thickness C *==========================================================* C | This routine is based on INI_CG2D, and simplified. It is C | used when the non-linear free surface mode is activated C | or when bottom depth is part of the control vector. C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" #include "CG2D.h" #ifdef ALLOW_SOLVE4_PS_AND_DRAG # include "DYNVARS.h" #endif /* ALLOW_SOLVE4_PS_AND_DRAG */ C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Current time in simulation C myIter :: Current iteration number in simulation C myThid :: Thread number for this instance of the routine. _RL myTime INTEGER myIter INTEGER myThid C !LOCAL VARIABLES: C-- Note : compared to "INI_CG2D", no needs to compute again C the solver normalisation factor or the solver tolerance C === Local variables === C bi,bj :: tile indices C i,j,k :: Loop counters C faceArea :: Temporary used to hold cell face areas. INTEGER bi, bj INTEGER i, j, k, ks _RL faceArea _RL pW_tmp, pS_tmp LOGICAL updatePreCond CEOP C-- Decide when to update cg2d Preconditioner : IF ( cg2dPreCondFreq.EQ.0 ) THEN updatePreCond = .FALSE. ELSE updatePreCond = ( myIter.EQ.nIter0 ) IF ( MOD(myIter,cg2dPreCondFreq).EQ.0 ) updatePreCond=.TRUE. ENDIF C-- Initialise laplace operator C aW2d: integral in Z Ax/dX C aS2d: integral in Z Ay/dY DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx aW2d(i,j,bi,bj) = 0. _d 0 aS2d(i,j,bi,bj) = 0. _d 0 #ifdef ALLOW_AUTODIFF aC2d(i,j,bi,bj) = 0. _d 0 #endif /* ALLOW_AUTODIFF */ ENDDO ENDDO #ifdef ALLOW_SOLVE4_PS_AND_DRAG IF ( selectImplicitDrag.EQ.2 ) THEN DO k=1,Nr DO j=1,sNy+1 DO i=1,sNx+1 faceArea = _dyG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k) & *_hFacW(i,j,k,bi,bj) aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj) & + faceArea*dU_psFacX(i,j,k,bi,bj) & *recip_dxC(i,j,bi,bj) faceArea = _dxG(i,j,bi,bj)*drF(k)*deepFacC(k)*rhoFacC(k) & *_hFacS(i,j,k,bi,bj) aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj) & + faceArea*dV_psFacY(i,j,k,bi,bj) & *recip_dyC(i,j,bi,bj) ENDDO ENDDO ENDDO ELSE #endif /* ALLOW_SOLVE4_PS_AND_DRAG */ DO k=1,Nr DO j=1,sNy+1 DO i=1,sNx+1 C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect faceArea = _dyG(i,j,bi,bj)*drF(k) & *_hFacW(i,j,k,bi,bj) aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj) & + faceArea*recip_dxC(i,j,bi,bj) faceArea = _dxG(i,j,bi,bj)*drF(k) & *_hFacS(i,j,k,bi,bj) aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj) & + faceArea*recip_dyC(i,j,bi,bj) ENDDO ENDDO ENDDO #ifdef ALLOW_SOLVE4_PS_AND_DRAG ENDIF #endif /* ALLOW_SOLVE4_PS_AND_DRAG */ DO j=1,sNy+1 DO i=1,sNx+1 aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm & *implicSurfPress*implicDiv2DFlow #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) #endif aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm & *implicSurfPress*implicDiv2DFlow #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) #endif ENDDO ENDDO C-- compute matrix main diagonal : IF ( deepAtmosphere ) THEN DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) aC2d(i,j,bi,bj) = -( & aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj) & +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj) & +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks) & *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf & ) ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx aC2d(i,j,bi,bj) = -( & aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj) & +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj) & +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj) & *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf & ) ENDDO ENDDO ENDIF C- end bi,bj loops ENDDO ENDDO IF ( updatePreCond ) THEN C-- Update overlap regions CALL EXCH_XY_RS(aC2d, myThid) C-- Initialise preconditioner DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy+1 DO i=1,sNx+1 IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN pC(i,j,bi,bj) = 1. _d 0 ELSE pC(i,j,bi,bj) = 1. _d 0 / aC2d(i,j,bi,bj) ENDIF pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj) IF ( pW_tmp .EQ. 0. ) THEN pW(i,j,bi,bj) = 0. ELSE pW(i,j,bi,bj) = & -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 ) ENDIF pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj) IF ( pS_tmp .EQ. 0. ) THEN pS(i,j,bi,bj) = 0. ELSE pS(i,j,bi,bj) = & -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 ) ENDIF ENDDO ENDDO ENDDO ENDDO C- if update Preconditioner : end ENDIF RETURN END