C $Header: /u/gcmpack/MITgcm/model/src/ini_depths.F,v 1.48 2010/01/21 00:15:31 jmc Exp $
C $Name: $
#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 ( debugLevel.GE.debLevB ) 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 )
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 ( debugLevel.GE.debLevB ) 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