/[MITgcm]/MITgcm/verification/aim.5l_Equatorial_Channel/code/ini_depths.F
ViewVC logotype

Diff of /MITgcm/verification/aim.5l_Equatorial_Channel/code/ini_depths.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.9 by jmc, Mon Mar 17 16:42:07 2003 UTC revision 1.11 by jmc, Wed Oct 18 20:44:37 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 9  C     !INTERFACE: Line 10  C     !INTERFACE:
10        SUBROUTINE INI_DEPTHS( myThid )        SUBROUTINE INI_DEPTHS( myThid )
11  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | SUBROUTINE INI_DEPTHS                                      C     | SUBROUTINE INI_DEPTHS
14  C     | o define R_position of Lower and Surface Boundaries        C     | o define R_position of Lower and Surface Boundaries
15  C     *==========================================================*  C     *==========================================================*
16  C     |atmosphere orography:                                        C     |atmosphere orography:
17  C     | define either in term of P_topo or converted from Z_topo    C     | define either in term of P_topo or converted from Z_topo
18  C     |ocean bathymetry:                                            C     |ocean bathymetry:
19  C     | The depths of the bottom of the model is specified in      C     | The depths of the bottom of the model is specified in
20  C     | terms of an XY map with one depth for each column of        C     | terms of an XY map with one depth for each column of
21  C     | grid cells. Depths do not have to coincide with the        C     | grid cells. Depths do not have to coincide with the
22  C     | model levels. The model lopping algorithm makes it          C     | model levels. The model lopping algorithm makes it
23  C     | possible to represent arbitrary depths.                    C     | possible to represent arbitrary depths.
24  C     | The mode depths map also influences the models topology    C     | The mode depths map also influences the models topology
25  C     | By default the model domain wraps around in X and Y.        C     | By default the model domain wraps around in X and Y.
26  C     | This default doubly periodic topology is "supressed"        C     | This default doubly periodic topology is "supressed"
27  C     | if a depth map is defined which closes off all wrap        C     | if a depth map is defined which closes off all wrap
28  C     | around flow.                                                C     | around flow.
29  C     *==========================================================*  C     *==========================================================*
30  C     \ev  C     \ev
31    
# Line 36  C     === Global variables === Line 37  C     === Global variables ===
37  #include "PARAMS.h"  #include "PARAMS.h"
38  #include "GRID.h"  #include "GRID.h"
39  #include "SURFACE.h"  #include "SURFACE.h"
40    #ifdef ALLOW_SHELFICE
41    # include "SHELFICE.h"
42    #endif /* ALLOW_SHELFICE */
43    
44  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
45  C     == Routine arguments ==  C     == Routine arguments ==
# Line 49  C     iG, jG - Global coordinate index Line 53  C     iG, jG - Global coordinate index
53  C     bi,bj  - Tile indices  C     bi,bj  - Tile indices
54  C     I,J,K  - Loop counters  C     I,J,K  - Loop counters
55  C     oldPrec - Temporary used in controlling binary input dataset precision  C     oldPrec - Temporary used in controlling binary input dataset precision
56  C     msgBuf    - Informational/error meesage buffer  C     msgBuf    - Informational/error meesage buffer
57        INTEGER iG, jG        INTEGER iG, jG
58        INTEGER bi, bj        INTEGER bi, bj
59        INTEGER  I, J        INTEGER  I, J
60        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
61  CEOP  CEOP
62    
63        IF (groundAtK1 .AND. bathyFile .NE. ' '        IF (usingPCoords .AND. bathyFile .NE. ' '
64       &               .AND. topoFile  .NE. ' ' ) THEN       &                 .AND. topoFile  .NE. ' ' ) THEN
65         WRITE(msgBuf,'(A,A)')         WRITE(msgBuf,'(A,A)')
66       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
67       &  ' select the right one !'       &  ' select the right one !'
# Line 72  C------ Line 76  C------
76         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
77          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
78           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
79            R_low(i,j,bi,bj) = 0.            R_low(i,j,bi,bj)   = 0. _d 0
80            Ro_surf(i,j,bi,bj) = 0.            Ro_surf(i,j,bi,bj) = 0. _d 0
81            topoZ(i,j,bi,bj) = 0.            topoZ(i,j,bi,bj)   = 0. _d 0
82    #ifdef ALLOW_SHELFICE
83              R_shelfIce(i,j,bi,bj) = 0. _d 0
84    #endif /* ALLOW_SHELFICE */
85           ENDDO           ENDDO
86          ENDDO          ENDDO
87         ENDDO         ENDDO
88        ENDDO        ENDDO
89    
90    C-    Need to synchronize here before doing master-thread IO
91          _BARRIER
92    
93  C------  C------
94  C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :  C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
95  C------  C------
96        IF (groundAtK1 .OR. bathyFile .EQ. ' ') THEN        IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
97  C- e.g., atmosphere : R_low = Top of atmosphere  C- e.g., atmosphere : R_low = Top of atmosphere
98  C-            ocean : R_low = Bottom  C-            ocean : R_low = Bottom
99         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
# Line 96  C-            ocean : R_low = Bottom Line 106  C-            ocean : R_low = Bottom
106          ENDDO          ENDDO
107         ENDDO         ENDDO
108        ELSE        ELSE
109          _BEGIN_MASTER( myThid )  C Read the bathymetry using the mid-level I/O package read_write_rec
 C Read the bathymetry using the mid-level I/O pacakage read_write_rec  
110  C The 0 is the "iteration" argument. The 1 is the record number.  C The 0 is the "iteration" argument. The 1 is the record number.
111          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
112  C Read the bathymetry using the mid-level I/O pacakage read_write_fld  C Read the bathymetry using the mid-level I/O package read_write_fld
113  C The 0 is the "iteration" argument. The ' ' is an empty suffix  C The 0 is the "iteration" argument. The ' ' is an empty suffix
114  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
115  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
116  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
117  c    &                     'RS', 1, R_low, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
         _END_MASTER(myThid)  
   
118        ENDIF        ENDIF
119  C- end setup R_low in the interior  C- end setup R_low in the interior
120    
121  C- fill in the overlap :  C- fill in the overlap (+ BARRIER):
122        _EXCH_XY_R4(R_low, myThid )        _EXCH_XY_R4(R_low, myThid )
123    
124  c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)        IF ( debugLevel.GE.debLevB ) THEN
125  c     _BEGIN_MASTER( myThid )  c      PRINT *, ' Calling plot field', myThid
126  c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)         CALL PLOT_FIELD_XYRS( R_low, 'Bottom depths (ini_depths)',
127  c     _END_MASTER(myThid)       &                       -1, myThid )
128          ENDIF
129    c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
130    
131  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
132    
# Line 125  C------ Line 134  C------
134  C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere  C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
135  C------  C------
136    
137        IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN        IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
138  C------ read directly Po_surf from bathyFile (only for backward compatibility)  C------ read directly Po_surf from bathyFile (only for backward compatibility)
139    
         _BEGIN_MASTER( myThid )  
140          CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )          CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
         _END_MASTER(myThid)  
         _BARRIER  
141    
142        ELSEIF ( topoFile.EQ.' ' ) THEN        ELSEIF ( topoFile.EQ.' ' ) THEN
143  C------ set default value:  C------ set default value:
# Line 150  C------ set default value: Line 156  C------ set default value:
156  C------ read from file:  C------ read from file:
157    
158  C- read surface topography (in m) from topoFile (case topoFile.NE.' '):  C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
159          _BEGIN_MASTER( myThid )  
160          CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )          CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
         _END_MASTER(myThid)  
161          _BARRIER          _BARRIER
162    
163          IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN          IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
164  C----  C----
165  C   Convert Surface Geopotential to (reference) Surface Pressure  C   Convert Surface Geopotential to (reference) Surface Pressure
166  C   according to Tref profile, using same discretisation as in calc_phi_hyd  C   according to Tref profile, using same discretisation as in calc_phi_hyd
167  C----  C----
168  c         _BEGIN_MASTER( myThid )  c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
 c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)  
 c         _END_MASTER(myThid)  
169    
170            CALL INI_P_GROUND( 2, topoZ,            CALL INI_P_GROUND( 2, topoZ,
171       O                       Ro_surf,       O                       Ro_surf,
172       I                       myThid )       I                       myThid )
173    
174            _BARRIER  c         _BARRIER
175            _BEGIN_MASTER( myThid )  C         This I/O is now done in write_grid.F
176            CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)  c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
           _END_MASTER(myThid)  
177    
178          ELSE          ELSE
179  C----  C----
180  C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary  C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
181  C    below an ice-shelf - NOTE - actually not yet implemented )  C    below an ice-shelf - NOTE - actually not yet implemented )
182            DO bj = myByLo(myThid), myByHi(myThid)            DO bj = myByLo(myThid), myByHi(myThid)
183             DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 192  C    below an ice-shelf - NOTE - actuall Line 194  C    below an ice-shelf - NOTE - actuall
194  C------ end case "read topoFile"  C------ end case "read topoFile"
195        ENDIF        ENDIF
196    
197  C----- fill in the overlap :  C----- fill in the overlap (+ BARRIER):
198        _EXCH_XY_R4(Ro_surf, myThid )        _EXCH_XY_R4(Ro_surf, myThid )
199    
200  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 200  c---+----1----+----2----+----3----+----4 Line 202  c---+----1----+----2----+----3----+----4
202  C------  C------
203  C   3) Close the Domain (special configuration).  C   3) Close the Domain (special configuration).
204  C------  C------
205        IF (groundAtK1) THEN        IF (usingPCoords) THEN
206         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
207          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
208           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 213  C          Test for northern edge Line 215  C          Test for northern edge
215  c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.  c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
216  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
217             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
218       &         abs(yC(I,J,bi,bj)).GE.-phiMin)       &         ABS(yC(I,J,bi,bj)).GE.-phiMin)
219       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
220             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. ABS(yC(I,J,bi,bj)).GE.90. )
221       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
222            ENDDO            ENDDO
223           ENDDO           ENDDO
# Line 234  C          Test for northern edge Line 236  C          Test for northern edge
236  c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.  c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
237  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
238             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
239       &         abs(yC(I,J,bi,bj)).GE.-phiMin)       &         ABS(yC(I,J,bi,bj)).GE.-phiMin)
240       &       R_low(I,J,bi,bj) = Ro_SeaLevel       &       R_low(I,J,bi,bj) = Ro_SeaLevel
241             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. ABS(yC(I,J,bi,bj)).GE.90. )
242       &       R_low(I,J,bi,bj) = Ro_SeaLevel       &       R_low(I,J,bi,bj) = Ro_SeaLevel
243            ENDDO            ENDDO
244           ENDDO           ENDDO
# Line 244  c- Domain : Symetric % Eq. & closed at N Line 246  c- Domain : Symetric % Eq. & closed at N
246         ENDDO         ENDDO
247        ENDIF        ENDIF
248    
249  c     _BEGIN_MASTER( myThid )        IF ( debugLevel.GE.debLevB ) THEN
250  c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)         CALL PLOT_FIELD_XYRS( Ro_surf,
251  c     _END_MASTER(myThid)       &           'Surface reference r-position (ini_depths)',
252         &           -1, myThid )
253          ENDIF
254    c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
255    
256    #ifdef ALLOW_SHELFICE
257    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
258          IF ( useShelfIce ) THEN
259    C------
260    C   4) Set R_shelfIce = the Lower (in r sense) boundary of floating shelfice :
261    C------
262          IF (usingPCoords .OR. shelfIceFile .EQ. ' ') THEN
263    C- e.g., atmosphere : R_low = Top of atmosphere
264    C-            ocean : R_low = Bottom
265           DO bj = myByLo(myThid), myByHi(myThid)
266            DO bi = myBxLo(myThid), myBxHi(myThid)
267             DO j=1,sNy
268              DO i=1,sNx
269               R_shelfIce(i,j,bi,bj) = 0. _d 0
270              ENDDO
271             ENDDO
272            ENDDO
273           ENDDO
274          ELSE
275    C Read the shelfIce draught using the mid-level I/O pacakage read_write_rec
276    C The 0 is the "iteration" argument. The 1 is the record number.
277            CALL READ_REC_XY_RS( shelfIceFile, R_shelfIce, 1, 0, myThid )
278    C Read the shelfIce draught using the mid-level I/O pacakage read_write_fld
279    C The 0 is the "iteration" argument. The ' ' is an empty suffix
280    C        CALL READ_FLD_XY_RS( shelfIceFile, ' ', R_shelfIce, 0, myThid )
281    c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
282    C Read the selfIce draught using the low-level I/O package
283    c       CALL MDSREADFIELD( shelfIceFile, readBinaryPrec,
284    c    &                     'RS', 1, R_selfIce, 1, myThid )
285    
286          ENDIF
287    C- end setup R_shelfIce in the interior
288    
289    C- fill in the overlap (+ BARRIER):
290          _EXCH_XY_R4(R_shelfIce, myThid )
291    
292    c     CALL PLOT_FIELD_XYRS(R_selfIce,'Shelf ice draught (ini_depths)',
293    c    &     1,myThid)
294    CML      CALL WRITE_FLD_XY_RS( 'R_shelfIce' ,' ', R_shelfIce, 0,myThid)
295          ENDIF
296    #endif /* ALLOW_SHELFICE */
297    
298    C--   Everyone else must wait for the depth to be loaded
299    C-    note: might not be necessary since all single-thread IO above
300    C           are followed by an EXCH (with BARRIER)
301          _BARRIER
302    
303        RETURN        RETURN
304        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22