/[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.4 by adcroft, Tue May 29 14:01:46 2001 UTC revision 1.5 by jmc, Fri Aug 24 11:57:43 2001 UTC
# Line 7  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 27  C     === Global variables === Line 30  C     === Global variables ===
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
 #include "INI_DEPTHS.h"  
33    
34  C     == Routine arguments ==  C     == Routine arguments ==
35  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid -  Number of this instance of INI_DEPTHS
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     H      - Depth of base of fluid from upper surface f[X,Y] (m).  
 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)  
 C--------------------  
 C NOTE: will change soon: 2 separed files for r_lower and r_surface boudaries  
 C        and for the atmosphere, topography will be defined in term of height  
 C--------------------  
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        INTEGER  I, J
54        _RL Hdefault        CHARACTER*(MAX_LEN_MBUF) msgBuf
55    
56          IF (groundAtK1 .AND. bathyFile .NE. ' '
57         &               .AND. topoFile  .NE. ' ' ) THEN
58           WRITE(msgBuf,'(A,A)')
59         &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
60         &  ' select the right one !'
61           CALL PRINT_ERROR( msgBuf , myThid)
62           STOP 'ABNORMAL END: S/R INI_DEPTHS'
63          ENDIF
64    
65        _BARRIER  C------
66        IF ( bathyFile .EQ. ' '  ) THEN  C   0) Initialize R_low and Ro_surf (define an empty domain)
67  C      Set up a flat bottom box with doubly periodic topology.  C------
68  C      H is the basic variable from which other terms are derived. It        DO bj = myByLo(myThid), myByHi(myThid)
69  C      is the term that would be set from an external file for a         DO bi = myBxLo(myThid), myBxHi(myThid)
70  C      realistic problem.          DO j=1-Oly,sNy+Oly
71         IF (groundAtK1) THEN           DO i=1-Olx,sNx+Olx
72           Hdefault = Ro_SeaLevel            R_low(i,j,bi,bj) = 0.
73         ELSE            Ro_surf(i,j,bi,bj) = 0.
74           Hdefault = rF(Nr+1)           ENDDO
75         ENDIF          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             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.  
90            ENDDO            ENDDO
91           ENDDO           ENDDO
92          ENDDO          ENDDO
93         ENDDO         ENDDO
94        ELSE        ELSE
95         _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
96  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
97  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.
98         CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )          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  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  C The 0 is the "iteration" argument. The ' ' is an empty suffix
101  C      CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
102  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
103  C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
104  C    &                    'RS', 1, H, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
105         _END_MASTER(myThid)          _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        ENDIF
187    
188        _EXCH_XY_R4(    H, myThid )  C----- fill in the overlap :
189          _EXCH_XY_R4(Ro_surf, myThid )
190    
191  C  c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192  C     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)  
193  C  C------
194    C   3) Close the Domain (special configuration).
195    C------
196        IF (groundAtK1) THEN        IF (groundAtK1) THEN
197         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
198          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
199           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
200            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
201             R_low(I,J,bi,bj) = rF(Nr+1)             iG = myXGlobalLo-1+(bi-1)*sNx+I
202             Ro_surf(I,J,bi,bj) = H(I,J,bi,bj)             jG = myYGlobalLo-1+(bj-1)*sNy+J
203    C          Test for eastern edge
204    c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
205    C          Test for northern edge
206    c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
207  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
208             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
209       &         abs(yC(I,J,bi,bj)).GE.-phiMin) Ro_surf(I,J,bi,bj) = 0.       &         abs(yC(I,J,bi,bj)).GE.-phiMin)
210         &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
211             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
212       &     THEN       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
             Ro_surf(I,J,bi,bj) = 0.  
            ENDIF  
213            ENDDO            ENDDO
214           ENDDO           ENDDO
215          ENDDO          ENDDO
# Line 121  c- Domain : Symetric % Eq. & closed at N Line 219  c- Domain : Symetric % Eq. & closed at N
219          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
220           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
221            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
222             R_low(I,J,bi,bj) = H(I,J,bi,bj)             iG = myXGlobalLo-1+(bi-1)*sNx+I
223             Ro_surf(I,J,bi,bj) = Ro_SeaLevel             jG = myYGlobalLo-1+(bj-1)*sNy+J
224    C          Test for eastern edge
225    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
226    C          Test for northern edge
227    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
228  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
229             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
230       &         abs(yC(I,J,bi,bj)).GE.-phiMin) R_low(I,J,bi,bj) = 0.       &         abs(yC(I,J,bi,bj)).GE.-phiMin)
231         &       R_low(I,J,bi,bj) = Ro_SeaLevel
232             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
233       &     THEN       &       R_low(I,J,bi,bj) = Ro_SeaLevel
             R_low(I,J,bi,bj) = 0.  
            ENDIF  
234            ENDDO            ENDDO
235           ENDDO           ENDDO
236          ENDDO          ENDDO
237         ENDDO         ENDDO
238        ENDIF        ENDIF
239    
240    c     _BEGIN_MASTER( myThid )
241    c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
242    c     _END_MASTER(myThid)
243    
244        RETURN        RETURN
245        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22