#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

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

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_DEPTHS
C     | o define R_position of Lower and Surface Boundaries
C     *==========================================================*
C     |atmosphere orography:
C     | define either in term of P_topo or converted from Z_topo
C     |ocean bathymetry:
C     | The depths of the bottom of the model is specified in
C     | terms of an XY map with one depth for each column of
C     | grid cells. Depths do not have to coincide with the
C     | model levels. The model lopping algorithm makes it
C     | possible to represent arbitrary depths.
C     | The mode depths map also influences the models topology
C     | By default the model domain wraps around in X and Y.
C     | This default doubly periodic topology is "supressed"
C     | if a depth map is defined which closes off all wrap
C     | around flow.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#ifdef ALLOW_MNC
# include "MNC_PARAMS.h"
#endif

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

C     !LOCAL VARIABLES:
C     == Local variables ==
C     iG, jG  :: Global coordinate index
C     bi, bj  :: Tile indices
C     i, j    :: Loop counters
C     oldPrec :: Temporary used in controlling binary input dataset precision
C     msgBuf  :: Informational/error message buffer
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER  i, j
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      IF (usingPCoords .AND. bathyFile .NE. ' '
     &                 .AND. topoFile  .NE. ' ' ) THEN
       WRITE(msgBuf,'(A,A)')
     &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
     &  ' select the right one !'
       CALL PRINT_ERROR( msgBuf , myThid)
       STOP 'ABNORMAL END: S/R INI_DEPTHS'
      ENDIF

C------
C   0) Initialize R_low and Ro_surf (define an empty domain)
C------
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          R_low(i,j,bi,bj)   = 0. _d 0
          Ro_surf(i,j,bi,bj) = 0. _d 0
          topoZ(i,j,bi,bj)   = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C-    Need to synchronize here before doing master-thread IO
C     this is done within IO routines => no longer needed
c     _BARRIER

C------
C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
C------
      IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
C- e.g., atmosphere : R_low = Top of atmosphere
C-            ocean : R_low = Bottom
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           R_low(i,j,bi,bj) = rF(Nr+1)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ELSE

#ifdef ALLOW_MNC
        IF (useMNC .AND. mnc_read_bathy) THEN
          CALL MNC_CW_ADD_VNAME('bathy', 'Cen_xy_Hn__-__-', 3,4, myThid)
          CALL MNC_FILE_CLOSE_ALL_MATCHING(bathyFile, myThid)
          CALL MNC_CW_SET_UDIM(bathyFile, 1, myThid)
          CALL MNC_CW_SET_CITER(bathyFile, 2, -1, -1, -1, myThid)
          CALL MNC_CW_SET_UDIM(bathyFile, 1, myThid)
          CALL MNC_CW_RS_R('D',bathyFile,0,0,'bathy',R_low, myThid)
          CALL MNC_FILE_CLOSE_ALL_MATCHING(bathyFile, myThid)
          CALL MNC_CW_DEL_VNAME('bathy', myThid)
        ELSE
#endif /*  ALLOW_MNC  */
C Read the bathymetry using the mid-level I/O package read_write_rec
C The 0 is the "iteration" argument. The 1 is the record number.
          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
C Read the bathymetry using the mid-level I/O package read_write_fld
C The 0 is the "iteration" argument. The ' ' is an empty suffix
c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
C Read the bathymetry using the low-level I/O package
c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
c    &                     'RS', 1, R_low, 1, myThid )

#ifdef ALLOW_MNC
        ENDIF
#endif /*  ALLOW_MNC  */

      ENDIF
C- end setup R_low in the interior

C- fill in the overlap (+ BARRIER):
      _EXCH_XY_RS(R_low, myThid )

      IF ( plotLevel.GE.debLevC ) THEN
c      PRINT *, ' Calling plot field', myThid
       CALL PLOT_FIELD_XYRS( R_low, 'Bottom depths (ini_depths)',
     &                       -1, myThid )
      ENDIF
c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)

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

C------
C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
C------

      IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
C------ read directly Po_surf from bathyFile (only for backward compatibility)

        CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )

      ELSEIF ( topoFile.EQ.' ' ) THEN
C------ set default value:

        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            Ro_surf(i,j,bi,bj) = rF(1)
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      ELSE
C------ read from file:

C- read surface topography (in m) from topoFile (case topoFile.NE.' '):

        CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
        _EXCH_XY_RS( topoZ, myThid )

        IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
C----
C   Convert Surface Geopotential to (reference) Surface Pressure
C   according to Tref profile, using same discretisation as in calc_phi_hyd
C----
c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)

          CALL INI_P_GROUND( 2, topoZ,
     O                       Ro_surf,
     I                       myThid )

C         This I/O is now done in write_grid.F
c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)

        ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN

          WRITE(msgBuf,'(A,A)') 'S/R INI_DEPTHS: ',
     &     'from topoFile (in m) to ref.bottom pressure: Not yet coded'
          CALL PRINT_ERROR( msgBuf , myThid)
          STOP 'ABNORMAL END: S/R INI_DEPTHS'

        ELSE
C----
C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
C    below an ice-shelf - NOTE - actually not yet implemented )
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1,sNy
             DO i=1,sNx
              Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
             ENDDO
            ENDDO
           ENDDO
          ENDDO

        ENDIF

C------ end case "read topoFile"
      ENDIF

C----- fill in the overlap (+ BARRIER):
      _EXCH_XY_RS(Ro_surf, myThid )

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

C------
C   3) Close the Domain (special configuration).
C------
      IF (usingPCoords) THEN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           iG = myXGlobalLo-1+(bi-1)*sNx+i
           jG = myYGlobalLo-1+(bj-1)*sNy+j
C          Test for eastern edge
c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
C          Test for northern edge
c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
C- Domain : Symetric % Eq. & closed at N & S boundaries:
c          IF ( usingSphericalPolarGrid .AND.
c    &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
c    &       Ro_surf(i,j,bi,bj) = rF(Nr+1)
           IF ( usingSphericalPolarGrid .AND.
     &          ABS(yC(i,j,bi,bj)).GE.90. )
     &       Ro_surf(i,j,bi,bj) = rF(Nr+1)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ELSE
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           iG = myXGlobalLo-1+(bi-1)*sNx+i
           jG = myYGlobalLo-1+(bj-1)*sNy+j
C          Test for eastern edge
c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
C          Test for northern edge
c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
C- Domain : Symetric % Eq. & closed at N & S boundaries:
c          IF ( usingSphericalPolarGrid .AND.
c    &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
c    &       R_low(i,j,bi,bj) = rF(1)
           IF ( usingSphericalPolarGrid .AND.
     &          ABS(yC(i,j,bi,bj)).GE.90. )
     &       R_low(i,j,bi,bj) = rF(1)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

      IF ( plotLevel.GE.debLevC ) THEN
       _BARRIER
       CALL PLOT_FIELD_XYRS( Ro_surf,
     &           'Surface reference r-position (ini_depths)',
     &           -1, myThid )
      ENDIF
c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)

C--   Everyone else must wait for the depth to be loaded
C-    note: not necessary since all single-thread IO above are followed
C           by an EXCH (with BARRIER) + BARRIER within IO routine
c     _BARRIER

#ifdef ALLOW_OBCS
      IF ( useOBCS ) THEN
C     check for inconsistent topography along boundaries and fix it
       CALL OBCS_CHECK_DEPTHS( myThid )
C     update the overlaps
       _EXCH_XY_RS( R_low, myThid )
      ENDIF
#endif /* ALLOW_OBCS */

#ifdef ALLOW_EXCH2
C     Check domain boundary (e.g., in case of missing tiles)
      CALL EXCH2_CHECK_DEPTHS( R_low, Ro_surf, myThid )
#endif

      RETURN
      END
