C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/bbl/code/Attic/mypackage_calc_rhs.F,v 1.1 2010/11/18 04:00:05 dimitri Exp $ C $Name: $ #include "BBL_OPTIONS.h" CBOP C !ROUTINE: BBL_CALC_RHS C !INTERFACE: SUBROUTINE MYPACKAGE_CALC_RHS( I bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: C Calculate tendency of tracers due to bottom boundary layer. C Scheme is currently coded for dTtracerLev(k) .EQ. deltaT. C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #include "BBL.h" C !INPUT PARAMETERS: C bi,bj :: Tile indices C myTime :: Current time in simulation C myIter :: Current time-step number C myThid :: my Thread Id number INTEGER bi, bj _RL myTime INTEGER myIter, myThid C !OUTPUT PARAMETERS: C !LOCAL VARIABLES: C i,j :: Loop indices C kBot :: k index of bottommost wet grid box C resThk :: thickness of bottommost wet grid box minus bbl_eta C resTheta :: temperature of this residual volume C resSalt :: salinity of this residual volume C deltaRho :: density change C deltaDpt :: depth change C bbl_tend :: temporary variable for tendency terms C sloc :: salinity of bottommost wet grid box C tloc :: temperature of bottommost wet grid box C rholoc :: Potential density of bottommost wet grid box C bbl_rho :: Potential density of bottom boundary layer C fZon :: Zonal flux C fMer :: Meridional flux INTEGER i, j, kBot _RL resThk, resTheta, resSalt _RL deltaRho, deltaDpt, bbl_tend _RL sloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) _RL tloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) _RL rholoc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) _RL bbl_rho ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) _RL fZon ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) _RL fMer ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) CEOP C Initialize tendency terms and C make local copy of bottomost temperature and salinity DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx bbl_TendTheta(i,j,bi,bj) = 0. _d 0 bbl_TendSalt (i,j,bi,bj) = 0. _d 0 kBot = max(1,kLowC(i,j,bi,bj)) tLoc(i,j) = theta(i,j,kBot,bi,bj) sLoc(i,j) = salt (i,j,kBot,bi,bj) ENDDO ENDDO C Calculate potential density of bottommost wet grid box CALL FIND_RHO_2D( I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, I tLoc, sLoc, O rhoLoc, I 1, bi, bj, myThid ) C Calculate potential density of bbl CALL FIND_RHO_2D( I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, I bbl_theta(1-OLx,1-OLy,bi,bj), bbl_salt(1-OLx,1-OLy,bi,bj), O bbl_rho, I 1, bi, bj, myThid ) C==== Compute and apply vertical exchange between BBL and C residual volume of botommost wet grid box. C This operation does not change total tracer quantity C in botommost wet grid box. DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx kBot = kLowC(i,j,bi,bj) IF ( kBot .GT. 0 ) THEN resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj) C If bbl occupies the complete bottom model grid box or C if model density is higher than BBL then mix instantly. IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.bbl_rho(i,j)) ) THEN bbl_theta(i,j,bi,bj) = tLoc(i,j) bbl_salt (i,j,bi,bj) = sLoc(i,j) C If model density is lower than BBL, slowly diffuse upward. ELSE resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) - & (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk resSalt = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) - & (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) + & deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) + & deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR ENDIF ENDIF ENDDO ENDDO C==== Compute zonal bbl exchange. DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx-1 IF ((kLowC(i,j,bi,bj).GT.0).AND.(kLowC(i+1,j,bi,bj).GT.0)) THEN deltaRho = bbl_rho(i+1,j) - bbl_rho(i,j) deltaDpt = R_low(i ,j,bi,bj) + bbl_eta(i ,j,bi,bj) - & R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj) C If heavy BBL water is higher than light BBL water, C exchange properties laterally. IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) + & ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) / & bbl_RelaxH bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) - & ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) * & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / & ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH ) bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) + & ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) / & bbl_RelaxH bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) - & ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) * & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / & ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH ) ENDIF ENDIF ENDDO ENDDO C==== Compute meridional bbl exchange. DO j=1-Oly,sNy+Oly-1 DO i=1-Olx,sNx+Olx IF ((kLowC(i,j,bi,bj).GT.0).AND.(kLowC(i,j+1,bi,bj).GT.0)) THEN deltaRho = bbl_rho(i,j+1) - bbl_rho(i,j) deltaDpt = R_low(i,j ,bi,bj) + bbl_eta(i,j ,bi,bj) - & R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj) C If heavy BBL water is higher than light BBL water, C exchange properties laterally. IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) + & ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) / & bbl_RelaxH bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) - & ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) * & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / & ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) / & bbl_RelaxH bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) + & ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) / & bbl_RelaxH bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) - & ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) * & ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / & ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH ) ENDIF ENDIF ENDDO ENDDO C==== Apply lateral BBL exchange then scale tendency term C for botommost wet grid box. DO j=1-Oly,sNy+Oly-1 DO i=1-Olx,sNx+Olx-1 kBot = kLowC(i,j,bi,bj) IF ( kBot .GT. 0 ) THEN bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) + & deltaT * bbl_TendTheta(i,j,bi,bj) bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) + & deltaT * bbl_TendSalt (i,j,bi,bj) bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) * & bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot)) bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) * & bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot)) ENDIF ENDDO ENDDO CALL EXCH_XY_RL( bbl_theta, myThid ) CALL EXCH_XY_RL( bbl_salt , myThid ) RETURN END