/[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.6 by jmc, Thu Sep 27 23:42:16 2001 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"
 #include "INI_DEPTHS.h"  
38    
39    C     !INPUT/OUTPUT PARAMETERS:
40  C     == Routine arguments ==  C     == Routine arguments ==
41  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid -  Number of this instance of INI_DEPTHS
42        INTEGER myThid        INTEGER myThid
43  CEndOfInterface  CEndOfInterface
44    
45    C     !LOCAL VARIABLES:
46    C     == Local variables in common ==
47    C     Hloc  - Temporary array used to read surface topography
48    C             has to be in common for multi threading    
49          COMMON / LOCAL_INI_DEPTHS / Hloc
50          _RS Hloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51  C     == Local variables ==  C     == Local variables ==
52  C     iG, jG - Global coordinate index  C     iG, jG - Global coordinate index
53  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
54  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--------------------  
55  C     oldPrec - Temporary used in controlling binary input dataset precision  C     oldPrec - Temporary used in controlling binary input dataset precision
56    C     msgBuf    - Informational/error meesage buffer
57        INTEGER iG, jG        INTEGER iG, jG
58        INTEGER bi, bj        INTEGER bi, bj
59        INTEGER  I, J        INTEGER  I, J
60        _RL Hdefault        CHARACTER*(MAX_LEN_MBUF) msgBuf
61    CEOP
62    
63        _BARRIER        IF (groundAtK1 .AND. bathyFile .NE. ' '
64        IF ( bathyFile .EQ. ' '  ) THEN       &               .AND. topoFile  .NE. ' ' ) THEN
65  C      Set up a flat bottom box with doubly periodic topology.         WRITE(msgBuf,'(A,A)')
66  C      H is the basic variable from which other terms are derived. It       &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
67  C      is the term that would be set from an external file for a       &  ' select the right one !'
68  C      realistic problem.         CALL PRINT_ERROR( msgBuf , myThid)
69         IF (groundAtK1) THEN         STOP 'ABNORMAL END: S/R INI_DEPTHS'
70           Hdefault = Ro_SeaLevel        ENDIF
71         ELSE  
72           Hdefault = rF(Nr+1)  C------
73         ENDIF  C   0) Initialize R_low and Ro_surf (define an empty domain)
74    C------
75          DO bj = myByLo(myThid), myByHi(myThid)
76           DO bi = myBxLo(myThid), myBxHi(myThid)
77            DO j=1-Oly,sNy+Oly
78             DO i=1-Olx,sNx+Olx
79              R_low(i,j,bi,bj) = 0.
80              Ro_surf(i,j,bi,bj) = 0.
81             ENDDO
82            ENDDO
83           ENDDO
84          ENDDO
85    
86    C------
87    C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
88    C------
89          IF (groundAtK1 .OR. bathyFile .EQ. ' ') THEN
90    C- e.g., atmosphere : R_low = Top of atmosphere
91    C-            ocean : R_low = Bottom
92         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
93          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
94           DO j=1,sNy           DO j=1,sNy
95            DO i=1,sNx            DO i=1,sNx
96             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.  
97            ENDDO            ENDDO
98           ENDDO           ENDDO
99          ENDDO          ENDDO
100         ENDDO         ENDDO
101        ELSE        ELSE
102         _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
103  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
104  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.
105         CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )          CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
106  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
107  C The 0 is the "iteration" argument. The ' ' is an empty suffix  C The 0 is the "iteration" argument. The ' ' is an empty suffix
108  C      CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )  c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
109  C Read the bathymetry using the low-level I/O package  C Read the bathymetry using the low-level I/O package
110  C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,  c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
111  C    &                    'RS', 1, H, 1, myThid )  c    &                     'RS', 1, R_low, 1, myThid )
112         _END_MASTER(myThid)          _END_MASTER(myThid)
113    
114        ENDIF        ENDIF
115    C- end setup R_low in the interior
116    
117        _EXCH_XY_R4(    H, myThid )  C- fill in the overlap :
118          _EXCH_XY_R4(R_low, myThid )
119    
120  C  c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
121  C     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)  c     _BEGIN_MASTER( myThid )
122  C  c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
123    c     _END_MASTER(myThid)
124    
125    c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
126    
127    C------
128    C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
129    C------
130    
131          IF ( groundAtK1 .AND. bathyFile.NE.' ' ) THEN
132    C------ read directly Po_surf from bathyFile (only for backward compatibility)
133    
134            _BEGIN_MASTER( myThid )
135            CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
136            _END_MASTER(myThid)
137            _BARRIER
138    
139          ELSEIF ( topoFile.EQ.' ' ) THEN
140    C------ set default value:
141    
142            DO bj = myByLo(myThid), myByHi(myThid)
143             DO bi = myBxLo(myThid), myBxHi(myThid)
144              DO j=1,sNy
145               DO i=1,sNx
146                Ro_surf(i,j,bi,bj) = Ro_SeaLevel
147               ENDDO
148              ENDDO
149             ENDDO
150            ENDDO
151    
152          ELSE
153    C------ read from file:
154    
155    C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
156            _BEGIN_MASTER( myThid )
157            CALL READ_REC_XY_RS( topoFile, Hloc, 1, 0, myThid )
158            _END_MASTER(myThid)
159            _BARRIER
160    
161            IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
162    C----
163    C   Convert Surface Geopotential to (reference) Surface Pressure
164    C   according to Tref profile, using same discretisation as in calc_phi_hyd
165    C----
166    c         _BEGIN_MASTER( myThid )
167    c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',Hloc,0,myThid)
168    c         _END_MASTER(myThid)
169    
170              CALL INI_P_GROUND( Hloc, Ro_surf, myThid )
171    
172              _BARRIER
173              _BEGIN_MASTER( myThid )
174              CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
175              _END_MASTER(myThid)
176    
177            ELSE
178    C----
179    C   Direct Transfer to Ro_surf :
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) = Hloc(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.6
  ViewVC Help
Powered by ViewVC 1.1.22