--- MITgcm/pkg/mom_fluxform/mom_fluxform.F 2002/11/07 21:51:15 1.7 +++ MITgcm/pkg/mom_fluxform/mom_fluxform.F 2003/01/26 21:18:50 1.8 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.7 2002/11/07 21:51:15 adcroft Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.8 2003/01/26 21:18:50 jmc Exp $ C $Name: $ CBOI @@ -35,7 +35,7 @@ I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown, I phi_hyd,KappaRU,KappaRV, U fVerU, fVerV, - I myCurrentTime,myIter,myThid) + I myTime,myIter,myThid) C !DESCRIPTION: C Calculates all the horizontal accelerations except for the implicit surface @@ -63,7 +63,7 @@ C KappaRV :: vertical viscosity C fVerU :: vertical flux of U, 2 1/2 dim for pipe-lining C fVerV :: vertical flux of V, 2 1/2 dim for pipe-lining -C myCurrentTime :: current time +C myTime :: current time C myIter :: current time-step number C myThid :: thread number INTEGER bi,bj,iMin,iMax,jMin,jMax @@ -73,7 +73,7 @@ _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) - _RL myCurrentTime + _RL myTime INTEGER myIter INTEGER myThid @@ -119,6 +119,8 @@ _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C I,J,K - Loop counters C rVelMaskOverride - Factor for imposing special surface boundary conditions C ( set according to free-surface condition ). @@ -174,6 +176,8 @@ pF(i,j) = 0. fZon(i,j) = 0. fMer(i,j) = 0. + rTransU(i,j) = 0. + rTransV(i,j) = 0. ENDDO ENDDO @@ -250,6 +254,41 @@ CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid) +C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp) + IF (momAdvection.AND.k.EQ.1) THEN + +C- Calculate vertical transports above U & V points (West & South face): + CALL MOM_CALC_RTRANS( k, bi, bj, + O rTransU, rTransV, + I myTime, myIter, myThid) + +C- Free surface correction term (flux at k=1) + CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + fVerU(i,j,kUp) = af(i,j) + ENDDO + ENDDO + + CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + fVerV(i,j,kUp) = af(i,j) + ENDDO + ENDDO + +C--- endif momAdvection & k=1 + ENDIF + + +C--- Calculate vertical transports (at k+1) below U & V points : + IF (momAdvection) THEN + CALL MOM_CALC_RTRANS( k+1, bi, bj, + O rTransU, rTransV, + I myTime, myIter, myThid) + ENDIF + + C---- Zonal momentum equation starts here C Bi-harmonic term del^2 U -> v4F @@ -294,18 +333,9 @@ C-- Vertical flux (fVer is at upper face of "u" cell) -C-- Free surface correction term (flux at k=1) - IF (momAdvection.AND.k.EQ.1) THEN - CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,af,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - fVerU(i,j,kUp) = af(i,j) - ENDDO - ENDDO - ENDIF C Mean flow component of vertical flux (at k+1) -> aF IF (momAdvection) - & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,af,myThid) + & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,rTransU,af,myThid) C Eddy component of vertical flux (interior component only) -> vrF IF (momViscosity.AND..NOT.implicitViscosity) @@ -347,6 +377,27 @@ ENDDO ENDDO +#ifdef NONLIN_FRSURF +C-- account for 3.D divergence of the flow in rStar coordinate: + IF ( momAdvection .AND. select_rStar.GT.0 ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + & - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf + & *uVel(i,j,k,bi,bj) + ENDDO + ENDDO + ENDIF + IF ( momAdvection .AND. select_rStar.LT.0 ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + & - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj) + ENDDO + ENDDO + ENDIF +#endif /* NONLIN_FRSURF */ + C-- No-slip and drag BCs appear as body forces in cell abutting topography IF (momViscosity.AND.no_slip_sides) THEN C- No-slip BCs impose a drag at walls... @@ -371,7 +422,7 @@ IF (momForcing) & CALL EXTERNAL_FORCING_U( I iMin,iMax,jMin,jMax,bi,bj,k, - I myCurrentTime,myThid) + I myTime,myThid) C-- Metric terms for curvilinear grid systems IF (useNHMTerms) THEN @@ -444,18 +495,9 @@ C-- Vertical flux (fVer is at upper face of "v" cell) -C-- Free surface correction term (flux at k=1) - IF (momAdvection.AND.k.EQ.1) THEN - CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,af,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - fVerV(i,j,kUp) = af(i,j) - ENDDO - ENDDO - ENDIF C o Mean flow component of vertical flux IF (momAdvection) - & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,af,myThid) + & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,rTransV,af,myThid) C Eddy component of vertical flux (interior component only) -> vrF IF (momViscosity.AND..NOT.implicitViscosity) @@ -497,6 +539,27 @@ ENDDO ENDDO +#ifdef NONLIN_FRSURF +C-- account for 3.D divergence of the flow in rStar coordinate: + IF ( momAdvection .AND. select_rStar.GT.0 ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + & - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf + & *vVel(i,j,k,bi,bj) + ENDDO + ENDDO + ENDIF + IF ( momAdvection .AND. select_rStar.LT.0 ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + & - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj) + ENDDO + ENDDO + ENDIF +#endif /* NONLIN_FRSURF */ + C-- No-slip and drag BCs appear as body forces in cell abutting topography IF (momViscosity.AND.no_slip_sides) THEN C- No-slip BCs impose a drag at walls... @@ -521,7 +584,7 @@ IF (momForcing) & CALL EXTERNAL_FORCING_V( I iMin,iMax,jMin,jMax,bi,bj,k, - I myCurrentTime,myThid) + I myTime,myThid) C-- Metric terms for curvilinear grid systems IF (useNHMTerms) THEN