/[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.17 by adcroft, Wed Dec 9 16:11:52 1998 UTC revision 1.25 by jmc, Mon Jul 9 16:46:29 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 6  CStartOfInterface Line 7  CStartOfInterface
7        SUBROUTINE INI_DEPTHS( myThid )        SUBROUTINE INI_DEPTHS( myThid )
8  C     /==========================================================\  C     /==========================================================\
9  C     | SUBROUTINE INI_DEPTHS                                    |  C     | SUBROUTINE INI_DEPTHS                                    |
10  C     | o Initialise map of model depths                         |  C     | o define R_position of Lower and Surface Boundaries      |
11  C     |==========================================================|  C     |==========================================================|
12    C     |atmosphere orography:                                     |
13    C     | define either in term of P_topo or converted from Z_topo |
14    C     |ocean bathymetry:                                         |
15  C     | The depths of the bottom of the model is specified in    |  C     | The depths of the bottom of the model is specified in    |
16  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     |
17  C     | grid cells. Depths do not have to coincide with the      |  C     | grid cells. Depths do not have to coincide with the      |
# Line 32  C     myThid -  Number of this instance Line 36  C     myThid -  Number of this instance
36        INTEGER myThid        INTEGER myThid
37  CEndOfInterface  CEndOfInterface
38    
39    C     == Local variables in common ==
40    C     Hloc  - Temporary array used to read surface topography
41    C             has to be in common for multi threading    
42          COMMON / LOCAL_INI_DEPTHS / Hloc
43          _RS Hloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44    
45  C     == Local variables ==  C     == Local variables ==
46  C     iG, jG - Global coordinate index  C     iG, jG - Global coordinate index
47  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
48  C     I,J,K  C     I,J,K
 C     phi    - total depth of model  
49  C     oldPrec - Temporary used in controlling binary input dataset precision  C     oldPrec - Temporary used in controlling binary input dataset precision
50    C     msgBuf    - Informational/error meesage buffer
51        INTEGER iG, jG        INTEGER iG, jG
52        INTEGER bi, bj        INTEGER bi, bj
53        INTEGER  I,  J, K        INTEGER  I, J
54        INTEGER oldPrec        CHARACTER*(MAX_LEN_MBUF) msgBuf
55        _RL phi  
56          IF (groundAtK1 .AND. bathyFile .NE. ' '
57        _BARRIER       &               .AND. topoFile  .NE. ' ' ) THEN
58        IF ( bathyFile .EQ. ' '  ) THEN         WRITE(msgBuf,'(A,A)')
59  C      Set up a flat bottom box with doubly periodic topology.       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
60  C      H is the basic variable from which other terms are derived. It       &  ' select the right one !'
61  C      is the term that would be set from an external file for a         CALL PRINT_ERROR( msgBuf , myThid)
62  C      realistic problem.         STOP 'ABNORMAL END: S/R INI_DEPTHS'
63         phi = rF(Nr+1)        ENDIF
64    
65    C------
66    C   0) Initialize R_low and Ro_surf (define an empty domain)
67    C------
68          DO bj = myByLo(myThid), myByHi(myThid)
69           DO bi = myBxLo(myThid), myBxHi(myThid)
70            DO j=1-Oly,sNy+Oly
71             DO i=1-Olx,sNx+Olx
72              R_low(i,j,bi,bj) = 0.
73              Ro_surf(i,j,bi,bj) = 0.
74             ENDDO
75            ENDDO
76           ENDDO
77          ENDDO
78    
79    C------
80    C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
81    C------
82          IF (groundAtK1 .OR. bathyFile .EQ. ' ') THEN
83    C- e.g., atmosphere : R_low = Top of atmosphere
84    C-            ocean : R_low = Bottom
85         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
86          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
87           DO j=1,sNy           DO j=1,sNy
88            DO i=1,sNx            DO i=1,sNx
89               R_low(i,j,bi,bj) = rF(Nr+1)
90              ENDDO
91             ENDDO
92            ENDDO
93           ENDDO
94          ELSE
95            _BEGIN_MASTER( myThid )
96    C Read the bathymetry using the mid-level I/O pacakage read_write_rec
97    C The 0 is the "iteration" argument. The 1 is the record number.
98            CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
99    C Read the bathymetry using the mid-level I/O pacakage read_write_fld
100    C The 0 is the "iteration" argument. The ' ' is an empty suffix
101    c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
102    C Read the bathymetry using the low-level I/O package
103    c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
104    c    &                     'RS', 1, R_low, 1, myThid )
105            _END_MASTER(myThid)
106    
107          ENDIF
108    C- end setup R_low in the interior
109    
110    C- fill in the overlap :
111          _EXCH_XY_R4(R_low, myThid )
112    
113    c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
114    c     _BEGIN_MASTER( myThid )
115    c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
116    c     _END_MASTER(myThid)
117    
118    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
119    
120    C------
121    C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
122    C------
123    
124          IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN
125    C------ read directly Po_surf from bathyFile (only for backward compatibility)
126    
127            _BEGIN_MASTER( myThid )
128            CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
129            _END_MASTER(myThid)
130            _BARRIER
131    
132          ELSEIF ( topoFile.EQ.' ' ) THEN
133    C------ set default value:
134    
135            DO bj = myByLo(myThid), myByHi(myThid)
136             DO bi = myBxLo(myThid), myBxHi(myThid)
137              DO j=1,sNy
138               DO i=1,sNx
139                Ro_surf(i,j,bi,bj) = Ro_SeaLevel
140               ENDDO
141              ENDDO
142             ENDDO
143            ENDDO
144    
145          ELSE
146    C------ read from file:
147    
148    C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
149            _BEGIN_MASTER( myThid )
150            CALL READ_REC_XY_RS( topoFile, Hloc, 1, 0, myThid )
151            _END_MASTER(myThid)
152            _BARRIER
153    
154            IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
155    C----
156    C   Convert Surface Geopotential to (reference) Surface Pressure
157    C   according to Tref profile, using same discretisation as in calc_phi_hyd
158    C----
159    c         _BEGIN_MASTER( myThid )
160    c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',Hloc,0,myThid)
161    c         _END_MASTER(myThid)
162    
163              CALL INI_P_GROUND( Hloc, Ro_surf, myThid )
164    
165              _BARRIER
166              _BEGIN_MASTER( myThid )
167              CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
168              _END_MASTER(myThid)
169    
170            ELSE
171    C----
172    C   Direct Transfer to Ro_surf :
173              DO bj = myByLo(myThid), myByHi(myThid)
174               DO bi = myBxLo(myThid), myBxHi(myThid)
175                DO j=1,sNy
176                 DO i=1,sNx
177                  Ro_surf(i,j,bi,bj) = Hloc(i,j,bi,bj)
178                 ENDDO
179                ENDDO
180               ENDDO
181              ENDDO
182    
183            ENDIF
184    
185    C------ end case "read topoFile"
186          ENDIF
187    
188    C----- fill in the overlap :
189          _EXCH_XY_R4(Ro_surf, myThid )
190    
191    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192    
193    C------
194    C   3) Close the Domain (special configuration).
195    C------
196          IF (groundAtK1) THEN
197           DO bj = myByLo(myThid), myByHi(myThid)
198            DO bi = myBxLo(myThid), myBxHi(myThid)
199             DO j=1-Oly,sNy+Oly
200              DO i=1-Olx,sNx+Olx
201             iG = myXGlobalLo-1+(bi-1)*sNx+I             iG = myXGlobalLo-1+(bi-1)*sNx+I
202             jG = myYGlobalLo-1+(bj-1)*sNy+J             jG = myYGlobalLo-1+(bj-1)*sNy+J
 C          Default depth of full domain  
            H(i,j,bi,bj) = phi  
203  C          Test for eastern edge  C          Test for eastern edge
204             IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.  c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
205  C          Test for northern edge  C          Test for northern edge
206             IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.  c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
207  C          Island             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
208             IF ( iG .EQ. 1 .AND.       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
      &          jG .EQ. 24 ) H(i,j,bi,bj) = 0.75*phi  
209            ENDDO            ENDDO
210           ENDDO           ENDDO
211          ENDDO          ENDDO
212         ENDDO         ENDDO
213        ELSE        ELSE
214         _BEGIN_MASTER( myThid )         DO bj = myByLo(myThid), myByHi(myThid)
215            DO bi = myBxLo(myThid), myBxHi(myThid)
216             DO j=1-Oly,sNy+Oly
217              DO i=1-Olx,sNx+Olx
218               iG = myXGlobalLo-1+(bi-1)*sNx+I
219               jG = myYGlobalLo-1+(bj-1)*sNy+J
220    C          Test for eastern edge
221    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
222    C          Test for northern edge
223    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
224               IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
225         &       R_low(I,J,bi,bj) = Ro_SeaLevel
226              ENDDO
227             ENDDO
228            ENDDO
229           ENDDO
230          ENDIF
231    
232  CcnhDebugStarts  c     _BEGIN_MASTER( myThid )
233  C       Force 64-bit IO  c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
234          oldPrec        = readBinaryPrec  c     _END_MASTER(myThid)
         readBinaryPrec = precFloat64  
 CcnhDEbugEnds  
        CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  
 CcnhDebugStarts  
        readBinaryPrec = oldPrec  
 CcnhdDebugEnds  
235    
        _END_MASTER(myThid)  
       ENDIF  
         
       _EXCH_XY_R4(    H, myThid )  
 C  
       CALL PLOT_FIELD_XYRS( H, 'Model depths' , 1, myThid )  
 C  
236        RETURN        RETURN
237        END        END

Legend:
Removed from v.1.17  
changed lines
  Added in v.1.25

  ViewVC Help
Powered by ViewVC 1.1.22