/[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.9 by jmc, Mon Mar 17 16:42:07 2003 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  - Tile indices
50  C     I,J,K  C     I,J,K  - Loop counters
 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)
            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.  
94            ENDDO            ENDDO
95           ENDDO           ENDDO
96          ENDDO          ENDDO
97         ENDDO         ENDDO
98        ELSE        ELSE
99         _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
100  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
101  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.
102         CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )          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  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  C The 0 is the "iteration" argument. The ' ' is an empty suffix
105  C      CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
106  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
107  C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
108  C    &                    'RS', 1, H, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
109         _END_MASTER(myThid)          _END_MASTER(myThid)
110    
111        ENDIF        ENDIF
112    C- end setup R_low in the interior
113    
114        _EXCH_XY_R4(    H, myThid )  C- fill in the overlap :
115          _EXCH_XY_R4(R_low, myThid )
116    
117  C  c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
118  C     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)  c     _BEGIN_MASTER( myThid )
119  C  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    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
199    
200    C------
201    C   3) Close the Domain (special configuration).
202    C------
203        IF (groundAtK1) THEN        IF (groundAtK1) THEN
204         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
205          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
206           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
207            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
208             R_low(I,J,bi,bj) = rF(Nr+1)             iG = myXGlobalLo-1+(bi-1)*sNx+I
209             Ro_surf(I,J,bi,bj) = H(I,J,bi,bj)             jG = myYGlobalLo-1+(bj-1)*sNy+J
210    C          Test for eastern edge
211    c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
212    C          Test for northern edge
213    c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
214  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
215             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
216       &         abs(yC(I,J,bi,bj)).GE.-phiMin) Ro_surf(I,J,bi,bj) = 0.       &         abs(yC(I,J,bi,bj)).GE.-phiMin)
217         &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
218             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
219       &     THEN       &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
             Ro_surf(I,J,bi,bj) = 0.  
            ENDIF  
220            ENDDO            ENDDO
221           ENDDO           ENDDO
222          ENDDO          ENDDO
# Line 121  c- Domain : Symetric % Eq. & closed at N Line 226  c- Domain : Symetric % Eq. & closed at N
226          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
227           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
228            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
229             R_low(I,J,bi,bj) = H(I,J,bi,bj)             iG = myXGlobalLo-1+(bi-1)*sNx+I
230             Ro_surf(I,J,bi,bj) = Ro_SeaLevel             jG = myYGlobalLo-1+(bj-1)*sNy+J
231    C          Test for eastern edge
232    c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
233    C          Test for northern edge
234    c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
235  c- Domain : Symetric % Eq. & closed at N & S boundaries:  c- Domain : Symetric % Eq. & closed at N & S boundaries:
236             IF (usingSphericalPolarGrid .AND.             IF (usingSphericalPolarGrid .AND.
237       &         abs(yC(I,J,bi,bj)).GE.-phiMin) R_low(I,J,bi,bj) = 0.       &         abs(yC(I,J,bi,bj)).GE.-phiMin)
238         &       R_low(I,J,bi,bj) = Ro_SeaLevel
239             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )             IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
240       &     THEN       &       R_low(I,J,bi,bj) = Ro_SeaLevel
             R_low(I,J,bi,bj) = 0.  
            ENDIF  
241            ENDDO            ENDDO
242           ENDDO           ENDDO
243          ENDDO          ENDDO
244         ENDDO         ENDDO
245        ENDIF        ENDIF
246    
247    c     _BEGIN_MASTER( myThid )
248    c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
249    c     _END_MASTER(myThid)
250    
251        RETURN        RETURN
252        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22