C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/obcs/obcs_init_variables.F,v 1.44 2012/11/15 20:46:52 dimitri Exp $
C $Name:  $

#include "OBCS_OPTIONS.h"

CBOP
C     !ROUTINE: OBCS_INIT_VARIABLES
C     !INTERFACE:
      SUBROUTINE OBCS_INIT_VARIABLES( myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE OBCS_INIT_VARIABLES
C     | o Initialise OBCs variable data
C     *==========================================================*
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "OBCS_PARAMS.h"
c#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "OBCS_SEAICE.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#include "PTRACERS_FIELDS.h"
#include "OBCS_PTRACERS.h"
#endif /* ALLOW_PTRACERS */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid :: my Thread Id Number
      INTEGER myThid
CEOP

#ifdef ALLOW_OBCS

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER bi, bj
      INTEGER i, j, k
      INTEGER fp
#ifdef ALLOW_PTRACERS
      INTEGER iTracer
#endif /* ALLOW_PTRACERS */

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('OBCS_INIT_VARIABLES',myThid)
#endif

      fp = readBinaryPrec

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

#ifdef ALLOW_OBCS_PRESCRIBE
        OBCS_ldRec(bi,bj) = 0
#endif
        DO k=1,Nr
         DO i=1-Olx,sNx+Olx
#ifdef ALLOW_OBCS_NORTH
          OBNu(i,k,bi,bj)=0. _d 0
          OBNv(i,k,bi,bj)=0. _d 0
          OBNt(i,k,bi,bj)=0. _d 0
          OBNs(i,k,bi,bj)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBNu0(i,k,bi,bj)=0. _d 0
          OBNv0(i,k,bi,bj)=0. _d 0
          OBNt0(i,k,bi,bj)=0. _d 0
          OBNs0(i,k,bi,bj)=0. _d 0
          OBNu1(i,k,bi,bj)=0. _d 0
          OBNv1(i,k,bi,bj)=0. _d 0
          OBNt1(i,k,bi,bj)=0. _d 0
          OBNs1(i,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_STEVENS
          OBNvStevens(i,k,bi,bj)=0. _d 0
          OBNtStevens(i,k,bi,bj)=0. _d 0
          OBNsStevens(i,k,bi,bj)=0. _d 0
# endif /* ALLOW_OBCS_STEVENS */
#endif /* ALLOW_OBCS_NORTH */

#ifdef ALLOW_OBCS_SOUTH
          OBSu(i,k,bi,bj)=0. _d 0
          OBSv(i,k,bi,bj)=0. _d 0
          OBSt(i,k,bi,bj)=0. _d 0
          OBSs(i,k,bi,bj)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBSu0(i,k,bi,bj)=0. _d 0
          OBSv0(i,k,bi,bj)=0. _d 0
          OBSt0(i,k,bi,bj)=0. _d 0
          OBSs0(i,k,bi,bj)=0. _d 0
          OBSu1(i,k,bi,bj)=0. _d 0
          OBSv1(i,k,bi,bj)=0. _d 0
          OBSt1(i,k,bi,bj)=0. _d 0
          OBSs1(i,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_STEVENS
          OBSvStevens(i,k,bi,bj)=0. _d 0
          OBStStevens(i,k,bi,bj)=0. _d 0
          OBSsStevens(i,k,bi,bj)=0. _d 0
# endif /* ALLOW_OBCS_STEVENS */
#endif /* ALLOW_OBCS_SOUTH */
         ENDDO

         DO j=1-Oly,sNy+Oly
#ifdef ALLOW_OBCS_EAST
          OBEu(j,k,bi,bj)=0. _d 0
          OBEv(j,k,bi,bj)=0. _d 0
          OBEt(j,k,bi,bj)=0. _d 0
          OBEs(j,k,bi,bj)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBEu0(j,k,bi,bj)=0. _d 0
          OBEv0(j,k,bi,bj)=0. _d 0
          OBEt0(j,k,bi,bj)=0. _d 0
          OBEs0(j,k,bi,bj)=0. _d 0
          OBEu1(j,k,bi,bj)=0. _d 0
          OBEv1(j,k,bi,bj)=0. _d 0
          OBEt1(j,k,bi,bj)=0. _d 0
          OBEs1(j,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_STEVENS
          OBEuStevens(j,k,bi,bj)=0. _d 0
          OBEtStevens(j,k,bi,bj)=0. _d 0
          OBEsStevens(j,k,bi,bj)=0. _d 0
# endif /* ALLOW_OBCS_STEVENS */
#endif /* ALLOW_OBCS_EAST */

#ifdef ALLOW_OBCS_WEST
          OBWu(j,k,bi,bj)=0. _d 0
          OBWv(j,k,bi,bj)=0. _d 0
          OBWt(j,k,bi,bj)=0. _d 0
          OBWs(j,k,bi,bj)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBWu0(j,k,bi,bj)=0. _d 0
          OBWv0(j,k,bi,bj)=0. _d 0
          OBWt0(j,k,bi,bj)=0. _d 0
          OBWs0(j,k,bi,bj)=0. _d 0
          OBWu1(j,k,bi,bj)=0. _d 0
          OBWv1(j,k,bi,bj)=0. _d 0
          OBWt1(j,k,bi,bj)=0. _d 0
          OBWs1(j,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_STEVENS
          OBWuStevens(j,k,bi,bj)=0. _d 0
          OBWtStevens(j,k,bi,bj)=0. _d 0
          OBWsStevens(j,k,bi,bj)=0. _d 0
# endif /* ALLOW_OBCS_STEVENS */
#endif /* ALLOW_OBCS_WEST */
         ENDDO
        ENDDO

#ifdef ALLOW_NONHYDROSTATIC
        DO k=1,Nr
         DO i=1-Olx,sNx+Olx
          OBNw (i,k,bi,bj) = 0. _d 0
          OBSw (i,k,bi,bj) = 0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBNw0(i,k,bi,bj) = 0. _d 0
          OBSw0(i,k,bi,bj) = 0. _d 0
          OBNw1(i,k,bi,bj) = 0. _d 0
          OBSw1(i,k,bi,bj) = 0. _d 0
# endif
         ENDDO
         DO j=1-Oly,sNy+Oly
          OBEw (j,k,bi,bj) = 0. _d 0
          OBWw (j,k,bi,bj) = 0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBEw0(j,k,bi,bj) = 0. _d 0
          OBWw0(j,k,bi,bj) = 0. _d 0
          OBEw1(j,k,bi,bj) = 0. _d 0
          OBWw1(j,k,bi,bj) = 0. _d 0
# endif
         ENDDO
        ENDDO
#endif /* ALLOW_NONHYDROSTATIC */

#ifdef NONLIN_FRSURF
        DO i=1-Olx,sNx+Olx
          OBNeta (i,bi,bj) = 0. _d 0
          OBSeta (i,bi,bj) = 0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBNeta0(i,bi,bj) = 0. _d 0
          OBSeta0(i,bi,bj) = 0. _d 0
          OBNeta1(i,bi,bj) = 0. _d 0
          OBSeta1(i,bi,bj) = 0. _d 0
# endif
        ENDDO
        DO j=1-Oly,sNy+Oly
          OBEeta (j,bi,bj) = 0. _d 0
          OBWeta (j,bi,bj) = 0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
          OBEeta0(j,bi,bj) = 0. _d 0
          OBWeta0(j,bi,bj) = 0. _d 0
          OBEeta1(j,bi,bj) = 0. _d 0
          OBWeta1(j,bi,bj) = 0. _d 0
# endif
        ENDDO
#endif /* NONLIN_FRSURF */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_SEAICE
        DO i=1-Olx,sNx+Olx
#ifdef ALLOW_OBCS_NORTH
         OBNa (i,bi,bj)=0. _d 0
         OBNh (i,bi,bj)=0. _d 0
         OBNa0(i,bi,bj)=0. _d 0
         OBNh0(i,bi,bj)=0. _d 0
         OBNa1(i,bi,bj)=0. _d 0
         OBNh1(i,bi,bj)=0. _d 0
         OBNsl (i,bi,bj)=0. _d 0
         OBNsn (i,bi,bj)=0. _d 0
         OBNsl0(i,bi,bj)=0. _d 0
         OBNsn0(i,bi,bj)=0. _d 0
         OBNsl1(i,bi,bj)=0. _d 0
         OBNsn1(i,bi,bj)=0. _d 0
         OBNuice (i,bi,bj)=0. _d 0
         OBNvice (i,bi,bj)=0. _d 0
         OBNuice0(i,bi,bj)=0. _d 0
         OBNvice0(i,bi,bj)=0. _d 0
         OBNuice1(i,bi,bj)=0. _d 0
         OBNvice1(i,bi,bj)=0. _d 0
#endif /* ALLOW_OBCS_NORTH */
#ifdef ALLOW_OBCS_SOUTH
         OBSa (i,bi,bj)=0. _d 0
         OBSh (i,bi,bj)=0. _d 0
         OBSa0(i,bi,bj)=0. _d 0
         OBSh0(i,bi,bj)=0. _d 0
         OBSa1(i,bi,bj)=0. _d 0
         OBSh1(i,bi,bj)=0. _d 0
         OBSsl (i,bi,bj)=0. _d 0
         OBSsn (i,bi,bj)=0. _d 0
         OBSsl0(i,bi,bj)=0. _d 0
         OBSsn0(i,bi,bj)=0. _d 0
         OBSsl1(i,bi,bj)=0. _d 0
         OBSsn1(i,bi,bj)=0. _d 0
         OBSuice (i,bi,bj)=0. _d 0
         OBSvice (i,bi,bj)=0. _d 0
         OBSuice0(i,bi,bj)=0. _d 0
         OBSvice0(i,bi,bj)=0. _d 0
         OBSuice1(i,bi,bj)=0. _d 0
         OBSvice1(i,bi,bj)=0. _d 0
#endif /* ALLOW_OBCS_SOUTH */
        ENDDO
        DO j=1-Oly,sNy+Oly
#ifdef ALLOW_OBCS_EAST
         OBEa (j,bi,bj)=0. _d 0
         OBEh (j,bi,bj)=0. _d 0
         OBEa0(j,bi,bj)=0. _d 0
         OBEh0(j,bi,bj)=0. _d 0
         OBEa1(j,bi,bj)=0. _d 0
         OBEh1(j,bi,bj)=0. _d 0
         OBEsl (j,bi,bj)=0. _d 0
         OBEsn (j,bi,bj)=0. _d 0
         OBEsl0(j,bi,bj)=0. _d 0
         OBEsn0(j,bi,bj)=0. _d 0
         OBEsl1(j,bi,bj)=0. _d 0
         OBEsn1(j,bi,bj)=0. _d 0
         OBEuice (j,bi,bj)=0. _d 0
         OBEvice (j,bi,bj)=0. _d 0
         OBEuice0(j,bi,bj)=0. _d 0
         OBEvice0(j,bi,bj)=0. _d 0
         OBEuice1(j,bi,bj)=0. _d 0
         OBEvice1(j,bi,bj)=0. _d 0
#endif /* ALLOW_OBCS_EAST */
#ifdef ALLOW_OBCS_WEST
         OBWa (j,bi,bj)=0. _d 0
         OBWh (j,bi,bj)=0. _d 0
         OBWa0(j,bi,bj)=0. _d 0
         OBWh0(j,bi,bj)=0. _d 0
         OBWa1(j,bi,bj)=0. _d 0
         OBWh1(j,bi,bj)=0. _d 0
         OBWsl (j,bi,bj)=0. _d 0
         OBWsn (j,bi,bj)=0. _d 0
         OBWsl0(j,bi,bj)=0. _d 0
         OBWsn0(j,bi,bj)=0. _d 0
         OBWsl1(j,bi,bj)=0. _d 0
         OBWsn1(j,bi,bj)=0. _d 0
         OBWuice (j,bi,bj)=0. _d 0
         OBWvice (j,bi,bj)=0. _d 0
         OBWuice0(j,bi,bj)=0. _d 0
         OBWvice0(j,bi,bj)=0. _d 0
         OBWuice1(j,bi,bj)=0. _d 0
         OBWvice1(j,bi,bj)=0. _d 0
#endif /* ALLOW_OBCS_WEST */
        ENDDO
#endif /* ALLOW_SEAICE */

#ifdef ALLOW_PTRACERS
#ifndef ALLOW_AUTODIFF_TAMC
        IF ( usePTRACERS ) THEN
#endif
         DO iTracer=1,PTRACERS_numInUse
          DO k=1,Nr
           DO i=1-Olx,sNx+Olx
#ifdef ALLOW_OBCS_NORTH
            OBNptr (i,k,bi,bj,iTracer)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
            OBNptr0(i,k,bi,bj,iTracer)=0. _d 0
            OBNptr1(i,k,bi,bj,iTracer)=0. _d 0
# endif
#endif /* ALLOW_OBCS_NORTH */

#ifdef ALLOW_OBCS_SOUTH
            OBSptr (i,k,bi,bj,iTracer)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
            OBSptr0(i,k,bi,bj,iTracer)=0. _d 0
            OBSptr1(i,k,bi,bj,iTracer)=0. _d 0
# endif
#endif /* ALLOW_OBCS_SOUTH */
           ENDDO

           DO j=1-Oly,sNy+Oly
#ifdef ALLOW_OBCS_EAST
            OBEptr (j,k,bi,bj,iTracer)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
            OBEptr0(j,k,bi,bj,iTracer)=0. _d 0
            OBEptr1(j,k,bi,bj,iTracer)=0. _d 0
# endif
#endif /* ALLOW_OBCS_EAST */

#ifdef ALLOW_OBCS_WEST
            OBWptr (j,k,bi,bj,iTracer)=0. _d 0
# ifdef ALLOW_OBCS_PRESCRIBE
            OBWptr0(j,k,bi,bj,iTracer)=0. _d 0
            OBWptr1(j,k,bi,bj,iTracer)=0. _d 0
# endif
#endif /* ALLOW_OBCS_WEST */
           ENDDO
          ENDDO
         ENDDO
#ifndef ALLOW_AUTODIFF_TAMC
        ENDIF
#endif
#endif /* ALLOW_PTRACERS */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_ORLANSKI
        IF (useOrlanskiNorth.OR.useOrlanskiSouth.OR.
     &      useOrlanskiEast.OR.useOrlanskiWest) THEN
#ifdef ALLOW_DEBUG
          IF (debugMode) CALL DEBUG_CALL('ORLANSKI_INIT',myThid)
#endif
          CALL ORLANSKI_INIT(bi, bj, myThid)
        ENDIF
#endif /* ALLOW_ORLANSKI */

       ENDDO
      ENDDO

#ifdef ALLOW_OBCS_TIDES
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,tidalComponents
         DO i=1-Olx,sNx+Olx
# ifdef ALLOW_OBCS_NORTH
          OBNam (i,k,bi,bj)=0. _d 0
          OBNph (i,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_SOUTH
          OBSam (i,k,bi,bj)=0. _d 0
          OBSph (i,k,bi,bj)=0. _d 0
# endif
         ENDDO
         DO j=1-Oly,sNy+Oly
# ifdef ALLOW_OBCS_EAST
          OBEam (j,k,bi,bj)=0. _d 0
          OBEph (j,k,bi,bj)=0. _d 0
# endif
# ifdef ALLOW_OBCS_WEST
          OBWam (j,k,bi,bj)=0. _d 0
          OBWph (j,k,bi,bj)=0. _d 0
# endif
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _BARRIER
# ifdef ALLOW_OBCS_NORTH
      IF ( OBNamFile .NE. ' '  ) CALL READ_REC_XZ_RL
     &     (OBNamFile,fp,tidalComponents,OBNam,1,nIter0,myThid )
      IF ( OBNphFile .NE. ' '  ) CALL READ_REC_XZ_RL
     &     (OBNphFile,fp,tidalComponents,OBNph,1,nIter0,myThid )
# endif
# ifdef ALLOW_OBCS_SOUTH
      IF ( OBSamFile .NE. ' '  ) CALL READ_REC_XZ_RL
     &     (OBSamFile,fp,tidalComponents,OBSam,1,nIter0,myThid )
      IF ( OBSphFile .NE. ' '  ) CALL READ_REC_XZ_RL
     &     (OBSphFile,fp,tidalComponents,OBSph,1,nIter0,myThid )
# endif
# ifdef ALLOW_OBCS_EAST
      IF ( OBEamFile .NE. ' '  ) CALL READ_REC_YZ_RL
     &     (OBEamFile,fp,tidalComponents,OBEam,1,nIter0,myThid )
      IF ( OBEphFile .NE. ' '  ) CALL READ_REC_YZ_RL
     &     (OBEphFile,fp,tidalComponents,OBEph,1,nIter0,myThid )
# endif
# ifdef ALLOW_OBCS_WEST
      IF ( OBWamFile .NE. ' '  ) CALL READ_REC_YZ_RL
     &     (OBWamFile,fp,tidalComponents,OBWam,1,nIter0,myThid )
      IF ( OBWphFile .NE. ' '  ) CALL READ_REC_YZ_RL
     &     (OBWphFile,fp,tidalComponents,OBWph,1,nIter0,myThid )
# endif
      _BARRIER
#endif /* ALLOW_OBCS_TIDES */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C     Only needed for Orlanski:
      IF ( nIter0.NE.0 .OR. pickupSuff.NE.' ' ) THEN
        CALL OBCS_READ_PICKUP( nIter0, myThid )
      ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Load/compute OBCS values (initial conditions only)
      IF ( nIter0.EQ.0 .AND. pickupSuff.EQ.' ' ) THEN
#ifdef ALLOW_DEBUG
       IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
#endif
       CALL OBCS_CALC( startTime, nIter0,
     &              uVel, vVel, wVel, theta, salt, myThid )

C--   Apply OBCS values to initial conditions for consistency
C      (but initial conditions only)
#ifdef ALLOW_DEBUG
       IF (debugMode)
     &    CALL DEBUG_CALL('OBCS_APPLY_UV + OBCS_APPLY_TS',myThid)
#endif
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
          CALL OBCS_APPLY_UV( bi, bj, 0, uVel, vVel, myThid )
          CALL OBCS_APPLY_TS( bi, bj, 0, theta, salt, myThid )
        ENDDO
       ENDDO
       IF (useOBCSprescribe) THEN
C     After applying the boundary conditions exchange the 3D-fields.
C     This is only necessary of the boundary values have been read
C     from a file.
#ifdef ALLOW_DEBUG
        IF (debugMode)
     &    CALL DEBUG_CALL('EXCHANGES in OBCS_INIT_VARIABLES',myThid)
#endif
        CALL EXCH_UV_XYZ_RL(uVel,vVel,.TRUE.,myThid)
        _EXCH_XYZ_RL( theta, myThid )
        _EXCH_XYZ_RL( salt , myThid )
       ENDIF
C     endif start from rest
#ifdef ALLOW_OBCS_PRESCRIBE
      ELSEIF ( useOBCSprescribe ) THEN
C-- Need this call (although loading U,V,T,S OB values is not really necessary)
C    a) OB values needed for etaH (NonLinFreeSurf) and for wVel (nonHydrostatic)
C    b) OB values needed for ptracers (in case nIter0 = PTRACERS_Iter0 <> 0)
C    c) with present implementation, only do initialisation when called
C       with myTime=startTime (S/R EXF_GETFFIELDREC, setting "first")
C       or with myIter=nIter0 (S/R OBCS_FIELDS_LOAD). And this cannot be changed
C       because of above call to OBCS_CALC(startTime,nIter0).
        CALL OBCS_PRESCRIBE_READ( startTime, nIter0, myThid )
#endif
      ENDIF
C--   calls to S/R OBCS_COPY_TRACER on theta & salt are no longer needed
C     with maskInC,W,S in pkg/generic_advdiff: removed

#ifdef ALLOW_PTRACERS
C     repeat everything for passive tracers
      IF ( usePTRACERS ) THEN
C     catch the case when we do start from a pickup for dynamics variables
C     but initialise ptracers differently
       IF ( nIter0.EQ.PTRACERS_Iter0 ) THEN
#ifdef ALLOW_DEBUG
        IF (debugMode) CALL DEBUG_CALL('OBCS_APPLY_PTRACER',myThid)
#endif
        DO iTracer=1,PTRACERS_numInUse
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
            CALL OBCS_APPLY_PTRACER(
     I           bi, bj, 0, iTracer,
     U           ptracer(1-Olx,1-Oly,1,bi,bj,iTracer),
     I           myThid )
          ENDDO
         ENDDO
        ENDDO
C     endif start from rest
       ENDIF
       IF ( nIter0.EQ.PTRACERS_Iter0 .AND. useOBCSprescribe ) THEN
C     After applying the boundary conditions exchange the 3D-fields.
C     This is only necessary of the boundary values have been read
C     from a file.
#ifdef ALLOW_DEBUG
         IF (debugMode) CALL DEBUG_CALL(
     &        'PTRACERS EXCHANGES in OBCS_INIT_VARIABLES',myThid)
#endif
         CALL PTRACERS_FIELDS_BLOCKING_EXCH( myThid )
c      ELSE
C--   Calls to S/R OBCS_COPY_TRACER on pTracers are no longer needed
C     with maskInC,W,S in pkg/generic_advdiff: removed
C-    This call was part of ptracers exchange S/R but was needed in all cases
       ENDIF
C     endif usePTRACERS
      ENDIF
#endif /* ALLOW_PTRACERS */

#endif /* ALLOW_OBCS */

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_LEAVE('OBCS_INIT_VARIABLES',myThid)
#endif
      RETURN
      END