C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.28 2011/12/22 19:00:58 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INTEGR_CONTINUITY C !INTERFACE: SUBROUTINE INTEGR_CONTINUITY( I uFld, vFld, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE INTEGR_CONTINUITY C | o Integrate the continuity Eq : compute vertical velocity C | and free surface "r-anomaly" (etaN,etaH) to satisfy C | exactly the conservation of the Total Volume C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" #include "SURFACE.h" #include "FFIELDS.h" C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C uFld :: Zonal velocity ( m/s ) C vFld :: Meridional velocity ( m/s ) C myTime :: Current time in simulation C myIter :: Current iteration number in simulation C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) C !LOCAL VARIABLES: C Local variables in common block C Local variables C bi,bj :: tile index C i,j,k :: Loop counters C uTrans :: Volume transports ( uVel.xA ) C vTrans :: Volume transports ( vVel.yA ) C hDivFlow :: Div. Barotropic Flow [transport unit m3/s] INTEGER k,bi,bj #ifdef EXACT_CONSERV INTEGER i, j INTEGER ks _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL hDivFlow(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL facEmP, facMass #else /* EXACT_CONSERV */ # ifdef ALLOW_OBCS INTEGER i, j # endif #endif /* EXACT_CONSERV */ #ifndef ALLOW_ADDFLUID _RL addMass(1) #endif /* ndef ALLOW_ADDFLUID */ #if (defined NONLIN_FRSURF) !(defined DISABLE_RSTAR_CODE) _RL rStarDhDt(1-OLx:sNx+OLx,1-OLy:sNy+OLy) #else _RL rStarDhDt(1) #endif CEOP C-- Start bi,bj loops DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef EXACT_CONSERV IF (exactConserv) THEN facEmP = 0. IF ( fluidIsWater.AND.useRealFreshWaterFlux ) facEmP=mass2rUnit facMass = 0. IF ( selectAddFluid.GE.1 ) facMass = mass2rUnit C-- Compute the Divergence of The Barotropic Flow : C- Initialise DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx hDivFlow(i,j) = 0. _d 0 utrans(i,j) = 0. _d 0 vtrans(i,j) = 0. _d 0 ENDDO ENDDO DO k=1,Nr C- Calculate velocity field "volume transports" through tracer cell faces C anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport). DO j=1,sNy+1 DO i=1,sNx+1 uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj) & *deepFacC(k)*rhoFacC(k) & *drF(k)*_hFacW(i,j,k,bi,bj) vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj) & *deepFacC(k)*rhoFacC(k) & *drF(k)*_hFacS(i,j,k,bi,bj) ENDDO ENDDO C- Integrate vertically the Horizontal Divergence DO j=1,sNy DO i=1,sNx hDivFlow(i,j) = hDivFlow(i,j) & +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j) & +vTrans(i,j+1)-vTrans(i,j) ) #ifdef ALLOW_ADDFLUID & -facMass*addMass(i,j,k,bi,bj) #endif /* ALLOW_ADDFLUID */ ENDDO ENDDO C- End DO k=1,Nr ENDDO C------------------------------------ C note: deep-model not implemented for P-coordinate + realFreshWaterFlux ; C anelastic: always assumes that rhoFacF(1) = 1 IF ( myIter.EQ.nIter0 .AND. myIter.NE.0 & .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN C needs previous time-step value of E-P-R, that has not been loaded C and was not in old pickup-file ; try to use etaN & etaH instead. IF ( usePickupBeforeC54 ) THEN DO j=1,sNy DO i=1,sNx dEtaHdt(i,j,bi,bj) = (etaN(i,j,bi,bj)-etaH(i,j,bi,bj)) & / (implicDiv2Dflow*deltaTfreesurf) ENDDO ENDDO ENDIF DO j=1,sNy DO i=1,sNx PmEpR(i,j,bi,bj) = dEtaHdt(i,j,bi,bj) & + hDivFlow(i,j)*recip_rA(i,j,bi,bj) & *recip_deepFac2F(1) PmEpR(i,j,bi,bj) = PmEpR(i,j,bi,bj)*rUnit2mass ENDDO ENDDO ELSEIF ( myIter.EQ.nIter0 ) THEN DO j=1,sNy DO i=1,sNx ks = kSurfC(I,J,bi,bj) PmEpR(i,j,bi,bj) = 0. _d 0 dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj) & *recip_deepFac2F(ks) ENDDO ENDDO ELSE C-- Needs to get valid PmEpR (for T,S forcing) also in overlap regions C (e.g., if using KPP) => set over full index range DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj) ENDDO ENDDO DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) dEtaHdt(i,j,bi,bj) = -hDivFlow(i,j)*recip_rA(i,j,bi,bj) & *recip_deepFac2F(ks) & -facEmP*EmPmR(i,j,bi,bj) ENDDO ENDDO ENDIF C------------------------------------ #ifdef ALLOW_OBCS C-- reset dEtaHdt to zero outside the OB interior region IF ( useOBCS ) THEN DO j=1,sNy DO i=1,sNx dEtaHdt(i,j,bi,bj) = dEtaHdt(i,j,bi,bj)*maskInC(i,j,bi,bj) ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS */ C-- end if exactConserv block ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( exactConserv .AND. myIter.NE.nIter0 ) THEN C-- Update etaN at the end of the time step : C Incorporate the Implicit part of -Divergence(Barotropic_Flow) IF (implicDiv2Dflow.EQ. 0. _d 0) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx etaN(i,j,bi,bj) = etaH(i,j,bi,bj) ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx etaN(i,j,bi,bj) = etaH(i,j,bi,bj) & + implicDiv2Dflow*dEtaHdt(i,j,bi,bj)*deltaTfreesurf ENDDO ENDDO ENDIF #ifdef ALLOW_OBCS C-- Was added on Dec 30, 2004 (to fix unrealistic etaN ?), but no longer C needed with proper masking in solver (matrix+cg2d_b,x) and masking C of dEtaHdt above. etaN next to OB does not enter present algorithm but C leave it commented out in case we would implement an aternative scheme. c IF ( useOBCS ) CALL OBCS_APPLY_ETA( bi, bj, etaN, myThid ) #endif /* ALLOW_OBCS */ C-- end if exactConserv and not myIter=nIter0 block ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| # ifdef NONLIN_FRSURF IF (select_rStar .NE. 0) THEN # ifndef DISABLE_RSTAR_CODE C-- note: rStarDhDt is similar to rStarDhCDt from S/R CALC_R_STAR C except for deep-factor and rho factor. DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) rStarDhDt(i,j) = dEtaHdt(i,j,bi,bj) & *deepFac2F(ks)*rhoFacF(ks) & *recip_Rcol(i,j,bi,bj) ENDDO ENDDO # endif /* DISABLE_RSTAR_CODE */ ENDIF # endif /* NONLIN_FRSURF */ #endif /* EXACT_CONSERV */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| DO k=Nr,1,-1 C-- Integrate continuity vertically for vertical velocity CALL INTEGRATE_FOR_W( I bi, bj, k, uFld, vFld, I addMass, rStarDhDt, O wVel, I myThid ) #ifdef EXACT_CONSERV IF ( k.EQ.Nr .AND. myIter.NE.0 .AND. usingPCoords & .AND. fluidIsWater .AND. useRealFreshWaterFlux ) THEN DO j=1,sNy DO i=1,sNx wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj) & +mass2rUnit*PmEpR(i,j,bi,bj)*maskC(i,j,k,bi,bj) ENDDO ENDDO ENDIF #endif /* EXACT_CONSERV */ #ifdef ALLOW_OBCS C-- reset W to zero outside the OB interior region IF ( useOBCS ) THEN DO j=1,sNy DO i=1,sNx wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)*maskInC(i,j,bi,bj) ENDDO ENDDO ENDIF C-- Apply OBC to W (non-hydrostatic): IF ( useOBCS.AND.nonHydrostatic ) & CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid ) #endif /* ALLOW_OBCS */ C- End DO k=Nr,1,-1 ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- End bi,bj loops ENDDO ENDDO IF ( exactConserv .AND. myIter.NE.nIter0 & .AND. implicDiv2Dflow .NE. 0. _d 0 ) & _EXCH_XY_RL( etaN , myThid ) IF ( implicitIntGravWave .OR. myIter.EQ.nIter0 ) & _EXCH_XYZ_RL( wVel , myThid ) #ifdef EXACT_CONSERV C-- Update etaH (from etaN and dEtaHdt): IF (exactConserv) THEN c IF ( useRealFreshWaterFlux .AND. myIter.EQ.nIter0 ) c & _EXCH_XY_RL( PmEpR, myThid ) #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('UPDATE_ETAH',myThid) #endif CALL UPDATE_ETAH( myTime, myIter, myThid ) ENDIF #ifdef NONLIN_FRSURF # ifndef DISABLE_SIGMA_CODE IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN CALL EXCH_XY_RL( dEtaHdt, myThid ) CALL UPDATE_ETAWS( myTime, myIter, myThid ) ENDIF # endif /* DISABLE_SIGMA_CODE */ #endif /* NONLIN_FRSURF */ #endif /* EXACT_CONSERV */ RETURN END