/[MITgcm]/MITgcm/verification/advect_xz/code/ini_depths.F
ViewVC logotype

Annotation of /MITgcm/verification/advect_xz/code/ini_depths.F

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


Revision 1.3 - (hide annotations) (download)
Wed Jun 22 00:18:48 2005 UTC (18 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint57m_post, checkpoint57v_post, checkpoint57s_post, checkpoint57j_post, checkpoint57y_post, checkpoint57x_post, checkpoint57y_pre, checkpoint57o_post, checkpoint57r_post, checkpoint57k_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, checkpoint57z_post, checkpoint57l_post
Changes since 1.2: +41 -39 lines
"usingPCoords" replaces "groundAtK1" (<- removed)

1 jmc 1.3 C $Header: /u/gcmpack/MITgcm/verification/advect_xz/code/ini_depths.F,v 1.2 2002/12/10 03:07:28 jmc Exp $
2 adcroft 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6 jmc 1.2 CBOP
7     C !ROUTINE: INI_DEPTHS
8     C !INTERFACE:
9 adcroft 1.1 SUBROUTINE INI_DEPTHS( myThid )
10 jmc 1.2 C !DESCRIPTION: \bv
11     C *==========================================================*
12 jmc 1.3 C | SUBROUTINE INI_DEPTHS
13     C | o define R_position of Lower and Surface Boundaries
14 jmc 1.2 C *==========================================================*
15 jmc 1.3 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 jmc 1.2 C *==========================================================*
29     C \ev
30    
31     C !USES:
32 adcroft 1.1 IMPLICIT NONE
33     C === Global variables ===
34     #include "SIZE.h"
35     #include "EEPARAMS.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38 jmc 1.2 #include "SURFACE.h"
39 adcroft 1.1
40 jmc 1.2 C !INPUT/OUTPUT PARAMETERS:
41 adcroft 1.1 C == Routine arguments ==
42     C myThid - Number of this instance of INI_DEPTHS
43     INTEGER myThid
44     CEndOfInterface
45    
46 jmc 1.2 C !LOCAL VARIABLES:
47 adcroft 1.1 C == Local variables ==
48     C iG, jG - Global coordinate index
49 jmc 1.3 C bi,bj - Tile indices
50     C I,J,K - Loop counters
51 adcroft 1.1 C oldPrec - Temporary used in controlling binary input dataset precision
52 jmc 1.3 C msgBuf - Informational/error meesage buffer
53 adcroft 1.1 INTEGER iG, jG
54     INTEGER bi, bj
55     INTEGER I, J
56 jmc 1.2 CHARACTER*(MAX_LEN_MBUF) msgBuf
57     CEOP
58 adcroft 1.1
59 jmc 1.3 IF (usingPCoords .AND. bathyFile .NE. ' '
60     & .AND. topoFile .NE. ' ' ) THEN
61 jmc 1.2 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 jmc 1.3 IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
87 jmc 1.2 C- e.g., atmosphere : R_low = Top of atmosphere
88     C- ocean : R_low = Bottom
89 adcroft 1.1 DO bj = myByLo(myThid), myByHi(myThid)
90     DO bi = myBxLo(myThid), myBxHi(myThid)
91     DO j=1,sNy
92     DO i=1,sNx
93 jmc 1.2 R_low(i,j,bi,bj) = rF(Nr+1)
94     C-- Specific modif for this experiment (advect_xz):
95 jmc 1.3 R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
96 adcroft 1.1 & *(1.-0.5*XC(I,J,bi,bj)/(float(Nx)*DelX(1)))
97 jmc 1.2 C-- end of modified part
98 adcroft 1.1 ENDDO
99     ENDDO
100     ENDDO
101     ENDDO
102     ELSE
103 jmc 1.2 _BEGIN_MASTER( myThid )
104 adcroft 1.1 C Read the bathymetry using the mid-level I/O pacakage read_write_rec
105     C The 0 is the "iteration" argument. The 1 is the record number.
106 jmc 1.2 CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
107 adcroft 1.1 C Read the bathymetry using the mid-level I/O pacakage read_write_fld
108     C The 0 is the "iteration" argument. The ' ' is an empty suffix
109 jmc 1.2 c CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
110 adcroft 1.1 C Read the bathymetry using the low-level I/O package
111 jmc 1.2 c CALL MDSREADFIELD( bathyFile, readBinaryPrec,
112     c & 'RS', 1, R_low, 1, myThid )
113     _END_MASTER(myThid)
114    
115 adcroft 1.1 ENDIF
116 jmc 1.2 C- end setup R_low in the interior
117 adcroft 1.1
118 jmc 1.2 C- fill in the overlap :
119     _EXCH_XY_R4(R_low, myThid )
120 adcroft 1.1
121 jmc 1.2 c CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
122     c _BEGIN_MASTER( myThid )
123 jmc 1.3 c CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
124 jmc 1.2 c _END_MASTER(myThid)
125    
126     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
127    
128     C------
129     C 2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
130     C------
131    
132 jmc 1.3 IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
133     C------ read directly Po_surf from bathyFile (only for backward compatibility)
134 jmc 1.2
135     _BEGIN_MASTER( myThid )
136     CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
137     _END_MASTER(myThid)
138     _BARRIER
139    
140     ELSEIF ( topoFile.EQ.' ' ) THEN
141     C------ set default value:
142    
143     DO bj = myByLo(myThid), myByHi(myThid)
144     DO bi = myBxLo(myThid), myBxHi(myThid)
145     DO j=1,sNy
146     DO i=1,sNx
147     Ro_surf(i,j,bi,bj) = Ro_SeaLevel
148     ENDDO
149     ENDDO
150     ENDDO
151     ENDDO
152    
153     ELSE
154     C------ read from file:
155    
156     C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
157     _BEGIN_MASTER( myThid )
158     CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
159     _END_MASTER(myThid)
160     _BARRIER
161    
162     IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
163     C----
164 jmc 1.3 C Convert Surface Geopotential to (reference) Surface Pressure
165 jmc 1.2 C according to Tref profile, using same discretisation as in calc_phi_hyd
166     C----
167     c _BEGIN_MASTER( myThid )
168 jmc 1.3 c CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
169 jmc 1.2 c _END_MASTER(myThid)
170    
171 jmc 1.3 CALL INI_P_GROUND( 2, topoZ,
172     O Ro_surf,
173 jmc 1.2 I myThid )
174    
175     _BARRIER
176 jmc 1.3 C This I/O is now done in write_grid.F
177     c _BEGIN_MASTER( myThid )
178     c CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
179     c _END_MASTER(myThid)
180 jmc 1.2
181     ELSE
182     C----
183 jmc 1.3 C Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
184     C below an ice-shelf - NOTE - actually not yet implemented )
185 jmc 1.2 DO bj = myByLo(myThid), myByHi(myThid)
186     DO bi = myBxLo(myThid), myBxHi(myThid)
187     DO j=1,sNy
188     DO i=1,sNx
189     Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
190     ENDDO
191     ENDDO
192     ENDDO
193     ENDDO
194    
195     ENDIF
196    
197     C------ end case "read topoFile"
198     ENDIF
199    
200 jmc 1.3 C----- fill in the overlap :
201 jmc 1.2 _EXCH_XY_R4(Ro_surf, myThid )
202    
203     c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
204    
205     C------
206     C 3) Close the Domain (special configuration).
207     C------
208 jmc 1.3 IF (usingPCoords) THEN
209 adcroft 1.1 DO bj = myByLo(myThid), myByHi(myThid)
210     DO bi = myBxLo(myThid), myBxHi(myThid)
211     DO j=1-Oly,sNy+Oly
212     DO i=1-Olx,sNx+Olx
213 jmc 1.2 iG = myXGlobalLo-1+(bi-1)*sNx+I
214     jG = myYGlobalLo-1+(bj-1)*sNy+J
215     C Test for eastern edge
216     c IF ( iG .EQ. Nx ) Ro_surf(i,j,bi,bj) = 0.
217     C Test for northern edge
218     c IF ( jG .EQ. Ny ) Ro_surf(i,j,bi,bj) = 0.
219 jmc 1.3 IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
220 jmc 1.2 & Ro_surf(I,J,bi,bj) = rF(Nr+1)
221 adcroft 1.1 ENDDO
222     ENDDO
223     ENDDO
224     ENDDO
225     ELSE
226     DO bj = myByLo(myThid), myByHi(myThid)
227     DO bi = myBxLo(myThid), myBxHi(myThid)
228     DO j=1-Oly,sNy+Oly
229     DO i=1-Olx,sNx+Olx
230 jmc 1.2 iG = myXGlobalLo-1+(bi-1)*sNx+I
231     jG = myYGlobalLo-1+(bj-1)*sNy+J
232     C Test for eastern edge
233     c IF ( iG .EQ. Nx ) R_low(i,j,bi,bj) = 0.
234     C Test for northern edge
235     c IF ( jG .EQ. Ny ) R_low(i,j,bi,bj) = 0.
236 jmc 1.3 IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
237 jmc 1.2 & R_low(I,J,bi,bj) = Ro_SeaLevel
238 adcroft 1.1 ENDDO
239     ENDDO
240     ENDDO
241     ENDDO
242     ENDIF
243 jmc 1.2
244     c _BEGIN_MASTER( myThid )
245 jmc 1.3 c CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
246 jmc 1.2 c _END_MASTER(myThid)
247 adcroft 1.1
248     RETURN
249     END

  ViewVC Help
Powered by ViewVC 1.1.22