/[MITgcm]/MITgcm_contrib/natl_12/code/ini_depths.F
ViewVC logotype

Annotation of /MITgcm_contrib/natl_12/code/ini_depths.F

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


Revision 1.2 - (hide annotations) (download)
Thu Aug 7 13:03:31 2003 UTC (21 years, 11 months ago) by cnh
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +5 -2 lines
Mods to allow laplacian as well as biharmonic viscosity with varying resolution

1 cnh 1.2 C $Header: /u/gcmpack/MITgcm_contrib/natl_12/code/ini_depths.F,v 1.1 2003/08/05 21:22:43 cnh Exp $
2 cnh 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: INI_DEPTHS
8     C !INTERFACE:
9     SUBROUTINE INI_DEPTHS( myThid )
10     C !DESCRIPTION: \bv
11     C *==========================================================*
12     C | SUBROUTINE INI_DEPTHS
13     C | o define R_position of Lower and Surface Boundaries
14     C *==========================================================*
15     C |atmosphere orography:
16     C | define either in term of P_topo or converted from Z_topo
17     C |ocean bathymetry:
18     C | The depths of the bottom of the model is specified in
19     C | terms of an XY map with one depth for each column of
20     C | grid cells. Depths do not have to coincide with the
21     C | model levels. The model lopping algorithm makes it
22     C | possible to represent arbitrary depths.
23     C | The mode depths map also influences the models topology
24     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 ===
34     #include "SIZE.h"
35     #include "EEPARAMS.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38     #include "SURFACE.h"
39    
40     C !INPUT/OUTPUT PARAMETERS:
41     C == Routine arguments ==
42     C myThid - Number of this instance of INI_DEPTHS
43     INTEGER myThid
44     CEndOfInterface
45    
46     C !LOCAL VARIABLES:
47     C == Local variables ==
48     C iG, jG - Global coordinate index
49     C bi,bj - Tile indices
50     C I,J,K - Loop counters
51     C oldPrec - Temporary used in controlling binary input dataset precision
52     C msgBuf - Informational/error meesage buffer
53     INTEGER iG, jG
54     INTEGER bi, bj
55     INTEGER I, J
56     CHARACTER*(MAX_LEN_MBUF) msgBuf
57     CEOP
58    
59     IF (groundAtK1 .AND. bathyFile .NE. ' '
60     & .AND. topoFile .NE. ' ' ) THEN
61     WRITE(msgBuf,'(A,A)')
62     & 'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
63     & ' select the right one !'
64     CALL PRINT_ERROR( msgBuf , myThid)
65     STOP 'ABNORMAL END: S/R INI_DEPTHS'
66     ENDIF
67    
68     C------
69     C 0) Initialize R_low and Ro_surf (define an empty domain)
70     C------
71     DO bj = myByLo(myThid), myByHi(myThid)
72     DO bi = myBxLo(myThid), myBxHi(myThid)
73     DO j=1-Oly,sNy+Oly
74     DO i=1-Olx,sNx+Olx
75     R_low(i,j,bi,bj) = 0.
76     Ro_surf(i,j,bi,bj) = 0.
77     topoZ(i,j,bi,bj) = 0.
78     ENDDO
79     ENDDO
80     ENDDO
81     ENDDO
82    
83     C------
84     C 1) Set R_low = the Lower (in r sense) boundary of the fluid column :
85     C------
86     IF (groundAtK1 .OR. bathyFile .EQ. ' ') THEN
87     C- e.g., atmosphere : R_low = Top of atmosphere
88     C- ocean : R_low = Bottom
89     DO bj = myByLo(myThid), myByHi(myThid)
90     DO bi = myBxLo(myThid), myBxHi(myThid)
91     DO j=1,sNy
92     DO i=1,sNx
93     R_low(i,j,bi,bj) = rF(Nr+1)
94     ENDDO
95     ENDDO
96     ENDDO
97     ENDDO
98     ELSE
99     _BEGIN_MASTER( myThid )
100     C Read the bathymetry using the mid-level I/O pacakage read_write_rec
101     C The 0 is the "iteration" argument. The 1 is the record number.
102     CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
103     C Read the bathymetry using the mid-level I/O pacakage read_write_fld
104     C The 0 is the "iteration" argument. The ' ' is an empty suffix
105     c CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
106     C Read the bathymetry using the low-level I/O package
107     c CALL MDSREADFIELD( bathyFile, readBinaryPrec,
108     c & 'RS', 1, R_low, 1, myThid )
109     _END_MASTER(myThid)
110    
111     ENDIF
112     C- end setup R_low in the interior
113    
114     C- fill in the overlap :
115     _EXCH_XY_R4(R_low, myThid )
116    
117 cnh 1.2 CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
118 cnh 1.1 c _BEGIN_MASTER( myThid )
119     c CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
120     c _END_MASTER(myThid)
121    
122     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
123    
124     C------
125     C 2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
126     C------
127    
128     IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN
129     C------ read directly Po_surf from bathyFile (only for backward compatibility)
130    
131     _BEGIN_MASTER( myThid )
132     CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
133     _END_MASTER(myThid)
134     _BARRIER
135    
136     ELSEIF ( topoFile.EQ.' ' ) THEN
137     C------ set default value:
138    
139     DO bj = myByLo(myThid), myByHi(myThid)
140     DO bi = myBxLo(myThid), myBxHi(myThid)
141     DO j=1,sNy
142     DO i=1,sNx
143     Ro_surf(i,j,bi,bj) = Ro_SeaLevel
144     ENDDO
145     ENDDO
146     ENDDO
147     ENDDO
148    
149     ELSE
150     C------ read from file:
151    
152     C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
153     _BEGIN_MASTER( myThid )
154     CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
155     _END_MASTER(myThid)
156     _BARRIER
157    
158     IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
159     C----
160     C Convert Surface Geopotential to (reference) Surface Pressure
161     C according to Tref profile, using same discretisation as in calc_phi_hyd
162     C----
163     c _BEGIN_MASTER( myThid )
164     c CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
165     c _END_MASTER(myThid)
166    
167     CALL INI_P_GROUND( 2, topoZ,
168     O Ro_surf,
169     I myThid )
170    
171     _BARRIER
172     _BEGIN_MASTER( myThid )
173     CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
174     _END_MASTER(myThid)
175    
176     ELSE
177     C----
178     C Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
179     C below an ice-shelf - NOTE - actually not yet implemented )
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) = topoZ(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     Ccnh Begin custom code for natl_12
199     DO bj = myByLo(myThid), myByHi(myThid)
200     DO bi = myBxLo(myThid), myBxHi(myThid)
201     DO j=1-Oly,sNy+Oly
202     DO i=1-Olx,sNx+Olx
203     iG = myXGlobalLo-1+(bi-1)*sNx+I
204     jG = myYGlobalLo-1+(bj-1)*sNy+J
205     C Test for northern or southern overlap
206     IF ( jG .GT. Ny ) Ro_surf(i,j,bi,bj) = 0.
207     IF ( jG .LT. 1 ) Ro_surf(i,j,bi,bj) = 0.
208     IF ( jG .GT. Ny ) R_Low( i,j,bi,bj) = 0.
209     IF ( jG .LT. 1 ) R_Low( i,j,bi,bj) = 0.
210     IF ( jG .EQ. Ny ) THEN
211     IF ( iG .GE. 347 .AND. iG .LE. 357 ) THEN
212     R_Low( i,j,bi,bj) = 0
213     ENDIF
214     ENDIF
215     IF ( jG .GE. 286 .AND. jG .LE. 311 .AND. iG .LE. 390 ) THEN
216     R_Low( i,j,bi,bj) = 0
217     ENDIF
218     IF ( jG .GT. 311 .AND. iG .LE. 350 ) THEN
219     R_Low( i,j,bi,bj) = 0
220     ENDIF
221     IF ( jG .GT. 308 .AND. iG .LE. 397 ) THEN
222     R_Low( i,j,bi,bj) = 0
223     ENDIF
224     ENDDO
225     ENDDO
226     ENDDO
227     ENDDO
228     Ccnh End custom code for natl_12
229 cnh 1.2 CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths after custom)',1,myThid)
230 cnh 1.1
231     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
232    
233     C------
234     C 3) Close the Domain (special configuration).
235     C------
236     IF (groundAtK1) THEN
237     DO bj = myByLo(myThid), myByHi(myThid)
238     DO bi = myBxLo(myThid), myBxHi(myThid)
239     DO j=1-Oly,sNy+Oly
240     DO i=1-Olx,sNx+Olx
241     iG = myXGlobalLo-1+(bi-1)*sNx+I
242     jG = myYGlobalLo-1+(bj-1)*sNy+J
243     C Test for eastern edge
244     c IF ( iG .EQ. Nx ) Ro_surf(i,j,bi,bj) = 0.
245     C Test for northern edge
246     c IF ( jG .EQ. Ny ) Ro_surf(i,j,bi,bj) = 0.
247     IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
248     & Ro_surf(I,J,bi,bj) = rF(Nr+1)
249     ENDDO
250     ENDDO
251     ENDDO
252     ENDDO
253     ELSE
254     DO bj = myByLo(myThid), myByHi(myThid)
255     DO bi = myBxLo(myThid), myBxHi(myThid)
256     DO j=1-Oly,sNy+Oly
257     DO i=1-Olx,sNx+Olx
258     iG = myXGlobalLo-1+(bi-1)*sNx+I
259     jG = myYGlobalLo-1+(bj-1)*sNy+J
260     C Test for eastern edge
261     c IF ( iG .EQ. Nx ) R_low(i,j,bi,bj) = 0.
262     C Test for northern edge
263     c IF ( jG .EQ. Ny ) R_low(i,j,bi,bj) = 0.
264     IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
265     & R_low(I,J,bi,bj) = Ro_SeaLevel
266     ENDDO
267     ENDDO
268     ENDDO
269     ENDDO
270     ENDIF
271 cnh 1.2
272     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths routine end)',1,myThid)
273 cnh 1.1
274     c _BEGIN_MASTER( myThid )
275     c CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
276     c _END_MASTER(myThid)
277    
278     RETURN
279     END

  ViewVC Help
Powered by ViewVC 1.1.22