C $Header: /u/gcmpack/MITgcm/model/src/pre_cg3d.F,v 1.6 2016/05/04 22:10:37 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: PRE_CG3D C !INTERFACE: SUBROUTINE PRE_CG3D( I oldFreeSurfTerm, I cg2d_x, U cg3d_b, I myTime, myIter, myThid ) C !DESCRIPTION: C Called from SOLVE_FOR_PRESSURE, before 3-D solver (cg3d): C Finish calculation of 3-D RHS after 2-D inversionis done. C !USES: IMPLICIT NONE C == Global variables #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" #include "FFIELDS.h" #include "DYNVARS.h" #ifdef ALLOW_NONHYDROSTATIC #include "NH_VARS.h" #endif C === Functions ==== c LOGICAL DIFFERENT_MULTIPLE c EXTERNAL DIFFERENT_MULTIPLE C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C oldFreeSurfTerm :: Treat free-surface term in the old way (no exactConserv) C cg2d_x :: Solution vector of the 2-D solver equation a.x=b C cg3d_b :: Right Hand side vector of the 3-D solver equation A.X=B C myTime :: Current time in simulation C myIter :: Current iteration number in simulation C myThid :: My Thread Id number LOGICAL oldFreeSurfTerm _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL myTime INTEGER myIter INTEGER myThid #ifdef ALLOW_NONHYDROSTATIC C !LOCAL VARIABLES: C == Local variables == C wSurfP2d :: surface vertical velocity after 2-D solver INTEGER i,j,k,bi,bj INTEGER ks, kp1 c CHARACTER*10 sufx c CHARACTER*(MAX_LEN_MBUF) msgBuf _RL locGamma, surfFac, tmpFac _RL wFacKm, wFacKp _RL uf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL wSurfP2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP c IF ( use3Dsolver ) THEN C-- Solve for a three-dimensional pressure term (NH or IGW or both ). C see CG3D.h for the interface to this routine. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C-- Calculate updated (after 2-D solver) vertical velocity at the surface IF ( oldFreeSurfTerm .OR. implicDiv2DFlow.EQ.zeroRL ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx wSurfP2d(i,j) = 0. _d 0 ENDDO ENDDO ELSE DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx wSurfP2d(i,j) = ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) ) & / ( implicDiv2DFlow*deltaTFreeSurf ) ENDDO ENDDO ENDIF C-- Add EmPmR contribution to top level cg3d_b or to wSurfP2d: C (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT) IF ( useRealFreshWaterFlux.AND.fluidIsWater ) THEN IF ( oldFreeSurfTerm .OR. usingPCoords ) THEN tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow/deltaTMom ks = 1 IF ( usingPCoords ) ks = Nr DO j=1,sNy DO i=1,sNx cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj) & + tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj) & *maskInC(i,j,bi,bj) ENDDO ENDDO ELSE DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx wSurfP2d(i,j) = wSurfP2d(i,j) & + EmPmR(i,j,bi,bj)*mass2rUnit & *maskInC(i,j,bi,bj) ENDDO ENDDO ENDIF ENDIF C-- Update or Add free-surface contribution to cg3d_b: surfFac = 0. IF ( selectNHfreeSurf.GE.1 ) THEN tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom DO j=1,sNy DO i=1,sNx locGamma = drC(1)*recip_Bo(i,j,bi,bj) & /( deltaTMom*deltaTFreeSurf & *implicitNHPress*implicDiv2DFlow ) ks = 1 c ks = kSurfC(i,j,bi,bj) c IF ( ks.LE.Nr ) THEN cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj) & + tmpFac*( wSurfP2d(i,j) & + locGamma*wVel(i,j,ks,bi,bj) ) & /( 1. _d 0 + locGamma ) & *_rA(i,j,bi,bj)*deepFac2F(ks) c ENDIF C- Save wSurfP2d (used in POST_CG3D) into dPhiNH : dPhiNH(i,j,bi,bj) = wSurfP2d(i,j) ENDDO ENDDO ELSEIF ( .NOT.oldFreeSurfTerm ) THEN tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) IF ( ks.LE.Nr ) THEN cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj) & + tmpFac*wSurfP2d(i,j) & *_rA(i,j,bi,bj)*deepFac2F(ks) ENDIF ENDDO ENDDO ELSEIF ( uniformFreeSurfLev ) THEN C- Z coordinate: assume surface @ level k=1 surfFac = freeSurfFac*deepFac2F(1) ELSE C- Other than Z coordinate: no assumption on surface level index DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) IF ( ks.LE.Nr ) THEN cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj) & +freeSurfFac*etaN(i,j,bi,bj)/deltaTFreeSurf & *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom ENDIF ENDDO ENDDO ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x C 2) add vertical velocity contribution. DO j=1,sNy+1 DO i=1,sNx+1 uf(i,j) = -_recip_dxC(i,j,bi,bj) & * implicSurfPress*implicDiv2DFlow & *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj)) #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) #endif vf(i,j) = -_recip_dyC(i,j,bi,bj) & * implicSurfPress*implicDiv2DFlow & *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj)) #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) #endif ENDDO ENDDO C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to C uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is C always recomputed from continuity eq (like eta when exactConserv=T) k=1 kp1 = MIN(k+1,Nr) wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1) IF (k.GE.Nr) wFacKp = 0. DO j=1,sNy DO i=1,sNx cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) & +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j) & -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j) & +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1) & -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j ) & +( surfFac*etaN(i,j,bi,bj)/deltaTFreeSurf & -wVel(i,j,kp1,bi,bj)*wFacKp & )*_rA(i,j,bi,bj)/deltaTMom ENDDO ENDDO DO k=2,Nr kp1 = MIN(k+1,Nr) C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ; C both appear in wVel term, but at 2 different levels wFacKm = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k ) wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1) IF (k.GE.Nr) wFacKp = 0. DO j=1,sNy DO i=1,sNx cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) & +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j) & -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j) & +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1) & -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j ) & +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj) & -wVel(i,j,kp1,bi,bj)*wFacKp & )*_rA(i,j,bi,bj)/deltaTMom ENDDO ENDDO ENDDO #ifdef ALLOW_OBCS C- Note: solver matrix is trivial outside OB region (main diagonal only) C => no real need to reset RHS (=cg3d_b) & cg3d_x, except that: C a) normalisation is fct of Max(RHS), which can be large ouside OB region C (would be different if we were solving for increment of phi_nh C instead of directly for phi_nh). C => need to reset RHS to ensure that interior solution does not depend C on ouside OB region. C b) provide directly the trivial solution cg3d_x == 0 for outside OB region C (=> no residual => no effect on solver convergence and interior solution) IF (useOBCS) THEN DO k=1,Nr DO j=1,sNy DO i=1,sNx cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) & *maskInC(i,j,bi,bj) phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj) & *maskInC(i,j,bi,bj) ENDDO ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS */ C- end bi,bj loops ENDDO ENDDO c ENDIF #endif /* ALLOW_NONHYDROSTATIC */ RETURN END