C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc.F,v 1.34 2012/11/15 15:55:42 dimitri Exp $ C $Name: $ #include "OBCS_OPTIONS.h" CBOP C !ROUTINE: OBCS_CALC C !INTERFACE: SUBROUTINE OBCS_CALC( futureTime, futureIter, & uVel, vVel, wVel, theta, salt, & myThid ) C !DESCRIPTION: C *==========================================================* C | SUBROUTINE OBCS_CALC C | o Calculate future boundary data at open boundaries C | at time = futureTime C *==========================================================* C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "OBCS_PARAMS.h" #include "OBCS_GRID.h" #include "OBCS_FIELDS.h" #ifdef ALLOW_PTRACERS # include "PTRACERS_SIZE.h" # include "PTRACERS_PARAMS.h" # include "PTRACERS_FIELDS.h" # include "OBCS_PTRACERS.h" #endif /* ALLOW_PTRACERS */ #ifdef ALLOW_NEST_CHILD # include "NEST_CHILD.h" #endif /* ALLOW_NEST_CHILD */ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == INTEGER futureIter _RL futureTime _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) INTEGER myThid #ifdef ALLOW_OBCS C !LOCAL VARIABLES: C bi, bj :: tile indices C i,j,k :: loop indices C I_obc, J_obc :: local index of open boundary C msgBuf :: Informational/error message buffer INTEGER bi, bj INTEGER i, j, k #ifdef ALLOW_PTRACERS INTEGER I_obc, J_obc CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER iTracer #endif /* ALLOW_PTRACERS */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid) #endif DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) #ifdef ALLOW_NEST_CHILD IF ( useNEST_CHILD ) THEN IF ( PASSI.LT.2 ) THEN CALL NEST_CHILD_RECV ( myThid ) ENDIF ENDIF #endif /* ALLOW_NEST_CHILD */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_OBCS_EAST C Eastern OB #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: East',myThid) #endif IF (useOrlanskiEast) THEN #ifdef ALLOW_ORLANSKI CALL ORLANSKI_EAST( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) #endif #ifdef ALLOW_NEST_CHILD ELSEIF ( useNEST_CHILD ) THEN DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN OBEu(j,k,bi,bj)= U_F1(j,k,2) OBEv(j,k,bi,bj)= V_F1(j,k,2) OBEt(j,k,bi,bj)= T_F1(j,k,2) OBEs(j,k,bi,bj)= S_F1(j,k,2) #ifdef NONLIN_FRSURF OBEeta(j,bi,bj)= ETA_F1(j,1,2) #endif ENDIF ENDDO ENDDO #endif /* ALLOW_NEST_CHILD */ ELSE DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN OBEu(j,k,bi,bj)=0. OBEv(j,k,bi,bj)=0. OBEt(j,k,bi,bj)=tRef(k) OBEs(j,k,bi,bj)=sRef(k) #ifdef ALLOW_NONHYDROSTATIC OBEw(j,k,bi,bj)=0. #endif #ifdef NONLIN_FRSURF OBEeta(j,bi,bj)=0. #endif ENDIF ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS_EAST */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_OBCS_WEST C Western OB #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: West',myThid) #endif IF (useOrlanskiWest) THEN #ifdef ALLOW_ORLANSKI CALL ORLANSKI_WEST( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) #endif #ifdef ALLOW_NEST_CHILD ELSEIF ( useNEST_CHILD ) THEN DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN OBWu(j,k,bi,bj)= U_F1(j,k,1) OBWv(j,k,bi,bj)= V_F1(j,k,1) OBWt(j,k,bi,bj)= T_F1(j,k,1) OBWs(j,k,bi,bj)= S_F1(j,k,1) #ifdef NONLIN_FRSURF OBWeta(j,bi,bj)= ETA_F1(j,1,1) #endif ENDIF ENDDO ENDDO #endif /* ALLOW_NEST_CHILD */ ELSE DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN OBWu(j,k,bi,bj)=0. OBWv(j,k,bi,bj)=0. OBWt(j,k,bi,bj)=tRef(k) OBWs(j,k,bi,bj)=sRef(k) #ifdef ALLOW_NONHYDROSTATIC OBWw(j,k,bi,bj)=0. #endif #ifdef NONLIN_FRSURF OBWeta(j,bi,bj)=0. #endif ENDIF ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS_WEST */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_OBCS_NORTH C Northern OB #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: North',myThid) #endif IF (useOrlanskiNorth) THEN #ifdef ALLOW_ORLANSKI CALL ORLANSKI_NORTH( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) #endif ELSE DO k=1,Nr DO i=1-OLx,sNx+OLx IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN OBNv(i,k,bi,bj)=0. OBNu(i,k,bi,bj)=0. OBNt(i,k,bi,bj)=tRef(k) OBNs(i,k,bi,bj)=sRef(k) #ifdef ALLOW_NONHYDROSTATIC OBNw(i,k,bi,bj)=0. #endif #ifdef NONLIN_FRSURF OBNeta(i,bi,bj)=0. #endif ENDIF ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS_NORTH */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_OBCS_SOUTH C Southern OB #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_MSG('OBCS_CALC: South',myThid) #endif IF (useOrlanskiSouth) THEN #ifdef ALLOW_ORLANSKI CALL ORLANSKI_SOUTH( & bi, bj, futureTime, & uVel, vVel, wVel, theta, salt, & myThid ) #endif ELSE DO k=1,Nr DO i=1-OLx,sNx+OLx IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN OBSu(i,k,bi,bj)=0. OBSv(i,k,bi,bj)=0. OBSt(i,k,bi,bj)=tRef(k) OBSs(i,k,bi,bj)=sRef(k) #ifdef ALLOW_NONHYDROSTATIC OBSw(i,k,bi,bj)=0. #endif #ifdef NONLIN_FRSURF OBSeta(i,bi,bj)=0. #endif ENDIF ENDDO ENDDO ENDIF #endif /* ALLOW_OBCS_SOUTH */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_PTRACERS IF ( usePTRACERS ) THEN C C Calculate some default open boundary conditions for passive tracers: C The default is a homogeneous v.Neumann conditions, that is, the C tracer gradient across the open boundary is nearly zero; C only nearly, because the boundary conditions are applied throughout C the time step during which the interior field does change; therefore C we have to use the values from the previous time step here. If you C really want exact v.Neumann conditions, you have to modify C obcs_apply_ptracer directly. C # ifdef ALLOW_OBCS_EAST C Eastern OB # ifdef ALLOW_DEBUG IF (debugMode) & CALL DEBUG_MSG('OBCS_CALC: East, pTracers',myThid) # endif IF (useOrlanskiEast) THEN WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: useOrlanskiEast Rad OBC with' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: pTracers not yet implemented' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R OBCS_CALC' ELSE DO iTracer=1,PTRACERS_numInUse DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN I_obc = OB_Ie(j,bi,bj) OBEptr(j,k,bi,bj,iTracer) = & pTracer(I_obc-1,j,k,bi,bj,iTracer) & *_maskW(I_obc,j,k,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDIF # endif /* ALLOW_OBCS_EAST */ C ------------------------------------------------------------------------------ # ifdef ALLOW_OBCS_WEST C Western OB # ifdef ALLOW_DEBUG IF (debugMode) & CALL DEBUG_MSG('OBCS_CALC: West, pTracers',myThid) # endif IF (useOrlanskiWest) THEN WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: useOrlanskiWest Rad OBC with' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: pTracers not yet implemented' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R OBCS_CALC' ELSE DO iTracer=1,PTRACERS_numInUse DO k=1,Nr DO j=1-OLy,sNy+OLy IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN I_obc = OB_Iw(j,bi,bj) OBWptr(j,k,bi,bj,iTracer) = & pTracer(I_obc+1,j,k,bi,bj,iTracer) & *_maskW(I_obc+1,j,k,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDIF # endif /* ALLOW_OBCS_WEST */ C ------------------------------------------------------------------------------ # ifdef ALLOW_OBCS_NORTH C Northern OB # ifdef ALLOW_DEBUG IF (debugMode) & CALL DEBUG_MSG('OBCS_CALC: North, pTracers',myThid) # endif IF (useOrlanskiNorth) THEN WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: useOrlanskiNorth Rad OBC with' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: pTracers not yet implemented' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R OBCS_CALC' ELSE DO iTracer=1,PTRACERS_numInUse DO k=1,Nr DO i=1-OLx,sNx+OLx IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN J_obc = OB_Jn(i,bi,bj) OBNptr(i,k,bi,bj,iTracer) = & pTracer(i,J_obc-1,k,bi,bj,iTracer) & *_maskS(i,J_obc,k,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDIF # endif /* ALLOW_OBCS_NORTH */ C ------------------------------------------------------------------------------ # ifdef ALLOW_OBCS_SOUTH C Southern OB # ifdef ALLOW_DEBUG IF (debugMode) & CALL DEBUG_MSG('OBCS_CALC: South, pTracers',myThid) #endif IF (useOrlanskiSouth) THEN WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: useOrlanskiSouth Rad OBC with' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') & 'OBCS_CALC: ERROR: pTracers not yet implemented' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R OBCS_CALC' ELSE DO iTracer=1,PTRACERS_numInUse DO k=1,Nr DO i=1-OLx,sNx+OLx IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN J_obc = OB_Js(i,bi,bj) OBSptr(i,k,bi,bj,iTracer) = & pTracer(i,J_obc+1,k,bi,bj,iTracer) & *_maskS(i,J_obc+1,k,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDIF # endif /* ALLOW_OBCS_SOUTH */ C end if (usePTracers) ENDIF #endif /* ALLOW_PTRACERS */ C-- end bi,bj loops. ENDDO ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_OBCS_PRESCRIBE IF (useOBCSprescribe) THEN C-- Calculate future values on open boundaries #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('OBCS_PRESCRIBE_READ',myThid) #endif CALL OBCS_PRESCRIBE_READ( futureTime, futureIter, myThid ) ENDIF #endif /* ALLOW_OBCS_PRESCRIBE */ C ------------------------------------------------------------------------------ #ifdef ALLOW_OBCS_STEVENS C The Stevens (1990) boundary conditions come after reading data from C files, because they use the data to compute a mix of simplified C Orlanski and prescribed boundary conditions IF (useStevensNorth.OR.useStevensSouth.OR. & useStevensEast.OR.useStevensWest) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('OBCS_CALC_STEVENS',myThid) #endif CALL OBCS_CALC_STEVENS( futureTime, futureIter, myThid ) ENDIF #endif /* ALLOW_OBCS_STEVENS */ #ifdef ALLOW_OBCS_BALANCE IF ( useOBCSbalance ) THEN CALL OBCS_BALANCE_FLOW( futureTime, futureIter, myThid ) ENDIF #endif /* ALLOW_OBCS_BALANCE */ #ifdef ALLOW_OBCS_TIDES IF ( useOBCStides ) THEN CALL OBCS_ADD_TIDES( futureTime, futureIter, myThid ) ENDIF #endif /* ALLOW_OBCS_TIDES */ #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid) #endif #endif /* ALLOW_OBCS */ RETURN END