aim.5l_Equatorial_Channel/code/ini_depths.F

717a660aff Jean*0001 #include "PACKAGES_CONFIG.h"
42c525bfb4 Alis*0002 #include "CPP_OPTIONS.h"
                0003
26b1c6e660 Jean*0004 CBOP
                0005 C     !ROUTINE: INI_DEPTHS
                0006 C     !INTERFACE:
42c525bfb4 Alis*0007       SUBROUTINE INI_DEPTHS( myThid )
486e38797b Jean*0008
26b1c6e660 Jean*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
8d70182ac2 Jean*0011 C     | SUBROUTINE INI_DEPTHS
                0012 C     | o define R_position of Lower and Surface Boundaries
26b1c6e660 Jean*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.
26b1c6e660 Jean*0027 C     *==========================================================*
                0028 C     \ev
                0029
                0030 C     !USES:
42c525bfb4 Alis*0031       IMPLICIT NONE
                0032 C     === Global variables ===
                0033 #include "SIZE.h"
                0034 #include "EEPARAMS.h"
                0035 #include "PARAMS.h"
                0036 #include "GRID.h"
62a035ab16 Jean*0037 #include "SURFACE.h"
12ffad7671 Jean*0038 #ifdef ALLOW_MNC
                0039 # include "MNC_PARAMS.h"
                0040 #endif
42c525bfb4 Alis*0041
26b1c6e660 Jean*0042 C     !INPUT/OUTPUT PARAMETERS:
42c525bfb4 Alis*0043 C     == Routine arguments ==
486e38797b Jean*0044 C     myThid  ::  my Thread Id number
42c525bfb4 Alis*0045       INTEGER myThid
                0046
26b1c6e660 Jean*0047 C     !LOCAL VARIABLES:
42c525bfb4 Alis*0048 C     == Local variables ==
486e38797b 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
42c525bfb4 Alis*0054       INTEGER iG, jG
                0055       INTEGER bi, bj
486e38797b Jean*0056       INTEGER  i, j
61c0ef52f3 Jean*0057       CHARACTER*(MAX_LEN_MBUF) msgBuf
26b1c6e660 Jean*0058 CEOP
61c0ef52f3 Jean*0059
8d70182ac2 Jean*0060       IF (usingPCoords .AND. bathyFile .NE. ' '
                0061      &                 .AND. topoFile  .NE. ' ' ) THEN
61c0ef52f3 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)
                0074         DO j=1-Oly,sNy+Oly
                0075          DO i=1-Olx,sNx+Olx
717a660aff 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
61c0ef52f3 Jean*0079          ENDDO
                0080         ENDDO
                0081        ENDDO
                0082       ENDDO
                0083
717a660aff Jean*0084 C-    Need to synchronize here before doing master-thread IO
486e38797b Jean*0085 C     this is done within IO routines => no longer needed
                0086 c     _BARRIER
717a660aff Jean*0087
61c0ef52f3 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
61c0ef52f3 Jean*0092 C- e.g., atmosphere : R_low = Top of atmosphere
                0093 C-            ocean : R_low = Bottom
42c525bfb4 Alis*0094        DO bj = myByLo(myThid), myByHi(myThid)
                0095         DO bi = myBxLo(myThid), myBxHi(myThid)
                0096          DO j=1,sNy
                0097           DO i=1,sNx
61c0ef52f3 Jean*0098            R_low(i,j,bi,bj) = rF(Nr+1)
42c525bfb4 Alis*0099           ENDDO
                0100          ENDDO
                0101         ENDDO
                0102        ENDDO
                0103       ELSE
12ffad7671 Jean*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)
486e38797b Jean*0112           CALL MNC_CW_RS_R('D',bathyFile,0,0,'bathy',R_low, myThid)
12ffad7671 Jean*0113           CALL MNC_FILE_CLOSE_ALL_MATCHING(bathyFile, myThid)
                0114           CALL MNC_CW_DEL_VNAME('bathy', myThid)
                0115         ELSE
                0116 #endif /*  ALLOW_MNC  */
717a660aff Jean*0117 C Read the bathymetry using the mid-level I/O package read_write_rec
42c525bfb4 Alis*0118 C The 0 is the "iteration" argument. The 1 is the record number.
12ffad7671 Jean*0119           CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
717a660aff Jean*0120 C Read the bathymetry using the mid-level I/O package read_write_fld
42c525bfb4 Alis*0121 C The 0 is the "iteration" argument. The ' ' is an empty suffix
61c0ef52f3 Jean*0122 c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
42c525bfb4 Alis*0123 C Read the bathymetry using the low-level I/O package
61c0ef52f3 Jean*0124 c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
                0125 c    &                     'RS', 1, R_low, 1, myThid )
12ffad7671 Jean*0126
                0127 #ifdef ALLOW_MNC
                0128         ENDIF
                0129 #endif /*  ALLOW_MNC  */
                0130
                0131
42c525bfb4 Alis*0132       ENDIF
61c0ef52f3 Jean*0133 C- end setup R_low in the interior
                0134
717a660aff Jean*0135 C- fill in the overlap (+ BARRIER):
12ffad7671 Jean*0136       _EXCH_XY_RS(R_low, myThid )
61c0ef52f3 Jean*0137
b090e8cd4a Jean*0138       IF ( debugLevel.GE.debLevC ) THEN
717a660aff Jean*0139 c      PRINT *, ' Calling plot field', myThid
                0140        CALL PLOT_FIELD_XYRS( R_low, 'Bottom depths (ini_depths)',
                0141      &                       -1, myThid )
                0142       ENDIF
8d70182ac2 Jean*0143 c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
61c0ef52f3 Jean*0144
                0145 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0146
                0147 C------
                0148 C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
                0149 C------
                0150
8d70182ac2 Jean*0151       IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
                0152 C------ read directly Po_surf from bathyFile (only for backward compatibility)
42c525bfb4 Alis*0153
61c0ef52f3 Jean*0154         CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
                0155
                0156       ELSEIF ( topoFile.EQ.' ' ) THEN
                0157 C------ set default value:
                0158
                0159         DO bj = myByLo(myThid), myByHi(myThid)
                0160          DO bi = myBxLo(myThid), myBxHi(myThid)
                0161           DO j=1,sNy
                0162            DO i=1,sNx
12ffad7671 Jean*0163             Ro_surf(i,j,bi,bj) = rF(1)
61c0ef52f3 Jean*0164            ENDDO
                0165           ENDDO
                0166          ENDDO
                0167         ENDDO
42c525bfb4 Alis*0168
61c0ef52f3 Jean*0169       ELSE
                0170 C------ read from file:
                0171
                0172 C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
717a660aff Jean*0173
62a035ab16 Jean*0174         CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
61c0ef52f3 Jean*0175
                0176         IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
                0177 C----
8d70182ac2 Jean*0178 C   Convert Surface Geopotential to (reference) Surface Pressure
61c0ef52f3 Jean*0179 C   according to Tref profile, using same discretisation as in calc_phi_hyd
                0180 C----
8d70182ac2 Jean*0181 c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
61c0ef52f3 Jean*0182
6177f7885b Jean*0183           CALL INI_P_GROUND( 2, topoZ,
c7734eb7a9 Jean*0184      O                       Ro_surf,
                0185      I                       myThid )
61c0ef52f3 Jean*0186
8d70182ac2 Jean*0187 C         This I/O is now done in write_grid.F
                0188 c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
61c0ef52f3 Jean*0189
12ffad7671 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
61c0ef52f3 Jean*0197         ELSE
                0198 C----
8d70182ac2 Jean*0199 C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
6177f7885b Jean*0200 C    below an ice-shelf - NOTE - actually not yet implemented )
61c0ef52f3 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
62a035ab16 Jean*0205               Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
61c0ef52f3 Jean*0206              ENDDO
                0207             ENDDO
                0208            ENDDO
                0209           ENDDO
                0210
                0211         ENDIF
                0212
                0213 C------ end case "read topoFile"
                0214       ENDIF
                0215
717a660aff Jean*0216 C----- fill in the overlap (+ BARRIER):
12ffad7671 Jean*0217       _EXCH_XY_RS(Ro_surf, myThid )
61c0ef52f3 Jean*0218
                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
42c525bfb4 Alis*0225        DO bj = myByLo(myThid), myByHi(myThid)
                0226         DO bi = myBxLo(myThid), myBxHi(myThid)
                0227          DO j=1-Oly,sNy+Oly
                0228           DO i=1-Olx,sNx+Olx
486e38797b Jean*0229            iG = myXGlobalLo-1+(bi-1)*sNx+i
                0230            jG = myYGlobalLo-1+(bj-1)*sNy+j
61c0ef52f3 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.
486e38797b Jean*0235 C- Domain : Symetric % Eq. & closed at N & S boundaries:
                0236            IF ( usingSphericalPolarGrid .AND.
                0237      &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
                0238      &       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)
                0249          DO j=1-Oly,sNy+Oly
                0250           DO i=1-Olx,sNx+Olx
486e38797b Jean*0251            iG = myXGlobalLo-1+(bi-1)*sNx+i
                0252            jG = myYGlobalLo-1+(bj-1)*sNy+j
61c0ef52f3 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.
486e38797b Jean*0257 C- Domain : Symetric % Eq. & closed at N & S boundaries:
                0258            IF ( usingSphericalPolarGrid .AND.
                0259      &          ABS(yC(i,j,bi,bj)).GE.ABS(ygOrigin) )
                0260      &       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)
42c525bfb4 Alis*0264           ENDDO
                0265          ENDDO
                0266         ENDDO
                0267        ENDDO
                0268       ENDIF
aea29c8517 Alis*0269
b090e8cd4a Jean*0270       IF ( debugLevel.GE.debLevC ) THEN
486e38797b Jean*0271        _BARRIER
717a660aff 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)
717a660aff Jean*0277
                0278 C--   Everyone else must wait for the depth to be loaded
486e38797b Jean*0279 C-    note: not necessary since all single-thread IO above are followed
                0280 C           by an EXCH (with BARRIER) + BARRIER within IO routine
                0281 c     _BARRIER
61c0ef52f3 Jean*0282
12ffad7671 Jean*0283 #ifdef ALLOW_OBCS
                0284       IF ( useOBCS ) THEN
                0285 C     check for inconsistent topography along boundaries and fix it
                0286        CALL OBCS_CHECK_DEPTHS( myThid )
                0287 C     update the overlaps
486e38797b Jean*0288        _EXCH_XY_RS( R_low, myThid )
12ffad7671 Jean*0289       ENDIF
                0290 #endif /* ALLOW_OBCS */
                0291
                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
42c525bfb4 Alis*0297       RETURN
                0298       END