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

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

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

revision 1.1 by adcroft, Fri Sep 28 02:28:10 2001 UTC revision 1.2 by jmc, Tue Dec 10 03:07:28 2002 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.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 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        IMPLICIT NONE  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  #include "INI_DEPTHS.h"  #include "SURFACE.h"
39    
40    C     !INPUT/OUTPUT PARAMETERS:
41  C     == Routine arguments ==  C     == Routine arguments ==
42  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid -  Number of this instance of INI_DEPTHS
43        INTEGER myThid        INTEGER myThid
44  CEndOfInterface  CEndOfInterface
45    
46    C     !LOCAL VARIABLES:
47  C     == Local variables ==  C     == Local variables ==
48  C     iG, jG - Global coordinate index  C     iG, jG - Global coordinate index
49  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
50  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--------------------  
51  C     oldPrec - Temporary used in controlling binary input dataset precision  C     oldPrec - Temporary used in controlling binary input dataset precision
52    C     msgBuf    - Informational/error meesage buffer
53        INTEGER iG, jG        INTEGER iG, jG
54        INTEGER bi, bj        INTEGER bi, bj
55        INTEGER  I, J        INTEGER  I, J
56        _RL Hdefault        CHARACTER*(MAX_LEN_MBUF) msgBuf
57    CEOP
58    
59        _BARRIER        IF (groundAtK1 .AND. bathyFile .NE. ' '
60        IF ( bathyFile .EQ. ' '  ) THEN       &               .AND. topoFile  .NE. ' ' ) THEN
61  C      Set up a flat bottom box with doubly periodic topology.         WRITE(msgBuf,'(A,A)')
62  C      H is the basic variable from which other terms are derived. It       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
63  C      is the term that would be set from an external file for a       &  ' select the right one !'
64  C      realistic problem.         CALL PRINT_ERROR( msgBuf , myThid)
65         IF (groundAtK1) THEN         STOP 'ABNORMAL END: S/R INI_DEPTHS'
66           Hdefault = Ro_SeaLevel        ENDIF
67         ELSE  
68           Hdefault = rF(Nr+1)  C------
69         ENDIF  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)         DO bj = myByLo(myThid), myByHi(myThid)
90          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
91           DO j=1,sNy           DO j=1,sNy
92            DO i=1,sNx            DO i=1,sNx
93             iG = myXGlobalLo-1+(bi-1)*sNx+I             R_low(i,j,bi,bj) = rF(Nr+1)
94             jG = myYGlobalLo-1+(bj-1)*sNy+J  C-- Specific modif for this experiment (advect_xz):
95  C          Default depth of full domain             R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
            H(i,j,bi,bj) = Hdefault  
96       &     *(1.-0.5*XC(I,J,bi,bj)/(float(Nx)*DelX(1)))       &     *(1.-0.5*XC(I,J,bi,bj)/(float(Nx)*DelX(1)))
97  C          Test for eastern edge  C-- end of modified part
 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.  
98            ENDDO            ENDDO
99           ENDDO           ENDDO
100          ENDDO          ENDDO
101         ENDDO         ENDDO
102        ELSE        ELSE
103         _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
104  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
105  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.
106         CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
107  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
108  C The 0 is the "iteration" argument. The ' ' is an empty suffix  C The 0 is the "iteration" argument. The ' ' is an empty suffix
109  C      CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
110  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
111  C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
112  C    &                    'RS', 1, H, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
113         _END_MASTER(myThid)          _END_MASTER(myThid)
114    
115        ENDIF        ENDIF
116    C- end setup R_low in the interior
117    
118        _EXCH_XY_R4(    H, myThid )  C- fill in the overlap :
119          _EXCH_XY_R4(R_low, myThid )
120    
121  C  c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
122  C     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)  c     _BEGIN_MASTER( myThid )
123  C  c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
124    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          IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN
133    C------ read directly Po_surf from bathyFile (only for backward compatibility)
134    
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    C   Convert Surface Geopotential to (reference) Surface Pressure
165    C   according to Tref profile, using same discretisation as in calc_phi_hyd
166    C----
167    c         _BEGIN_MASTER( myThid )
168    c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
169    c         _END_MASTER(myThid)
170    
171              CALL INI_P_GROUND( selectFindRoSurf, topoZ,
172         O                       Ro_surf,
173         I                       myThid )
174    
175              _BARRIER
176              _BEGIN_MASTER( myThid )
177              CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
178              _END_MASTER(myThid)
179    
180            ELSE
181    C----
182    C   Direct Transfer to Ro_surf :
183              DO bj = myByLo(myThid), myByHi(myThid)
184               DO bi = myBxLo(myThid), myBxHi(myThid)
185                DO j=1,sNy
186                 DO i=1,sNx
187                  Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
188                 ENDDO
189                ENDDO
190               ENDDO
191              ENDDO
192    
193            ENDIF
194    
195    C------ end case "read topoFile"
196          ENDIF
197    
198    C----- fill in the overlap :
199          _EXCH_XY_R4(Ro_surf, myThid )
200    
201    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
202    
203    C------
204    C   3) Close the Domain (special configuration).
205    C------
206        IF (groundAtK1) THEN        IF (groundAtK1) THEN
207         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
208          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
209           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
210            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
211             R_low(I,J,bi,bj) = rF(Nr+1)             iG = myXGlobalLo-1+(bi-1)*sNx+I
212             Ro_surf(I,J,bi,bj) = H(I,J,bi,bj)             jG = myYGlobalLo-1+(bj-1)*sNy+J
213    C          Test for eastern edge
214    c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
215    C          Test for northern edge
216    c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
217             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
218       &     THEN       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
             Ro_surf(I,J,bi,bj) = 0.  
            ENDIF  
219            ENDDO            ENDDO
220           ENDDO           ENDDO
221          ENDDO          ENDDO
# Line 119  C Line 225  C
225          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
226           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
227            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
228             R_low(I,J,bi,bj) = H(I,J,bi,bj)             iG = myXGlobalLo-1+(bi-1)*sNx+I
229             Ro_surf(I,J,bi,bj) = Ro_SeaLevel             jG = myYGlobalLo-1+(bj-1)*sNy+J
230    C          Test for eastern edge
231    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
232    C          Test for northern edge
233    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
234             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
235       &     THEN       &       R_low(I,J,bi,bj) = Ro_SeaLevel
             R_low(I,J,bi,bj) = 0.  
            ENDIF  
236            ENDDO            ENDDO
237           ENDDO           ENDDO
238          ENDDO          ENDDO
239         ENDDO         ENDDO
240        ENDIF        ENDIF
241    
242    c     _BEGIN_MASTER( myThid )
243    c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
244    c     _END_MASTER(myThid)
245    
246        RETURN        RETURN
247        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22