/[MITgcm]/MITgcm/model/src/ini_depths.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_depths.F

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

revision 1.22 by cnh, Sun Feb 4 14:38:47 2001 UTC revision 1.35 by mlosch, Tue Feb 7 11:47:48 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  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_DEPTHS
9    C     !INTERFACE:
10        SUBROUTINE INI_DEPTHS( myThid )        SUBROUTINE INI_DEPTHS( myThid )
11  C     /==========================================================\  C     !DESCRIPTION: \bv
12  C     | SUBROUTINE INI_DEPTHS                                    |  C     *==========================================================*
13  C     | o Initialise map of model depths                         |  C     | SUBROUTINE INI_DEPTHS
14  C     |==========================================================|  C     | o define R_position of Lower and Surface Boundaries
15  C     | The depths of the bottom of the model is specified in    |  C     *==========================================================*
16  C     | terms of an XY map with one depth for each column of     |  C     |atmosphere orography:
17  C     | grid cells. Depths do not have to coincide with the      |  C     | define either in term of P_topo or converted from Z_topo
18  C     | model levels. The model lopping algorithm makes it       |  C     |ocean bathymetry:
19  C     | possible to represent arbitrary depths.                  |  C     | The depths of the bottom of the model is specified in
20  C     | The mode depths map also influences the models topology  |  C     | terms of an XY map with one depth for each column of
21  C     | By default the model domain wraps around in X and Y.     |  C     | grid cells. Depths do not have to coincide with the
22  C     | This default doubly periodic topology is "supressed"     |  C     | model levels. The model lopping algorithm makes it
23  C     | if a depth map is defined which closes off all wrap      |  C     | possible to represent arbitrary depths.
24  C     | around flow.                                             |  C     | The mode depths map also influences the models topology
25  C     \==========================================================/  C     | By default the model domain wraps around in X and Y.
26        IMPLICIT NONE  C     | This default doubly periodic topology is "supressed"
27    C     | if a depth map is defined which closes off all wrap
28    C     | around flow.
29    C     *==========================================================*
30    C     \ev
31    
32    C     !USES:
33          IMPLICIT NONE
34  C     === Global variables ===  C     === Global variables ===
35  #include "SIZE.h"  #include "SIZE.h"
36  #include "EEPARAMS.h"  #include "EEPARAMS.h"
37  #include "PARAMS.h"  #include "PARAMS.h"
38  #include "GRID.h"  #include "GRID.h"
39    #include "SURFACE.h"
40    #ifdef ALLOW_SHELFICE
41    # include "SHELFICE.h"
42    #endif /* ALLOW_SHELFICE */
43    
44    C     !INPUT/OUTPUT PARAMETERS:
45  C     == Routine arguments ==  C     == Routine arguments ==
46  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid -  Number of this instance of INI_DEPTHS
47        INTEGER myThid        INTEGER myThid
48  CEndOfInterface  CEndOfInterface
49    
50    C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52  C     iG, jG - Global coordinate index  C     iG, jG - Global coordinate index
53  C     bi,bj  - Loop counters  C     bi,bj  - Tile indices
54  C     I,J,K  C     I,J,K  - Loop counters
 C     Hdefault - default r-coordinate of the lower boundary (=ground)  
 C                (=minus(Total_depth) in the ocean model)  
 C                (=Total Pressure at Sea Level in the atmos model)  
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
57        INTEGER iG, jG        INTEGER iG, jG
58        INTEGER bi, bj        INTEGER bi, bj
59        INTEGER  I, J, K        INTEGER  I, J
60        _RL Hdefault        CHARACTER*(MAX_LEN_MBUF) msgBuf
61    CEOP
62    
63        _BARRIER        IF (usingPCoords .AND. bathyFile .NE. ' '
64        IF ( bathyFile .EQ. ' '  ) THEN       &                 .AND. topoFile  .NE. ' ' ) THEN
65  C      Set up a flat bottom box with doubly periodic topology.         WRITE(msgBuf,'(A,A)')
66  C      H is the basic variable from which other terms are derived. It       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
67  C      is the term that would be set from an external file for a       &  ' select the right one !'
68  C      realistic problem.         CALL PRINT_ERROR( msgBuf , myThid)
69         Hdefault = rF(1)         STOP 'ABNORMAL END: S/R INI_DEPTHS'
70         IF (.NOT. groundAtK1) THEN        ENDIF
71          DO K = 1, Nr  
72           Hdefault = Hdefault - rkfac*drF(K)  C------
73    C   0) Initialize R_low and Ro_surf (define an empty domain)
74    C------
75          DO bj = myByLo(myThid), myByHi(myThid)
76           DO bi = myBxLo(myThid), myBxHi(myThid)
77            DO j=1-Oly,sNy+Oly
78             DO i=1-Olx,sNx+Olx
79              R_low(i,j,bi,bj) = 0.
80              Ro_surf(i,j,bi,bj) = 0.
81              topoZ(i,j,bi,bj) = 0.
82             ENDDO
83          ENDDO          ENDDO
84         ENDIF         ENDDO
85          ENDDO
86    
87    C------
88    C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
89    C------
90          IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
91    C- e.g., atmosphere : R_low = Top of atmosphere
92    C-            ocean : R_low = Bottom
93         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
94          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
95           DO j=1,sNy           DO j=1,sNy
96            DO i=1,sNx            DO i=1,sNx
97             iG = myXGlobalLo-1+(bi-1)*sNx+I             R_low(i,j,bi,bj) = rF(Nr+1)
            jG = myYGlobalLo-1+(bj-1)*sNy+J  
 C          Default depth of full domain  
            H(i,j,bi,bj) = Hdefault  
 C          Test for eastern edge  
 C          IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.  
 C          Test for northern edge  
 C          IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.  
98            ENDDO            ENDDO
99           ENDDO           ENDDO
100          ENDDO          ENDDO
101         ENDDO         ENDDO
102        ELSE        ELSE
103         _BEGIN_MASTER( myThid )          _BARRIER
104    C       _BEGIN_MASTER( myThid )
105  C Read the bathymetry using the mid-level I/O pacakage read_write_rec  C Read the bathymetry using the mid-level I/O pacakage read_write_rec
106  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.
107         CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
108  C Read the bathymetry using the mid-level I/O pacakage read_write_fld  C Read the bathymetry using the mid-level I/O pacakage read_write_fld
109  C The 0 is the "iteration" argument. The ' ' is an empty suffix  C The 0 is the "iteration" argument. The ' ' is an empty suffix
110  C      CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
111  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
112  C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
113  C    &                    'RS', 1, H, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
114         _END_MASTER(myThid)  C       _END_MASTER(myThid)
115            _BARRIER
116        ENDIF        ENDIF
117    C- end setup R_low in the interior
118    
119    C- fill in the overlap :
120          _EXCH_XY_R4(R_low, myThid )
121    
122    c      PRINT *, ' Calling plot field', myThid
123          CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
124    c     _BEGIN_MASTER( myThid )
125    c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
126    c     _END_MASTER(myThid)
127    
128    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
129    
130        _EXCH_XY_R4(    H, myThid )  C------
131    C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
132    C------
133    
134        IF (usingSphericalPolarGrid) THEN        IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
135    C------ read directly Po_surf from bathyFile (only for backward compatibility)
136    
137            _BEGIN_MASTER( myThid )
138            CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
139            _END_MASTER(myThid)
140            _BARRIER
141    
142          ELSEIF ( topoFile.EQ.' ' ) THEN
143    C------ set default value:
144    
145            DO bj = myByLo(myThid), myByHi(myThid)
146             DO bi = myBxLo(myThid), myBxHi(myThid)
147              DO j=1,sNy
148               DO i=1,sNx
149                Ro_surf(i,j,bi,bj) = Ro_SeaLevel
150               ENDDO
151              ENDDO
152             ENDDO
153            ENDDO
154    
155          ELSE
156    C------ read from file:
157    
158    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 )
161            _END_MASTER(myThid)
162            _BARRIER
163    
164            IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
165    C----
166    C   Convert Surface Geopotential to (reference) Surface Pressure
167    C   according to Tref profile, using same discretisation as in calc_phi_hyd
168    C----
169    c         _BEGIN_MASTER( myThid )
170    c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
171    c         _END_MASTER(myThid)
172    
173              CALL INI_P_GROUND( 2, topoZ,
174         O                       Ro_surf,
175         I                       myThid )
176    
177              _BARRIER
178    C         This I/O is now done in write_grid.F
179    c         _BEGIN_MASTER( myThid )
180    c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
181    c         _END_MASTER(myThid)
182    
183            ELSE
184    C----
185    C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
186    C    below an ice-shelf - NOTE - actually not yet implemented )
187              DO bj = myByLo(myThid), myByHi(myThid)
188               DO bi = myBxLo(myThid), myBxHi(myThid)
189                DO j=1,sNy
190                 DO i=1,sNx
191                  Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
192                 ENDDO
193                ENDDO
194               ENDDO
195              ENDDO
196    
197            ENDIF
198    
199    C------ end case "read topoFile"
200          ENDIF
201    
202    C----- fill in the overlap :
203          _EXCH_XY_R4(Ro_surf, myThid )
204    
205    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
206    
207    C------
208    C   3) Close the Domain (special configuration).
209    C------
210          IF (usingPCoords) THEN
211           DO bj = myByLo(myThid), myByHi(myThid)
212            DO bi = myBxLo(myThid), myBxHi(myThid)
213             DO j=1-Oly,sNy+Oly
214              DO i=1-Olx,sNx+Olx
215               iG = myXGlobalLo-1+(bi-1)*sNx+I
216               jG = myYGlobalLo-1+(bj-1)*sNy+J
217    C          Test for eastern edge
218    c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
219    C          Test for northern edge
220    c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
221               IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
222         &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
223              ENDDO
224             ENDDO
225            ENDDO
226           ENDDO
227          ELSE
228         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
229          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
230           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
231            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
232             IF (abs(yC(I,J,bi,bj)).GE.90.) H(I,J,bi,bj)=0.             iG = myXGlobalLo-1+(bi-1)*sNx+I
233               jG = myYGlobalLo-1+(bj-1)*sNy+J
234    C          Test for eastern edge
235    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
236    C          Test for northern edge
237    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
238               IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
239         &       R_low(I,J,bi,bj) = Ro_SeaLevel
240            ENDDO            ENDDO
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244        ENDIF        ENDIF
245  C  
246        CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)  c     _BEGIN_MASTER( myThid )
247  C  c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
248    c     _END_MASTER(myThid)
249    
250    #ifdef ALLOW_SHELFICE
251    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
252          IF ( useShelfIce ) THEN
253    C------
254    C   4) Set R_shelfIce = the Lower (in r sense) boundary of floating shelfice :
255    C------
256          IF (usingPCoords .OR. shelfIceFile .EQ. ' ') THEN
257    C- e.g., atmosphere : R_low = Top of atmosphere
258    C-            ocean : R_low = Bottom
259           DO bj = myByLo(myThid), myByHi(myThid)
260            DO bi = myBxLo(myThid), myBxHi(myThid)
261             DO j=1,sNy
262              DO i=1,sNx
263               R_shelfIce(i,j,bi,bj) = 0. _d 0
264              ENDDO
265             ENDDO
266            ENDDO
267           ENDDO
268          ELSE
269            _BEGIN_MASTER( myThid )
270    C Read the shelfIce draught using the mid-level I/O pacakage read_write_rec
271    C The 0 is the "iteration" argument. The 1 is the record number.
272            CALL READ_REC_XY_RS( shelfIceFile, R_shelfIce, 1, 0, myThid )
273    C Read the shelfIce draught using the mid-level I/O pacakage read_write_fld
274    C The 0 is the "iteration" argument. The ' ' is an empty suffix
275    C        CALL READ_FLD_XY_RS( shelfIceFile, ' ', R_shelfIce, 0, myThid )
276    c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
277    C Read the selfIce draught using the low-level I/O package
278    c       CALL MDSREADFIELD( shelfIceFile, readBinaryPrec,
279    c    &                     'RS', 1, R_selfIce, 1, myThid )
280            _END_MASTER(myThid)
281    
282          ENDIF
283    C- end setup R_shelfIce in the interior
284    
285    C- fill in the overlap :
286          _EXCH_XY_R4(R_shelfIce, myThid )
287    
288    c     CALL PLOT_FIELD_XYRS(R_selfIce,'Shelf ice draught (ini_depths)',
289    c    &     1,myThid)
290    CML      _BEGIN_MASTER( myThid )
291    CML      CALL WRITE_FLD_XY_RS( 'R_shelfIce' ,' ', R_shelfIce, 0,myThid)
292    CML      _END_MASTER(myThid)
293          ENDIF
294    #endif /* ALLOW_SHELFICE */
295    
296        RETURN        RETURN
297        END        END
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.22