/[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.9 by adcroft, Thu Jul 2 15:46:21 1998 UTC revision 1.26 by cnh, Wed Sep 26 18:09:15 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7    C     !ROUTINE: INI_DEPTHS
8    C     !INTERFACE:
9        SUBROUTINE INI_DEPTHS( myThid )        SUBROUTINE INI_DEPTHS( myThid )
10  C     /==========================================================\  C     !DESCRIPTION: \bv
11  C     | SUBROUTINE INI_DEPTHS                                    |  C     *==========================================================*
12  C     | o Initialise map of model depths                         |  C     | SUBROUTINE INI_DEPTHS                                    
13  C     |==========================================================|  C     | o define R_position of Lower and Surface Boundaries      
14  C     | The depths of the bottom of the model is specified in    |  C     *==========================================================*
15  C     | terms of an XY map with one depth for each column of     |  C     |atmosphere orography:                                      
16  C     | grid cells. Depths do not have to coincide with the      |  C     | define either in term of P_topo or converted from Z_topo  
17  C     | model levels. The model's lopping algorithm makes it     |  C     |ocean bathymetry:                                          
18  C     | possible to represent arbitrary depths.                  |  C     | The depths of the bottom of the model is specified in    
19  C     | The mode depths map also influences the models topology  |  C     | terms of an XY map with one depth for each column of      
20  C     | By default the model domain wraps around in X and Y.     |  C     | grid cells. Depths do not have to coincide with the      
21  C     | This default doubly periodic topology is "supressed"     |  C     | model levels. The model lopping algorithm makes it        
22  C     | if a depth map is defined which closes off all wrap      |  C     | possible to represent arbitrary depths.                  
23  C     | around flow.                                             |  C     | The mode depths map also influences the models topology  
24  C     \==========================================================/  C     | By default the model domain wraps around in X and Y.      
25    C     | This default doubly periodic topology is "supressed"      
26    C     | if a depth map is defined which closes off all wrap      
27    C     | around flow.                                              
28    C     *==========================================================*
29    C     \ev
30    
31    C     !USES:
32          IMPLICIT NONE
33  C     === Global variables ===  C     === Global variables ===
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "GRID.h"  #include "GRID.h"
38    
39    C     !INPUT/OUTPUT PARAMETERS:
40  C     == Routine arguments ==  C     == Routine arguments ==
41  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid -  Number of this instance of INI_DEPTHS
42        INTEGER myThid        INTEGER myThid
43  CEndOfInterface  CEndOfInterface
44    
45    C     !LOCAL VARIABLES:
46    C     == Local variables in common ==
47    C     Hloc  - Temporary array used to read surface topography
48    C             has to be in common for multi threading    
49          COMMON / LOCAL_INI_DEPTHS / Hloc
50          _RS Hloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
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  - Loop counters
54  C     I,J,K  C     I,J,K
55  C     phi    - total depth of model  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 phi        CHARACTER*(MAX_LEN_MBUF) msgBuf
61    CEOP
62    
63          IF (groundAtK1 .AND. bathyFile .NE. ' '
64         &               .AND. topoFile  .NE. ' ' ) THEN
65           WRITE(msgBuf,'(A,A)')
66         &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
67         &  ' select the right one !'
68           CALL PRINT_ERROR( msgBuf , myThid)
69           STOP 'ABNORMAL END: S/R INI_DEPTHS'
70          ENDIF
71    
72    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             ENDDO
82            ENDDO
83           ENDDO
84          ENDDO
85    
86        _BARRIER  C------
87        IF ( bathyFile .EQ. ' '  ) THEN  C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
88  C      Set up a flat bottom box with doubly periodic topology.  C------
89  C      H is the basic variable from which other terms are derived. It        IF (groundAtK1 .OR. bathyFile .EQ. ' ') THEN
90  C      is the term that would be set from an external file for a  C- e.g., atmosphere : R_low = Top of atmosphere
91  C      realistic problem.  C-            ocean : R_low = Bottom
        phi = zFace(nZ+1)  
92         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
93          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
94           DO j=1,sNy           DO j=1,sNy
95            DO i=1,sNx            DO i=1,sNx
96               R_low(i,j,bi,bj) = rF(Nr+1)
97              ENDDO
98             ENDDO
99            ENDDO
100           ENDDO
101          ELSE
102            _BEGIN_MASTER( myThid )
103    C Read the bathymetry using the mid-level I/O pacakage read_write_rec
104    C The 0 is the "iteration" argument. The 1 is the record number.
105            CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
106    C Read the bathymetry using the mid-level I/O pacakage read_write_fld
107    C The 0 is the "iteration" argument. The ' ' is an empty suffix
108    c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
109    C Read the bathymetry using the low-level I/O package
110    c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
111    c    &                     'RS', 1, R_low, 1, myThid )
112            _END_MASTER(myThid)
113    
114          ENDIF
115    C- end setup R_low in the interior
116    
117    C- fill in the overlap :
118          _EXCH_XY_R4(R_low, myThid )
119    
120    c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
121    c     _BEGIN_MASTER( myThid )
122    c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
123    c     _END_MASTER(myThid)
124    
125    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
126    
127    C------
128    C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
129    C------
130    
131          IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN
132    C------ read directly Po_surf from bathyFile (only for backward compatibility)
133    
134            _BEGIN_MASTER( myThid )
135            CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
136            _END_MASTER(myThid)
137            _BARRIER
138    
139          ELSEIF ( topoFile.EQ.' ' ) THEN
140    C------ set default value:
141    
142            DO bj = myByLo(myThid), myByHi(myThid)
143             DO bi = myBxLo(myThid), myBxHi(myThid)
144              DO j=1,sNy
145               DO i=1,sNx
146                Ro_surf(i,j,bi,bj) = Ro_SeaLevel
147               ENDDO
148              ENDDO
149             ENDDO
150            ENDDO
151    
152          ELSE
153    C------ read from file:
154    
155    C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
156            _BEGIN_MASTER( myThid )
157            CALL READ_REC_XY_RS( topoFile, Hloc, 1, 0, myThid )
158            _END_MASTER(myThid)
159            _BARRIER
160    
161            IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
162    C----
163    C   Convert Surface Geopotential to (reference) Surface Pressure
164    C   according to Tref profile, using same discretisation as in calc_phi_hyd
165    C----
166    c         _BEGIN_MASTER( myThid )
167    c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',Hloc,0,myThid)
168    c         _END_MASTER(myThid)
169    
170              CALL INI_P_GROUND( Hloc, Ro_surf, myThid )
171    
172              _BARRIER
173              _BEGIN_MASTER( myThid )
174              CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
175              _END_MASTER(myThid)
176    
177            ELSE
178    C----
179    C   Direct Transfer to Ro_surf :
180              DO bj = myByLo(myThid), myByHi(myThid)
181               DO bi = myBxLo(myThid), myBxHi(myThid)
182                DO j=1,sNy
183                 DO i=1,sNx
184                  Ro_surf(i,j,bi,bj) = Hloc(i,j,bi,bj)
185                 ENDDO
186                ENDDO
187               ENDDO
188              ENDDO
189    
190            ENDIF
191    
192    C------ end case "read topoFile"
193          ENDIF
194    
195    C----- fill in the overlap :
196          _EXCH_XY_R4(Ro_surf, myThid )
197    
198    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
199    
200    C------
201    C   3) Close the Domain (special configuration).
202    C------
203          IF (groundAtK1) THEN
204           DO bj = myByLo(myThid), myByHi(myThid)
205            DO bi = myBxLo(myThid), myBxHi(myThid)
206             DO j=1-Oly,sNy+Oly
207              DO i=1-Olx,sNx+Olx
208             iG = myXGlobalLo-1+(bi-1)*sNx+I             iG = myXGlobalLo-1+(bi-1)*sNx+I
209             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  
210  C          Test for eastern edge  C          Test for eastern edge
211             IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.  c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
212  C          Test for northern edge  C          Test for northern edge
213             IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.  c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
214  C          Island             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
215             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  
216            ENDDO            ENDDO
217           ENDDO           ENDDO
218          ENDDO          ENDDO
219         ENDDO         ENDDO
220        ELSE        ELSE
221         _BEGIN_MASTER( myThid )         DO bj = myByLo(myThid), myByHi(myThid)
222         CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )          DO bi = myBxLo(myThid), myBxHi(myThid)
223         _END_MASTER(myThid)           DO j=1-Oly,sNy+Oly
224              DO i=1-Olx,sNx+Olx
225               iG = myXGlobalLo-1+(bi-1)*sNx+I
226               jG = myYGlobalLo-1+(bj-1)*sNy+J
227    C          Test for eastern edge
228    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
229    C          Test for northern edge
230    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
231               IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
232         &       R_low(I,J,bi,bj) = Ro_SeaLevel
233              ENDDO
234             ENDDO
235            ENDDO
236           ENDDO
237        ENDIF        ENDIF
238        _EXCH_XY_R4(    H, myThid )  
239  C  c     _BEGIN_MASTER( myThid )
240        CALL PLOT_FIELD_XYRS( H, 'Model depths' , 1, myThid )  c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
241  C  c     _END_MASTER(myThid)
242    
243        RETURN        RETURN
244        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22