/[MITgcm]/MITgcm/verification/aim.5l_Equatorial_Channel/code/ini_depths.F
ViewVC logotype

Annotation of /MITgcm/verification/aim.5l_Equatorial_Channel/code/ini_depths.F

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


Revision 1.5 - (hide annotations) (download)
Fri Aug 24 11:57:43 2001 UTC (22 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint40pre9, checkpoint42, checkpoint40, checkpoint41
Changes since 1.4: +158 -53 lines
Updated after changes (05-29-01) in AIM package (multi-threaded)
CV: ----------------------------------------------------------------------

1 jmc 1.5 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_depths.F,v 1.25 2001/07/09 16:46:29 jmc Exp $
2 adcroft 1.4 C $Name: $
3 adcroft 1.2
4     #include "CPP_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE INI_DEPTHS( myThid )
8     C /==========================================================\
9     C | SUBROUTINE INI_DEPTHS |
10 jmc 1.5 C | o define R_position of Lower and Surface Boundaries |
11 adcroft 1.2 C |==========================================================|
12 jmc 1.5 C |atmosphere orography: |
13     C | define either in term of P_topo or converted from Z_topo |
14     C |ocean bathymetry: |
15 adcroft 1.2 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 |
17     C | grid cells. Depths do not have to coincide with the |
18     C | model levels. The model lopping algorithm makes it |
19     C | possible to represent arbitrary depths. |
20     C | The mode depths map also influences the models topology |
21     C | By default the model domain wraps around in X and Y. |
22     C | This default doubly periodic topology is "supressed" |
23     C | if a depth map is defined which closes off all wrap |
24     C | around flow. |
25     C \==========================================================/
26     IMPLICIT NONE
27    
28     C === Global variables ===
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33    
34     C == Routine arguments ==
35     C myThid - Number of this instance of INI_DEPTHS
36     INTEGER myThid
37     CEndOfInterface
38    
39 jmc 1.5 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 adcroft 1.2 C == Local variables ==
46     C iG, jG - Global coordinate index
47     C bi,bj - Loop counters
48     C I,J,K
49     C oldPrec - Temporary used in controlling binary input dataset precision
50 jmc 1.5 C msgBuf - Informational/error meesage buffer
51 adcroft 1.2 INTEGER iG, jG
52     INTEGER bi, bj
53 adcroft 1.4 INTEGER I, J
54 jmc 1.5 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 adcroft 1.2
65 jmc 1.5 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 adcroft 1.2 DO bj = myByLo(myThid), myByHi(myThid)
86     DO bi = myBxLo(myThid), myBxHi(myThid)
87     DO j=1,sNy
88     DO i=1,sNx
89 jmc 1.5 R_low(i,j,bi,bj) = rF(Nr+1)
90 adcroft 1.2 ENDDO
91     ENDDO
92     ENDDO
93     ENDDO
94     ELSE
95 jmc 1.5 _BEGIN_MASTER( myThid )
96 adcroft 1.2 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 jmc 1.5 CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
99 adcroft 1.2 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 jmc 1.5 c CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
102 adcroft 1.2 C Read the bathymetry using the low-level I/O package
103 jmc 1.5 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 adcroft 1.2 ENDIF
187    
188 jmc 1.5 C----- fill in the overlap :
189     _EXCH_XY_R4(Ro_surf, myThid )
190 adcroft 1.2
191 jmc 1.5 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
192    
193     C------
194     C 3) Close the Domain (special configuration).
195     C------
196 adcroft 1.4 IF (groundAtK1) THEN
197 adcroft 1.2 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 jmc 1.5 iG = myXGlobalLo-1+(bi-1)*sNx+I
202     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 adcroft 1.4 c- Domain : Symetric % Eq. & closed at N & S boundaries:
208     IF (usingSphericalPolarGrid .AND.
209 jmc 1.5 & abs(yC(I,J,bi,bj)).GE.-phiMin)
210     & Ro_surf(I,J,bi,bj) = rF(Nr+1)
211 adcroft 1.4 IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
212 jmc 1.5 & Ro_surf(I,J,bi,bj) = rF(Nr+1)
213 adcroft 1.4 ENDDO
214     ENDDO
215     ENDDO
216     ENDDO
217     ELSE
218     DO bj = myByLo(myThid), myByHi(myThid)
219     DO bi = myBxLo(myThid), myBxHi(myThid)
220     DO j=1-Oly,sNy+Oly
221     DO i=1-Olx,sNx+Olx
222 jmc 1.5 iG = myXGlobalLo-1+(bi-1)*sNx+I
223     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 adcroft 1.4 c- Domain : Symetric % Eq. & closed at N & S boundaries:
229     IF (usingSphericalPolarGrid .AND.
230 jmc 1.5 & abs(yC(I,J,bi,bj)).GE.-phiMin)
231     & R_low(I,J,bi,bj) = Ro_SeaLevel
232 adcroft 1.4 IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
233 jmc 1.5 & R_low(I,J,bi,bj) = Ro_SeaLevel
234 adcroft 1.2 ENDDO
235     ENDDO
236     ENDDO
237     ENDDO
238     ENDIF
239 jmc 1.5
240     c _BEGIN_MASTER( myThid )
241     c CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
242     c _END_MASTER(myThid)
243 adcroft 1.4
244 adcroft 1.2 RETURN
245     END

  ViewVC Help
Powered by ViewVC 1.1.22