model/src/ini_depths.F

a6cbc7a360 Mart*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
924557e60a Chri*0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_DEPTHS
                0006 C     !INTERFACE:
924557e60a Chri*0007       SUBROUTINE INI_DEPTHS( myThid )
fb2f11e499 Jean*0008
9366854e02 Chri*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
8d70182ac2 Jean*0011 C     | SUBROUTINE INI_DEPTHS
                0012 C     | o define R_position of Lower and Surface Boundaries
9366854e02 Chri*0013 C     *==========================================================*
8d70182ac2 Jean*0014 C     |atmosphere orography:
                0015 C     | define either in term of P_topo or converted from Z_topo
                0016 C     |ocean bathymetry:
                0017 C     | The depths of the bottom of the model is specified in
                0018 C     | terms of an XY map with one depth for each column of
                0019 C     | grid cells. Depths do not have to coincide with the
                0020 C     | model levels. The model lopping algorithm makes it
                0021 C     | possible to represent arbitrary depths.
                0022 C     | The mode depths map also influences the models topology
                0023 C     | By default the model domain wraps around in X and Y.
                0024 C     | This default doubly periodic topology is "supressed"
                0025 C     | if a depth map is defined which closes off all wrap
                0026 C     | around flow.
9366854e02 Chri*0027 C     *==========================================================*
                0028 C     \ev
924557e60a Chri*0029
9366854e02 Chri*0030 C     !USES:
                0031       IMPLICIT NONE
924557e60a Chri*0032 C     === Global variables ===
                0033 #include "SIZE.h"
                0034 #include "EEPARAMS.h"
                0035 #include "PARAMS.h"
                0036 #include "GRID.h"
681eca3696 Jean*0037 #include "SURFACE.h"
5f4df5533c Ed H*0038 #ifdef ALLOW_MNC
76da4a5e0b Mart*0039 # include "MNC_PARAMS.h"
5f4df5533c Ed H*0040 #endif
924557e60a Chri*0041
9366854e02 Chri*0042 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0043 C     == Routine arguments ==
fb2f11e499 Jean*0044 C     myThid  ::  my Thread Id number
924557e60a Chri*0045       INTEGER myThid
                0046
9366854e02 Chri*0047 C     !LOCAL VARIABLES:
924557e60a Chri*0048 C     == Local variables ==
fb2f11e499 Jean*0049 C     iG, jG  :: Global coordinate index
                0050 C     bi, bj  :: Tile indices
                0051 C     i, j    :: Loop counters
                0052 C     oldPrec :: Temporary used in controlling binary input dataset precision
                0053 C     msgBuf  :: Informational/error message buffer
924557e60a Chri*0054       INTEGER iG, jG
                0055       INTEGER bi, bj
fb2f11e499 Jean*0056       INTEGER  i, j
6c9bcec054 Jean*0057       CHARACTER*(MAX_LEN_MBUF) msgBuf
9366854e02 Chri*0058 CEOP
6c9bcec054 Jean*0059
8d70182ac2 Jean*0060       IF (usingPCoords .AND. bathyFile .NE. ' '
                0061      &                 .AND. topoFile  .NE. ' ' ) THEN
6c9bcec054 Jean*0062        WRITE(msgBuf,'(A,A)')
                0063      &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
                0064      &  ' select the right one !'
                0065        CALL PRINT_ERROR( msgBuf , myThid)
                0066        STOP 'ABNORMAL END: S/R INI_DEPTHS'
                0067       ENDIF
                0068
                0069 C------
                0070 C   0) Initialize R_low and Ro_surf (define an empty domain)
                0071 C------
                0072       DO bj = myByLo(myThid), myByHi(myThid)
                0073        DO bi = myBxLo(myThid), myBxHi(myThid)
522c728681 Jean*0074         DO j=1-OLy,sNy+OLy
                0075          DO i=1-OLx,sNx+OLx
c2953005c8 Jean*0076           R_low(i,j,bi,bj)   = 0. _d 0
                0077           Ro_surf(i,j,bi,bj) = 0. _d 0
                0078           topoZ(i,j,bi,bj)   = 0. _d 0
6c9bcec054 Jean*0079          ENDDO
                0080         ENDDO
                0081        ENDDO
                0082       ENDDO
                0083
c2953005c8 Jean*0084 C-    Need to synchronize here before doing master-thread IO
3365bdc872 Jean*0085 C     this is done within IO routines => no longer needed
                0086 c     _BARRIER
c2953005c8 Jean*0087
6c9bcec054 Jean*0088 C------
                0089 C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
                0090 C------
8d70182ac2 Jean*0091       IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
6c9bcec054 Jean*0092 C- e.g., atmosphere : R_low = Top of atmosphere
6bf17245c3 Jean*0093 C-            ocean : R_low = Bottom
81bc00c2f0 Chri*0094        DO bj = myByLo(myThid), myByHi(myThid)
                0095         DO bi = myBxLo(myThid), myBxHi(myThid)
6bf17245c3 Jean*0096          DO j=1,sNy
                0097           DO i=1,sNx
6c9bcec054 Jean*0098            R_low(i,j,bi,bj) = rF(Nr+1)
81bc00c2f0 Chri*0099           ENDDO
924557e60a Chri*0100          ENDDO
                0101         ENDDO
                0102        ENDDO
81bc00c2f0 Chri*0103       ELSE
5f4df5533c Ed H*0104
                0105 #ifdef ALLOW_MNC
                0106         IF (useMNC .AND. mnc_read_bathy) THEN
                0107           CALL MNC_CW_ADD_VNAME('bathy', 'Cen_xy_Hn__-__-', 3,4, myThid)
                0108           CALL MNC_FILE_CLOSE_ALL_MATCHING(bathyFile, myThid)
                0109           CALL MNC_CW_SET_UDIM(bathyFile, 1, myThid)
                0110           CALL MNC_CW_SET_CITER(bathyFile, 2, -1, -1, -1, myThid)
                0111           CALL MNC_CW_SET_UDIM(bathyFile, 1, myThid)
fa6c69f783 Jean*0112           CALL MNC_CW_RS_R('D',bathyFile,0,0,'bathy',R_low, myThid)
5f4df5533c Ed H*0113           CALL MNC_FILE_CLOSE_ALL_MATCHING(bathyFile, myThid)
                0114           CALL MNC_CW_DEL_VNAME('bathy', myThid)
                0115         ELSE
                0116 #endif /*  ALLOW_MNC  */
c2953005c8 Jean*0117 C Read the bathymetry using the mid-level I/O package read_write_rec
ab42872a05 Alis*0118 C The 0 is the "iteration" argument. The 1 is the record number.
5f4df5533c Ed H*0119           CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
c2953005c8 Jean*0120 C Read the bathymetry using the mid-level I/O package read_write_fld
ab42872a05 Alis*0121 C The 0 is the "iteration" argument. The ' ' is an empty suffix
6c9bcec054 Jean*0122 c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
ab42872a05 Alis*0123 C Read the bathymetry using the low-level I/O package
6c9bcec054 Jean*0124 c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
                0125 c    &                     'RS', 1, R_low, 1, myThid )
5f4df5533c Ed H*0126
                0127 #ifdef ALLOW_MNC
                0128         ENDIF
                0129 #endif /*  ALLOW_MNC  */
                0130
81bc00c2f0 Chri*0131       ENDIF
6bf17245c3 Jean*0132 C- end setup R_low in the interior
                0133
c2953005c8 Jean*0134 C- fill in the overlap (+ BARRIER):
12c8b75709 Jean*0135       _EXCH_XY_RS(R_low, myThid )
fb481a83c2 Alis*0136
522c728681 Jean*0137       IF ( plotLevel.GE.debLevC ) THEN
f6070d4bb4 Patr*0138 c      PRINT *, ' Calling plot field', myThid
9f60bc5a7c Jean*0139        CALL PLOT_FIELD_XYRS( R_low, 'Bottom depths (ini_depths)',
                0140      &                       -1, myThid )
                0141       ENDIF
8d70182ac2 Jean*0142 c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
6c9bcec054 Jean*0143
                0144 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
fb481a83c2 Alis*0145
6c9bcec054 Jean*0146 C------
                0147 C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
                0148 C------
                0149
8d70182ac2 Jean*0150       IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
                0151 C------ read directly Po_surf from bathyFile (only for backward compatibility)
6bf17245c3 Jean*0152
                0153         CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
                0154
                0155       ELSEIF ( topoFile.EQ.' ' ) THEN
                0156 C------ set default value:
6c9bcec054 Jean*0157
                0158         DO bj = myByLo(myThid), myByHi(myThid)
                0159          DO bi = myBxLo(myThid), myBxHi(myThid)
6bf17245c3 Jean*0160           DO j=1,sNy
                0161            DO i=1,sNx
ccf58cd9c3 Jean*0162             Ro_surf(i,j,bi,bj) = rF(1)
6c9bcec054 Jean*0163            ENDDO
                0164           ENDDO
                0165          ENDDO
                0166         ENDDO
                0167
                0168       ELSE
6bf17245c3 Jean*0169 C------ read from file:
6c9bcec054 Jean*0170
                0171 C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
c2953005c8 Jean*0172
681eca3696 Jean*0173         CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
32faf9b967 Jean*0174         _EXCH_XY_RS( topoZ, myThid )
6c9bcec054 Jean*0175
                0176         IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
6bf17245c3 Jean*0177 C----
8d70182ac2 Jean*0178 C   Convert Surface Geopotential to (reference) Surface Pressure
6c9bcec054 Jean*0179 C   according to Tref profile, using same discretisation as in calc_phi_hyd
6bf17245c3 Jean*0180 C----
8d70182ac2 Jean*0181 c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
6c9bcec054 Jean*0182
46ac67bdf3 Jean*0183           CALL INI_P_GROUND( 2, topoZ,
                0184      O                       Ro_surf,
463053c692 Jean*0185      I                       myThid )
6c9bcec054 Jean*0186
f4a7634227 Alis*0187 C         This I/O is now done in write_grid.F
8d70182ac2 Jean*0188 c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
6c9bcec054 Jean*0189
ccf58cd9c3 Jean*0190         ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
                0191
                0192           WRITE(msgBuf,'(A,A)') 'S/R INI_DEPTHS: ',
                0193      &     'from topoFile (in m) to ref.bottom pressure: Not yet coded'
                0194           CALL PRINT_ERROR( msgBuf , myThid)
                0195           STOP 'ABNORMAL END: S/R INI_DEPTHS'
                0196
6c9bcec054 Jean*0197         ELSE
6bf17245c3 Jean*0198 C----
8d70182ac2 Jean*0199 C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
46ac67bdf3 Jean*0200 C    below an ice-shelf - NOTE - actually not yet implemented )
6c9bcec054 Jean*0201           DO bj = myByLo(myThid), myByHi(myThid)
                0202            DO bi = myBxLo(myThid), myBxHi(myThid)
                0203             DO j=1,sNy
                0204              DO i=1,sNx
681eca3696 Jean*0205               Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
6c9bcec054 Jean*0206              ENDDO
                0207             ENDDO
                0208            ENDDO
                0209           ENDDO
                0210
6bf17245c3 Jean*0211         ENDIF
6c9bcec054 Jean*0212
6bf17245c3 Jean*0213 C------ end case "read topoFile"
6c9bcec054 Jean*0214       ENDIF
                0215
c2953005c8 Jean*0216 C----- fill in the overlap (+ BARRIER):
12c8b75709 Jean*0217       _EXCH_XY_RS(Ro_surf, myThid )
6bf17245c3 Jean*0218
6c9bcec054 Jean*0219 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0220
                0221 C------
                0222 C   3) Close the Domain (special configuration).
                0223 C------
8d70182ac2 Jean*0224       IF (usingPCoords) THEN
fb481a83c2 Alis*0225        DO bj = myByLo(myThid), myByHi(myThid)
                0226         DO bi = myBxLo(myThid), myBxHi(myThid)
522c728681 Jean*0227          DO j=1-OLy,sNy+OLy
                0228           DO i=1-OLx,sNx+OLx
fb2f11e499 Jean*0229            iG = myXGlobalLo-1+(bi-1)*sNx+i
                0230            jG = myYGlobalLo-1+(bj-1)*sNy+j
6c9bcec054 Jean*0231 C          Test for eastern edge
                0232 c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
                0233 C          Test for northern edge
                0234 c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
fb2f11e499 Jean*0235 C- Domain : Symetric % Eq. & closed at N & S boundaries:
                0236 c          IF ( usingSphericalPolarGrid .AND.
                0237 c    &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
                0238 c    &       Ro_surf(i,j,bi,bj) = rF(Nr+1)
                0239            IF ( usingSphericalPolarGrid .AND.
                0240      &          ABS(yC(i,j,bi,bj)).GE.90. )
                0241      &       Ro_surf(i,j,bi,bj) = rF(Nr+1)
aea29c8517 Alis*0242           ENDDO
                0243          ENDDO
                0244         ENDDO
                0245        ENDDO
                0246       ELSE
                0247        DO bj = myByLo(myThid), myByHi(myThid)
                0248         DO bi = myBxLo(myThid), myBxHi(myThid)
522c728681 Jean*0249          DO j=1-OLy,sNy+OLy
                0250           DO i=1-OLx,sNx+OLx
fb2f11e499 Jean*0251            iG = myXGlobalLo-1+(bi-1)*sNx+i
                0252            jG = myYGlobalLo-1+(bj-1)*sNy+j
6c9bcec054 Jean*0253 C          Test for eastern edge
                0254 c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
                0255 C          Test for northern edge
                0256 c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
fb2f11e499 Jean*0257 C- Domain : Symetric % Eq. & closed at N & S boundaries:
                0258 c          IF ( usingSphericalPolarGrid .AND.
                0259 c    &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
                0260 c    &       R_low(i,j,bi,bj) = rF(1)
                0261            IF ( usingSphericalPolarGrid .AND.
                0262      &          ABS(yC(i,j,bi,bj)).GE.90. )
                0263      &       R_low(i,j,bi,bj) = rF(1)
fb481a83c2 Alis*0264           ENDDO
                0265          ENDDO
                0266         ENDDO
                0267        ENDDO
                0268       ENDIF
aea29c8517 Alis*0269
522c728681 Jean*0270       IF ( plotLevel.GE.debLevC ) THEN
3365bdc872 Jean*0271        _BARRIER
9f60bc5a7c Jean*0272        CALL PLOT_FIELD_XYRS( Ro_surf,
                0273      &           'Surface reference r-position (ini_depths)',
                0274      &           -1, myThid )
                0275       ENDIF
8d70182ac2 Jean*0276 c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
6c9bcec054 Jean*0277
06bb0cec77 Jean*0278 C--   Everyone else must wait for the depth to be loaded
fb2f11e499 Jean*0279 C-    note: not necessary since all single-thread IO above are followed
3365bdc872 Jean*0280 C           by an EXCH (with BARRIER) + BARRIER within IO routine
                0281 c     _BARRIER
06bb0cec77 Jean*0282
76da4a5e0b Mart*0283 #ifdef ALLOW_OBCS
                0284       IF ( useOBCS ) THEN
                0285 C     check for inconsistent topography along boundaries and fix it
0916a7672c Jean*0286        CALL OBCS_CHECK_DEPTHS( myThid )
76da4a5e0b Mart*0287 C     update the overlaps
fb2f11e499 Jean*0288        _EXCH_XY_RS( R_low, myThid )
76da4a5e0b Mart*0289       ENDIF
                0290 #endif /* ALLOW_OBCS */
                0291
0916a7672c Jean*0292 #ifdef ALLOW_EXCH2
                0293 C     Check domain boundary (e.g., in case of missing tiles)
                0294       CALL EXCH2_CHECK_DEPTHS( R_low, Ro_surf, myThid )
                0295 #endif
                0296
924557e60a Chri*0297       RETURN
                0298       END